Changeset 2258 for trunk/MagicSoft/Mars

Timestamp:

07/01/03 16:44:15 (21 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mhistmc/MHMcCT1CollectionArea.cc (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -6 +6 @@
      - changed "const TArrayD xed(10,xedge);"  to:
                "const TArrayD xed; xed.Set(10,xedge);"
-       (an the same for yed). Otherwise, a funny Error message was 
+       (and the same for yed). Otherwise, a funny Error message was 
        printed about TArrayD, although everything worked.
 

trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc

@@ -254 +254 @@
     {
       // This theta is not exactly the one of the MC events, just about 
-      // the same:
+      // the same (bins have been selected so):
+
       Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin);
 
-      Float_t emin1, emax1, emin2, emax2;
-      Float_t index, expo, k1, k2;
-      Float_t numevts1, numevts2;
-      Float_t r1, r2;        // Impact parameter range (on ground).
-
-      emin1 = 0;  emax1 = 0; emin2 = 0;  emax2 = 0;
-      expo = 0.; k1 = 0.; k2 = 0.; r1 = 0.; r2 = 0.;
-      numevts1 = 0; numevts2 = 0;
+      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).
+
+      memset(emin,    0, 4*sizeof(Float_t));
+      memset(emax,    0, 4*sizeof(Float_t));
+      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));
+
+      //
+      // rmin and rmax are the minimum and maximum values of the impact 
+      // parameter of the shower on the ground (horizontal plane).
+      //
+
+      Int_t num_MC_samples = 0;
 
       if (theta > 8 && theta < 9)   // 8.75 deg
         {
-          r1 = 0.;
-          r2 = 250.;     //meters
-          emin1 = 300.;
-          emax1 = 400.;  // Energies in GeV.
-          emin2 = 400.;
-          emax2 = 30000.;
-          numevts1 = 4000.;
-          numevts2 = 25740.;
+          rmin[0] = 0.;
+          rmax[0] = 250.;     //meters
+          emin[0] = 300.;
+          emax[0] = 400.;  // Energies in GeV.
+          index[0] = 1.5;
+          numevts[0] = 4000.;
+
+          rmin[1] = 0.;
+          rmax[1] = 250.;     //meters
+          emin[1] = 400.;
+          emax[1] = 30000.;
+          index[1] = 1.5;
+          numevts[1] = 25740.;
+
+          num_MC_samples = 2;
         }
       else if (theta > 20 && theta < 21)  // 20.5 deg 
         {
-          r1 = 0.;
-          r2 = 263.;     //meters
-          emin1 = 300.;
-          emax1 = 400.;  // Energies in GeV.
-          emin2 = 400.;
-          emax2 = 30000.;
-          numevts1 = 6611.;
-          numevts2 = 24448.;
+          rmin[0] = 0.;
+          rmax[0] = 263.;     //meters
+          emin[0] = 300.;
+          emax[0] = 400.;  // Energies in GeV.
+          index[0] = 1.5;
+          numevts[0] = 6611.;
+
+          rmin[1] = 0.;
+          rmax[1] = 263.;     //meters
+          emin[1] = 400.;
+          emax[1] = 30000.;
+          index[1] = 1.5;
+          numevts[1] = 24448.;
+
+          num_MC_samples = 2;
         }
       else if (theta > 26 && theta < 27)  // 26.5 degrees
         {
-          r1 = 0.;
-          r2 = 290.;     //meters
-          emin1 = 300.;
-          emax1 = 400.;  // Energies in GeV.
-          emax2 = emax1;          emin2 = 400.;
-          emax2 = 30000.;
-          numevts1 = 4000.;
-          numevts2 = 26316.;
+          rmin[0] = 0.;
+          rmax[0] = 290.;     //meters
+          emin[0] = 300.;
+          emax[0] = 400.;  // Energies in GeV.
+          index[0] = 1.5;
+          numevts[0] = 4000.;
+
+          rmin[1] = 0.;
+          rmax[1] = 290.;     //meters
+          emin[1] = 400.;
+          emax[1] = 30000.;
+          index[1] = 1.5;
+          numevts[1] = 26316.;
+
+          num_MC_samples = 2;
         }
       else if (theta > 32 && theta < 33)  // 32.5 degrees
         {
-          r1 = 0.;
-          r2 = 350.;     //meters
-          emin1 = 300.;
-          emax1 = 30000.;  // Energies in GeV.
-          emax2 = emax1;
-          numevts1 = 33646.;
+          rmin[0] = 0.;
+          rmax[0] = 350.;     //meters
+          emin[0] = 300.;
+          emax[0] = 30000.;  // Energies in GeV.
+          index[0] = 1.5;
+          numevts[0] = 33646.;
+
+          num_MC_samples = 1;
         }
       else if (theta > 38 && theta < 39)  // 38.75 degrees
         {
-          r1 = 0.;
-          r2 = 380.;     //meters
-          emin1 = 300.;
-          emax1 = 30000.;  // Energies in GeV.
-          emax2 = emax1;
-          numevts1 = 38415.;
+          rmin[0] = 0.;
+          rmax[0] = 380.;     //meters
+          emin[0] = 300.;
+          emax[0] = 30000.;  // Energies in GeV.
+          index[0] = 1.5;
+          numevts[0] = 38415.;
+
+          num_MC_samples = 1;
         }
       else if (theta > 45 && theta < 47)  // 46 degrees
         {
-          r1 = 0.;
-          r2 = 565.;     //meters
-          emin1 = 300.;
-          emax1 = 50000.;  // Energies in GeV.
-          emax2 = emax1;
-          numevts1 = 30197.;
-        }
-
-      index = 1.5;     // Differential spectral Index.
-      expo = 1.-index;
-      k1 = numevts1 / (pow(emax1,expo) - pow(emin1,expo));
-      k2 = numevts2 / (pow(emax2,expo) - pow(emin2,expo));
+          rmin[0] = 0.;
+          rmax[0] = 565.;     //meters
+          emin[0] = 300.;
+          emax[0] = 50000.;  // Energies in GeV.
+          index[0] = 1.5;
+          numevts[0] = 30197.;
+
+          num_MC_samples = 1;
+        }
+
+      //
+      // TO BE FIXED: Add the cases for 55 and 65 degrees zenith angle.
+      //
 
       for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++)
         {
-          const Float_t e1 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i));
-          const Float_t e2 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i+1));
-
-          if (e1 < emin1 || e2 > emax2)
-            continue;
-
-          Float_t events;
-
-          if (e2 <= emax1)
-            events = k1 * (pow(e2, expo) - pow(e1, expo));
-          else if (e1 >= emin2)
-            events = k2 * (pow(e2, expo) - pow(e1, expo));
-          else
-            events = 
-              k1 * (pow(emax1, expo) - pow(e1, expo))+
-              k2 * (pow(e2, expo) - pow(emin2, expo));
+          Float_t e1 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i));
+          Float_t e2 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i+1));
+
+          Float_t events = 0.;
+
+          for (Int_t sample = 0; sample < num_MC_samples; sample++)
+            {
+              Float_t expo = 1.-index[sample];
+              Float_t k = numevts[sample] / (pow(emax[sample],expo) - pow(emin[sample],expo));
+
+              if (e2 < emin[sample] || e1 > emax[sample])
+                continue;
+
+              if (emin[sample] > e1) 
+                e1 = emin[sample];
+
+              if (emax[sample] < e2) 
+                e2 = emax[sample];
+
+              events += k * (pow(e2, expo) - pow(e1, expo));
+            }
 
           fHistAll->SetBinContent(i, thetabin, events);
@@ -357 +401 @@
       // -----------------------------------------------------------
 
-      const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);
+      const Float_t dr = TMath::Pi() * (rmax[0]*rmax[0] - rmin[0]*rmin[0]);
 
       for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
@@ -392 +436 @@
           const Double_t efferr = sqrt((1.-eff)*Ns)/Na;
 
-
+          //
+          // 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.);