Changeset 2250

Timestamp:

06/27/03 17:29:49 (22 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mhistmc/MHMcCollectionArea.cc (modified) (4 diffs)
• mhistmc/MHMcCollectionArea.h (modified) (1 diff)
• mmontecarlo/MMcCollectionAreaCalc.cc (modified) (5 diffs)
• mmontecarlo/MMcCollectionAreaCalc.h (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+
+ 2003/06/27: Abelardo Moralejo
+
+   * mmontecarlo/MMcCollectionAreaCalc.[h,cc]:
+   * mhistmc/MHMcCollectionAreaCalc.[h,cc]:
+
+     - Adapted to allow their use with multiple files containing
+       MC data generated with diffferent energy spectra, even with
+       camera files which have only triggered events inside. Now the
+       histogram containing all showers (before trigger) is filled
+       in the ReInit function, and calculation of collection area
+       is done by CalcEfficiency2(). Some simplifications and cleaning 
+       are still possible.
 
  2003/06/27: Thomas Bretz

trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc

@@ -237 +237 @@
 //  flag
 //
-void MHMcCollectionArea::CalcEfficiency(UInt_t numevts, Float_t emin, Float_t emax, Float_t index)
+void MHMcCollectionArea::CalcEfficiency2()
 {
     // Here we estimate the total number of showers in each energy bin
@@ -244 +244 @@
 
     TH1D &histsel = *fHistSel->ProjectionX();
-
-    TH1D histall;
+    TH1D &histall = *fHistAll->ProjectionX();
 
     TAxis &xaxis = *histsel.GetXaxis();
-
-    MH::SetBinning(&histall, &xaxis);
     MH::SetBinning(fHistCol, &xaxis);
-
-    const Float_t expo = 1+index;
-
-    const Float_t k = (Float_t)numevts / (pow(emax,expo) - pow(emin,expo));
-
     const Int_t nbinx = xaxis.GetNbins();
-
-    for (Int_t i=1; i<=nbinx; i++)
-    {
-        const Float_t e1 = histall.GetBinLowEdge(i);
-        const Float_t e2 = histall.GetBinLowEdge(i+1);
-
-        if (e1 < emin || e2 > emax)
-            continue;
-
-        const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
-
-        histall.SetBinContent(i, events);
-        histall.SetBinError(i, sqrt(events));
-
-    }
 
     // -----------------------------------------------------------
@@ -282 +259 @@
 
     *fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of 300 meters for the MC events. Check that this is the true one!" <<endl<<endl;
-    const Float_t area = TMath::Pi() * (r2*r2 - r1*r1);
+    const Float_t total_area = TMath::Pi() * (r2*r2 - r1*r1);
 
     for (Int_t ix=1; ix<=nbinx; ix++)
@@ -303 +280 @@
 
         const Double_t efferr = sqrt((1.-eff)*Ns)/Na;
-
-        fHistCol->SetBinContent(ix, eff*area);
-        fHistCol->SetBinError(ix, efferr*area);
+        
+        const Float_t col_area  =  eff * total_area;
+        const Float_t col_area_error =  efferr * total_area;
+
+        fHistCol->SetBinContent(ix, col_area);
+        fHistCol->SetBinError(ix, col_area_error);
     }
 

trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h

@@ -35 +35 @@
     const TH1D *GetHist() const { return fHistCol; }
 
+    TH2D *GetHistAll()    { return fHistAll; }
+
     void Draw(Option_t *option="");
 
     void CalcEfficiency();
-    void CalcEfficiency(UInt_t allevts, Float_t emin, Float_t emax, Float_t index);
+    void CalcEfficiency2();
 
     void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall);

trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc

@@ -63 +63 @@
     AddToBranchList("MMcEvt.fEnergy");
     AddToBranchList("MMcEvt.fImpact");
+
+    fCorsikaVersion           =  0;
+    fAllEvtsTriggered         =  kFALSE;
+    fTotalNumSimulatedShowers =  0;
+    fTheta                    = -1.;
+
 } 
 
@@ -76 +82 @@
     }
 
-    fCollArea = (MHMcCollectionArea*)pList->FindCreateObj("MHMcCollectionArea");
-    if (!fCollArea)
-        return kFALSE;
-
-    fTheta                    =  -1;
-    fEmin                     =  -1;
-    fEmax                     =  -1;
-    fSlope                    =   0;
-    fTotalNumSimulatedShowers =   0;
-    fCorsikaVersion           =   0;
-    fAllEvtsTriggered         = kFALSE;
-
     return kTRUE;
 }
@@ -108 +102 @@
 
     fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
-
     *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;
+
 
     if (fTheta>=0 && fTheta!=runheader->GetTelesTheta())
@@ -116 +110 @@
         *fLog << fTheta << ", " << runheader->GetTelesTheta() << ")..." << endl;
     }
-
     fTheta = runheader->GetTelesTheta();
+
 
     if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
         *fLog << warn << dbginf << "Warning - Read files have different Corsika versions..." << endl;
-
     fCorsikaVersion = runheader->GetCorsikaVersion();
 
-    if ( fEmin > 0 &&
-         (fEmin  != corrunheader->GetELowLim() ||
-          fEmax  != corrunheader->GetEUppLim() ||
-          fSlope != corrunheader->GetSlopeSpec()) )
-      *fLog << warn << dbginf << "Warning - Read files have different energy distribution..." << endl;
-
-    fEmin  = corrunheader->GetELowLim();
-    fEmax  = corrunheader->GetEUppLim();
-    fSlope = corrunheader->GetSlopeSpec();
 
     fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
-
     *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;
 
+
+    fCollArea = (MHMcCollectionArea*)plist->FindCreateObj("MHMcCollectionArea");
+    if (!fCollArea)
+        return kFALSE;
+
     if (fAllEvtsTriggered)
+      {
+        //
+        // Calculate approximately the original number of events in each 
+        // energy bin:
+        //
+
+        const Float_t emin = corrunheader->GetELowLim();
+        const Float_t emax = corrunheader->GetEUppLim();
+        const Float_t expo = 1 + corrunheader->GetSlopeSpec();
+        const Float_t k = runheader->GetNumSimulatedShowers() / 
+          (pow(emax,expo) - pow(emin,expo)) ;
+
+        const Int_t nbinx = fCollArea->GetHistAll()->GetNbinsX();
+
+        for (Int_t i = 1; i <= nbinx; i++)
+          {
+            const Float_t e1 = fCollArea->GetHistAll()->GetXaxis()->GetBinLowEdge(i);
+            const Float_t e2 = fCollArea->GetHistAll()->GetXaxis()->GetBinLowEdge(i+1);
+
+            if (e1 < emin || e2 > emax)
+              continue;
+
+            const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
+            //
+            // We fill the ith energy bin, with the total number of events
+            // Second argument of Fill would be impact parameter of each
+            // event, but we don't really need it for the collection area,
+            // so we just put a dummy value (1.)
+            //
+
+            const Float_t energy = (e1+e2)/2.;
+            fCollArea->GetHistAll()->Fill(energy, 1., events);
+          }
+
         return kTRUE;
+      }
 
     fMcTrig = (MMcTrig*)plist->FindObject(fObjName);
@@ -190 +213 @@
     {
         *fLog << inf << "Total number of showers: " << fTotalNumSimulatedShowers << endl;
-        fCollArea->CalcEfficiency(fTotalNumSimulatedShowers, fEmin, fEmax, fSlope);
-    }
-
-    return kTRUE;
-}
-
+        fCollArea->CalcEfficiency2();
+    }
+
+    return kTRUE;
+}
+

trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h

@@ -5 +5 @@
 #include "MTask.h"
 #endif
+
+#include <TH2.h>
 
 class MParList;
@@ -23 +25 @@
     UInt_t fTotalNumSimulatedShowers;
     UInt_t fCorsikaVersion;
+    Bool_t fAllEvtsTriggered;
     Float_t fTheta;
-    Float_t fEmin;
-    Float_t fEmax;
-    Float_t fSlope;
-
-    Bool_t fAllEvtsTriggered;
 
     Bool_t ReInit(MParList *plist);