Changeset 1823 for trunk/MagicSoft/Mars

Timestamp:

03/13/03 22:25:00 (22 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• macros/CT1collarea.C (added)
• mhist/MHMcCT1CollectionArea.cc (modified) (8 diffs)
• mhist/MHMcCT1CollectionArea.h (modified) (3 diffs)
• mmontecarlo/MMcCT1CollectionAreaCalc.cc (modified) (5 diffs)
• mmontecarlo/MMcCT1CollectionAreaCalc.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -4 +4 @@
 
     * mhist/MHMcCT1CollectionArea.[h,cc]
-      - Added for calculations of collection area for CT1.
+      - Added for calculations of collection area for CT1.Contains three 2-d 
+        histograms with axis energy vs theta angle: one histogram for all 
+        events, one for analyzed events, one for the collection area.
 
     * mmontecarlo/MMcCT1CollectionAreaCalc.[h,cc]
-      - Added for the same reason.
+      - Added for the same reason. 
 
     * macros/CT1collarea.C

trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc

@@ -52 +52 @@
     //   we set the energy range from 10 Gev to 30000 GeV (in log 4.5 orders
     //   of magnitude) and for each order we take 10 subdivision --> 45 xbins
-    //
-    //   we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins
+    //   we set the theta range from 12.5 to 48 deg, with 6 bins.
     //
     fName  = name  ? name  : "MHMcCT1CollectionArea";
@@ -60 +59 @@
     MBinning binsx;
     MBinning binsy;
+
     binsx.SetEdgesLog(45, 10, 30000);
-    binsy.SetEdges   (50,  0,   500);
+    const Double_t yedge[7] = {12.5, 17.5, 23.5, 29.5, 35.5, 42., 48.};
+    const TArrayD yed(7,yedge);
+    binsy.SetEdges(yed);
 
     fHistAll = new TH2D;
     fHistSel = new TH2D;
-    fHistCol = new TH1D;
- 
+    fHistCol = new TH2D;
+
     MH::SetBinning(fHistAll, &binsx, &binsy);
     MH::SetBinning(fHistSel, &binsx, &binsy);
+    MH::SetBinning(fHistCol, &binsx, &binsy);
+
 
     fHistCol->SetName(fName);
@@ -75 +79 @@
 
     fHistCol->SetTitle(fTitle);
-    fHistAll->SetTitle("All showers - Radius vs Energy distribution");
-    fHistSel->SetTitle("Selected showers - Radius vs Energy distribution");
+    fHistAll->SetTitle("All showers - Theta vs Energy distribution");
+    fHistSel->SetTitle("Selected showers - Theta vs Energy distribution");
 
     fHistAll->SetDirectory(NULL);
@@ -83 +87 @@
 
     fHistAll->SetXTitle("E [GeV]");
-    fHistAll->SetYTitle("r [m]");
+    fHistAll->SetYTitle("theta [deg]");
     fHistAll->SetZTitle("N");
 
     fHistSel->SetXTitle("E [GeV]");
-    fHistSel->SetYTitle("r [m]");
+    fHistSel->SetYTitle("theta [deg]");
     fHistSel->SetYTitle("N");
 
     fHistCol->SetXTitle("E [GeV]");
-    fHistCol->SetYTitle("A [m^{2}]");
+    fHistCol->SetYTitle("theta [deg]");
+    fHistCol->SetZTitle("A [m^{2}]");
 }
 
@@ -109 +114 @@
 // Fill data into the histogram which contains the selected showers
 //
-void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t radius)
-{
-    fHistSel->Fill(energy, radius);
+void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t theta)
+{
+    fHistSel->Fill(energy, theta);
 }
 
@@ -159 +164 @@
 
     //
-    // This is necessary to get the expected bahviour of DrawClone
+    // This is necessary to get the expected behaviour of DrawClone
     //
     gROOT->SetSelectedPad(NULL);
@@ -190 +195 @@
 //  and set the 'ReadyToSave' flag
 //
-void MHMcCT1CollectionArea::CalcEfficiency(Float_t theta)
+void MHMcCT1CollectionArea::CalcEfficiency()
 {
   //
@@ -205 +210 @@
   //
 
-    TH1D &histsel = *fHistSel->ProjectionX();
-
-    TH1D histall;
-
-    TAxis &xaxis = *histsel.GetXaxis();
-
-    MH::SetBinning(&histall, &xaxis);
-    MH::SetBinning(fHistCol, &xaxis);
-
-    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;
-
-    if (theta > 0.26 && theta < 0.27)  // 15 degrees
-      {
-        r1 = 0.;
-        r2 = 250.;     //meters
-        emin1 = 300.;
-        emax1 = 400.;  // Energies in GeV.
-        emin2 = 400.;
-        emax2 = 30000.;
-        numevts1 = 4000.;
-        numevts2 = 25740.;
-      }
-    else if (theta > 0.34 && theta < 0.35)  // 20 degrees
-      {
-        r1 = 0.;
-        r2 = 263.;     //meters
-        emin1 = 300.;
-        emax1 = 400.;  // Energies in GeV.
-        emin2 = 400.;
-        emax2 = 30000.;
-        numevts1 = 6611.;
-        numevts2 = 24448.;
-      }
-    else if (theta > 0.47 && theta < 0.48)  // 27 degrees
-      {
-        r1 = 0.;
-        r2 = 290.;     //meters
-        emin1 = 300.;
-        emax1 = 400.;  // Energies in GeV.
-        emin2 = 400.;
-        emax2 = 30000.;
-        numevts1 = 4000.;
-        numevts2 = 26316.;
-      }
-    else if (theta > 0.55 && theta < 0.56)  // 32 degrees
-      {
-        r1 = 0.;
-        r2 = 350.;     //meters
-        emin1 = 300.;
-        emax1 = 30000.;  // Energies in GeV.
-        numevts1 = 33646.;
-      }
-    else if (theta > 0.68 && theta < 0.69)  // 39 degrees
-      {
-        r1 = 0.;
-        r2 = 380.;     //meters
-        emin1 = 300.;
-        emax1 = 30000.;  // Energies in GeV.
-        numevts1 = 38415.;
-      }
-    else if (theta > 0.78 && theta < 0.79)  // 45 degrees
-      {
-        r1 = 0.;
-        r2 = 565.;     //meters
-        emin1 = 300.;
-        emax1 = 50000.;  // Energies in GeV.
-        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));
-
-    const Int_t nbinx = xaxis.GetNbins();
-
-    for (Int_t i=1; i<=nbinx; i++)
+  for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++)
     {
-        const Float_t e1 = histall.GetBinLowEdge(i);
-        const Float_t e2 = histall.GetBinLowEdge(i+1);
-
-        if (e1 < emin1 || e2 > emax2)
+      // This theta is not exactly the one of the MC events, just about 
+      // the same:
+      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;
+
+      if (theta > 14 && theta < 16)   // 15 deg
+        {
+          r1 = 0.;
+          r2 = 250.;     //meters
+          emin1 = 300.;
+          emax1 = 400.;  // Energies in GeV.
+          emin2 = 400.;
+          emax2 = 30000.;
+          numevts1 = 4000.;
+          numevts2 = 25740.;
+        }
+      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.;
+        }
+      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.;
+        }
+      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.;
+        }
+      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.;
+        }
+      else if (theta > 44 && theta < 46)  // 45 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));
+
+      for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++)
+        {
+          const Float_t e1 = fHistAll->GetXaxis()->GetBinLowEdge(i);
+          const Float_t e2 = 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));
-
-        histall.SetBinContent(i, events);
-        histall.SetBinError(i, sqrt(events));
+          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));
+
+          fHistAll->SetBinContent(i, thetabin, events);
+          fHistAll->SetBinError(i, thetabin, sqrt(events));
+        }
+
+      // -----------------------------------------------------------
+
+      const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);
+
+      for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
+        {
+          const Float_t Na = fHistAll->GetBinContent(ix,thetabin);
+
+          if (Na <= 0)
+            continue;
+
+          const Float_t Ns = fHistSel->GetBinContent(ix,thetabin);
+
+          // Since Na is an estimate of the total number of showers generated
+          // in the energy bin, it may happen that Ns (triggered showers) is
+          // larger than Na. In that case, the bin is skipped:
+
+          if (Na < Ns)
+            continue;
+
+          const Double_t eff = Ns/Na;
+          const Double_t err = sqrt((1.-eff)*Ns)/Na;
+
+          const Float_t area = dr * cos(theta*TMath::Pi()/180.);
+
+          fHistCol->SetBinContent(ix, thetabin, eff*area);
+          fHistCol->SetBinError(ix, thetabin, err*area);
+        }
     }
 
-    // -----------------------------------------------------------
-
-    const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);
-
-    for (Int_t ix=1; ix<=nbinx; ix++)
-      {
-        const Float_t Na = histall.GetBinContent(ix);
-
-        if (Na <= 0)
-          continue;
-
-        const Float_t Ns = histsel.GetBinContent(ix);
-
-        // Since Na is an estimate of the total number of showers generated
-        // in the energy bin, it may happen that Ns (triggered showers) is
-        // larger than Na. In that case, the bin is skipped:
-
-        if (Na < Ns)
-          continue;
-
-        const Double_t eff = Ns/Na;
-
-        const Double_t err = sqrt((1.-eff)*Ns)/Na;
-
-        const Float_t area = dr * cos(theta);
-
-        fHistCol->SetBinContent(ix, eff*area);
-        fHistCol->SetBinError(ix, err*area);
-      }
-
-    delete &histsel;
-
     SetReadyToSave();
 }
+

trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h

@@ -15 +15 @@
     TH2D *fHistAll; //! all simulated showers
     TH2D *fHistSel; //! the selected showers
-
-    TH1D *fHistCol; //  the collection area
+    TH2D *fHistCol; //  the collection area
 
 public:
@@ -27 +26 @@
     void DrawSel(Option_t *option="");
 
-    const TH1D *GetHist() const { return fHistCol; }
+    const TH2D *GetHist() const { return fHistCol; }
 
     void Draw(Option_t *option="");
     TObject *DrawClone(Option_t *option="") const;
 
-    void CalcEfficiency(Float_t theta);
+    void CalcEfficiency();
 
     ClassDef(MHMcCT1CollectionArea, 1)  // Data Container to calculate Collection Area
@@ -38 +37 @@
 
 #endif
+
+

trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.cc

@@ -51 +51 @@
 
     AddToBranchList("MMcEvt.fEnergy");
-    AddToBranchList("MMcEvt.fImpact");
     AddToBranchList("MMcEvt.fTelescopeTheta");
 } 
@@ -72 +71 @@
     fTotalNumSimulatedShowers =  0;
     fAllEvtsTriggered         = kTRUE;
-    fAllEvtsAtSameTheta       = kTRUE;
-    fTelescopeTheta           = -1.;
 
     return kTRUE;
@@ -81 +78 @@
 {
     const Double_t energy = fMcEvt->GetEnergy();
-    const Double_t impact = fMcEvt->GetImpact()/100.;
 
-    if (fTelescopeTheta < 0.)
-      fTelescopeTheta = fMcEvt->GetTelescopeTheta();
+    Double_t TelescopeTheta = 180.*fMcEvt->GetTelescopeTheta()/TMath::Pi();
 
-    if (fTelescopeTheta != fMcEvt->GetTelescopeTheta())
-      {
-        fAllEvtsAtSameTheta = kFALSE;
-      }
-
-    fCollArea->FillSel(energy, impact);
+    fCollArea->FillSel(energy, TelescopeTheta);
 
     return kTRUE;
@@ -98 +88 @@
 Bool_t MMcCT1CollectionAreaCalc::PostProcess()
 { 
-  if ( ! fAllEvtsAtSameTheta )
-    {
-        *fLog << err << dbginf << "ERROR: input data contain events at different zenith angles (not supported)... exiting." << endl;
-        return kFALSE;
-    }
   //
   //   do the calculation of the effective area
@@ -108 +93 @@
   *fLog << inf << "Calculation Collection Area..." << endl;
   
-  fCollArea->CalcEfficiency(fTelescopeTheta);
+  fCollArea->CalcEfficiency();
 
   return kTRUE;

trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.h

@@ -21 +21 @@
     UInt_t fTotalNumSimulatedShowers;
     Bool_t fAllEvtsTriggered;
-    Bool_t fAllEvtsAtSameTheta;
-    Float_t fTelescopeTheta;
 
 public: