Changeset 2261 for trunk/MagicSoft

Timestamp:

07/03/03 15:09:11 (21 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• macros/CT1collarea.C (modified) (3 diffs)
• mhistmc/MHMcCT1CollectionArea.cc (modified) (18 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -4 +4 @@
  2003/07/03: Abelardo Moralejo
 
-  * mhistmc/MHMcCT1CollectionArea.cc
-    - Added code to allow the calculation of CT1 collection areas 
-      at 55 and 65 degrees (from the events in DK's MC library)
+   * mhistmc/MHMcCT1CollectionArea.cc
+     - Added code to allow the calculation of CT1 collection areas 
+       at 55 and 65 degrees (from the events in DK's MC library)
+
+   * macros/CT1collarea.C
+     - Changed binning in theta to include high ZAs
 
 

trunk/MagicSoft/Mars/macros/CT1collarea.C

@@ -25 +25 @@
 
 
-void CT1collarea(TString filename="MC1.root", TString outname="")
+void CT1collarea(TString filename="MC_SC4.root", TString outname="")
 { 
     //
@@ -46 +46 @@
 
     MBinning binsx("BinningE");
+
     //    binsx.SetEdges(30,2.,5.);
+    //    Double_t xedge[10] =
+    //    {2.50515, 2.69897, 2.89763, 3.09691, 3.30103, 3.49969, 3.69984, 3.89982, 4.10003, 4.29994};
 
-    Double_t xedge[10] =
-    {2.50515, 2.69897, 2.89763, 3.09691, 3.30103, 3.49969, 3.69984, 3.89982, 4.10003, 4.29994};
+
+    Double_t xedge[15] = {2.47712, 2.64345, 2.82607, 3., 3.17609, 3.35218, 3.52892, 3.70415, 3.88024, 4.05652, 4.23274, 4.40875, 4.58478, 4.76080, 4.90309};
+
     const TArrayD xed;
-    xed.Set(10,xedge);
+    xed.Set(15,xedge);
     binsx.SetEdges(xed);
 
     MBinning binsy("BinningTheta");
-    const Double_t yedge[7] = {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50.};
+    const Double_t yedge[9] = {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
     const TArrayD yed;
-    yed.Set(7,yedge);
+    yed.Set(9,yedge);
     binsy.SetEdges(yed);
 
@@ -66 +70 @@
     tasklist.AddToList(&reader);   
 
-    /*
-    MF filterhadrons("HadNN.fHadronness<0.25");
+    MF filterhadrons("HadSC.fHadronness<0.3 && {abs(MHillasSrc.fAlpha)} < 13.1");
     tasklist.AddToList(&filterhadrons);
-    */
+
 
     MFillH filler("MHMcCT1CollectionArea","MMcEvt");
-    //    filler.SetFilter(&filterhadrons);
+    filler.SetFilter(&filterhadrons);
     tasklist.AddToList(&filler);
 

trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc

@@ -250 +250 @@
   //
   //
+  //
+  // Only for the binning taken from D. Kranich :
+  //
 
   for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++)
@@ -258 +261 @@
       Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin);
 
-      Float_t emin[4];
-      Float_t emax[4];
-      Float_t index[4];
-      Float_t numevts[4];
-      Float_t rmin[4];
-      Float_t rmax[4];     // Impact parameter range (on ground).
+      Float_t emin[4];       // Minimum energy in MC sample
+      Float_t emax[4];       // Maximum energy in MC sample
+      Float_t index[4];      // Spectral index
+      Float_t numevts[4];    // Number of events
+      Float_t multfactor[4]; // Factor by which the original number of events in an MC 
+                             // sample has been multiplied to account for the differences
+                             // in the generation areas of the various samples.
+      Float_t rmax;          // Maximum impact parameter range (on ground up to 45 degrees, 
+                             // on a plane perpendicular to Shower axis for 55 and 65 deg).
 
       memset(emin,    0, 4*sizeof(Float_t));
@@ -269 +275 @@
       memset(index,   0, 4*sizeof(Float_t));
       memset(numevts, 0, 4*sizeof(Float_t));
-      memset(rmin,    0, 4*sizeof(Float_t));
-      memset(rmax,    0, 4*sizeof(Float_t));
+      rmax = 0.;
+
+      multfactor[0] = 1.;
+      multfactor[1] = 1.;
+      multfactor[2] = 1.;
+      multfactor[3] = 1.;
 
       //
@@ -281 +291 @@
       if (theta > 8 && theta < 9)   // 8.75 deg
         {
-          rmin[0] = 0.;
-          rmax[0] = 250.;     //meters
           emin[0] = 300.;
           emax[0] = 400.;  // Energies in GeV.
@@ -288 +296 @@
           numevts[0] = 4000.;
 
-          rmin[1] = 0.;
-          rmax[1] = 250.;     //meters
           emin[1] = 400.;
           emax[1] = 30000.;
@@ -295 +301 @@
           numevts[1] = 25740.;
 
+          rmax = 250.;     //meters
           num_MC_samples = 2;
         }
       else if (theta > 20 && theta < 21)  // 20.5 deg 
         {
-          rmin[0] = 0.;
-          rmax[0] = 263.;     //meters
           emin[0] = 300.;
           emax[0] = 400.;  // Energies in GeV.
@@ -306 +311 @@
           numevts[0] = 6611.;
 
-          rmin[1] = 0.;
-          rmax[1] = 263.;     //meters
           emin[1] = 400.;
           emax[1] = 30000.;
@@ -313 +316 @@
           numevts[1] = 24448.;
 
+          rmax = 263.;
           num_MC_samples = 2;
         }
       else if (theta > 26 && theta < 27)  // 26.5 degrees
         {
-          rmin[0] = 0.;
-          rmax[0] = 290.;     //meters
           emin[0] = 300.;
           emax[0] = 400.;  // Energies in GeV.
@@ -324 +326 @@
           numevts[0] = 4000.;
 
-          rmin[1] = 0.;
-          rmax[1] = 290.;     //meters
           emin[1] = 400.;
           emax[1] = 30000.;
@@ -331 +331 @@
           numevts[1] = 26316.;
 
+          rmax = 290.;     //meters
           num_MC_samples = 2;
         }
       else if (theta > 32 && theta < 33)  // 32.5 degrees
         {
-          rmin[0] = 0.;
-          rmax[0] = 350.;     //meters
           emin[0] = 300.;
           emax[0] = 30000.;  // Energies in GeV.
@@ -342 +341 @@
           numevts[0] = 33646.;
 
+          rmax = 350.;     //meters
           num_MC_samples = 1;
         }
       else if (theta > 38 && theta < 39)  // 38.75 degrees
         {
-          rmin[0] = 0.;
-          rmax[0] = 380.;     //meters
           emin[0] = 300.;
           emax[0] = 30000.;  // Energies in GeV.
@@ -353 +351 @@
           numevts[0] = 38415.;
 
+          rmax = 380.;     //meters
           num_MC_samples = 1;
         }
       else if (theta > 45 && theta < 47)  // 46 degrees
         {
-          rmin[0] = 0.;
-          rmax[0] = 565.;     //meters
           emin[0] = 300.;
           emax[0] = 50000.;  // Energies in GeV.
@@ -364 +361 @@
           numevts[0] = 30197.;
 
+          rmax = 565.;     //meters
           num_MC_samples = 1;
         }
-
-      //
-      // TO BE FIXED: Add the cases for 55 and 65 degrees zenith angle.
-      //
+      else if (theta > 54 && theta < 56)  // 55 degrees
+        {
+          //
+          // The value of numevts in the first sample (below) has been
+          // changed to simplify calculations. We have multiplied it
+          // times 1.2808997 to convert it to the number it would be if 
+          // the generation area was equal to that of the other samples 
+          // at 55 degrees (pi*600**2 m2). This has to be taken into account 
+          // in the error in the number of events.
+          //
+
+          emin[0] = 500.;
+          emax[0] = 50000.;  // Energies in GeV.
+          index[0] = 1.5;
+          numevts[0] = 3298.;
+          multfactor[0] = 1.2808997;
+
+          emin[1] = 1500.;
+          emax[1] = 50000.;  // Energies in GeV.
+          index[1] = 1.5;
+          numevts[1] = 22229.;
+
+          emin[2] = 1500.;
+          emax[2] = 50000.;  // Energies in GeV.
+          index[2] = 1.7;
+          numevts[2] = 7553.;
+
+          rmax = 600;     //meters
+          num_MC_samples = 3;
+        }
+
+      else if (theta > 64 && theta < 66)  // 65 degrees
+        {
+          emin[0] = 2000.;
+          emax[0] = 50000.;  // Energies in GeV.
+          index[0] = 1.5;
+          numevts[0] = 16310.;
+
+          emin[1] = 2000.;
+          emax[1] = 50000.;  // Energies in GeV.
+          index[1] = 1.7;
+          numevts[1] = 3000.;
+
+          //
+          // The value of numevts in the next two samples (below) has been
+          // changed to simplify calculations. We have converted them to the
+          // number it would be if the generation area was equal to that of 
+          // the first two samples at 65 degrees (pi*800**2 m2) (four times 
+          // as many, since the original maximum impact parameter was 400
+          // instead of 800. This is taken into account in the error too.
+          //
+
+          emin[2] = 5000.;
+          emax[2] = 50000.;  // Energies in GeV.
+          index[2] = 1.5;
+          numevts[2] = 56584.;
+          multfactor[2] = 4;
+
+          emin[3] = 5000.;
+          emax[3] = 50000.;  // Energies in GeV.
+          index[3] = 1.7;
+          numevts[3] = 11464;
+          multfactor[3] = 4;
+
+          rmax = 800;     // meters
+          num_MC_samples = 4;
+        }
+
 
       for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++)
@@ -377 +439 @@
 
           Float_t events = 0.;
+          Float_t errevents = 0.;
 
           for (Int_t sample = 0; sample < num_MC_samples; sample++)
@@ -393 +456 @@
 
               events += k * (pow(e2, expo) - pow(e1, expo));
+              errevents += multfactor[sample] * events;
             }
 
+          errevents= sqrt(errevents);
+
           fHistAll->SetBinContent(i, thetabin, events);
-          fHistAll->SetBinError(i, thetabin, sqrt(events));
+          fHistAll->SetBinError(i, thetabin, errevents);
         }
 
       // -----------------------------------------------------------
 
-      const Float_t dr = TMath::Pi() * (rmax[0]*rmax[0] - rmin[0]*rmin[0]);
+      const Float_t dr = TMath::Pi() * rmax * rmax;
 
       for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
@@ -415 +481 @@
               // below 1E4, it means that no events triggered -> eff area = 0
               //
-
+              // NOW DISABLED: because collection area after analysis does not
+              // saturate at high E!
+              //
+
+              /*
               if (fHistSel->GetXaxis()->GetBinLowEdge(ix) > 4.)
                 {
@@ -421 +491 @@
                   fHistCol->SetBinError(ix, thetabin, fHistCol->GetBinError(ix-1, thetabin));
                 }
+              */
               continue;
             }
@@ -437 +508 @@
 
           //
-          // Total area, perpendicular to the observation direction
+          // Now we get the total area, perpendicular to the observation direction
           // in which the events were generated (correct for cos theta):
           //
-          const Float_t area = dr * cos(theta*TMath::Pi()/180.);
+
+          Float_t area;
+
+          if (theta < 50)
+            area = dr * cos(theta*TMath::Pi()/180.);
+          else
+            area = dr;
+
+          // Above 50 degrees MC was generated with Corsika 6.xx, and the cores
+          // were distributed on a circle perpendicular to the observation direction, 
+          // and not on ground, hence the correction for cos(theta) is not necessary.
+          //
+
 
           fHistCol->SetBinContent(ix, thetabin, eff*area);