Changeset 2261 for trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc

Timestamp:

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

Author:

moralejo

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc (modified) (18 diffs)

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);