Changeset 6491 for trunk/MagicSoft/Mars/mmontecarlo

Timestamp:

02/15/05 17:20:40 (20 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mmontecarlo

Files:

• MMcCollectionAreaCalc.cc (modified) (4 diffs)
• MMcCollectionAreaCalc.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc

@@ -17 +17 @@
 !
 !   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
-!   Author(s): Harald Kornmayer 1/2001
+!   Author(s): Harald Kornmayer  1/2001
+!   Author(s): Abelardo Moralejo 2/2005 <mailto:moralejo@pd.infn.it>
 !
 !   Copyright: MAGIC Software Development, 2000-2002
@@ -26 +27 @@
 //////////////////////////////////////////////////////////////////////////////
 //
-//  MHMcCollectionAreaCalc
-//
-//  Remark: The initialization is mainly done in the ReInit function.
-//          Please make sure, that you don't use MReadTree when processing
-//          a file. Use a 'ReInit'-calling task like MReadMarsFile
+//  MMcCollectionAreaCalc
 //
 //////////////////////////////////////////////////////////////////////////////
@@ -42 +39 @@
 
 #include "MMcEvt.hxx"
-#include "MMcTrig.hxx"
+#include "MMcEvtBasic.h"
+
 #include "MMcRunHeader.hxx"
 #include "MMcCorsikaRunHeader.h"
@@ -52 +50 @@
 using namespace std;
 
-MMcCollectionAreaCalc::MMcCollectionAreaCalc(const char *input,
-                                             const char *name, const char *title)
+////////////////////////////////////////////////////////////////////////////////
+//
+//  Constructor
+//
+MMcCollectionAreaCalc::MMcCollectionAreaCalc(const char *input, const char *name, 
+                                             const char *title): 
+  fBinsTheta(0), fBinsEnergy(0), fSpectrum(0)
 {
     fName  = name  ? name  : "MMcCollectionAreaCalc";
     fTitle = title ? title : "Task to calculate the collection area";
-
-    fObjName = input ? input : "MMcTrig";
-    AddToBranchList(Form("%s.fNumFirstLevel", input));
-
-    AddToBranchList("MMcEvt.fEnergy");
-    AddToBranchList("MMcEvt.fImpact");
-
-    fCorsikaVersion           =  0;
-    fAllEvtsTriggered         =  kFALSE;
-    fTotalNumSimulatedShowers =  0;
-    fThetaMin                 = -1.;
-    fThetaMax                 = -1.;
-
 } 
 
+////////////////////////////////////////////////////////////////////////////////
+//
+// PreProcess. 
+// Create MHMcCollectionArea object if necessary. 
+// Search in parameter list for MBinning objects with binning in theta and E.
+// These contain the coarse binning to be used in the analysis. Then search 
+// for other necessary input containers: 
+// if MMcEvt is found, it means we are in the loop over the Events tree, 
+// and so we must fill the histogram MHMcCollectionArea::fHistSel (using 
+// MHMcCollectionArea::FillSel). If MMcEvt is not found, it means that we are in 
+// the loop over the "OriginalMC" tree, containing all the original Corsika events 
+// produced, and hence we must fill the histogram  fHistAll through 
+// MHMcCollectionArea::FillAll.
+//
 Int_t MMcCollectionAreaCalc::PreProcess (MParList *pList)
 {
-    // connect the raw data with this task
+  fCollArea = (MHMcCollectionArea*)pList->FindCreateObj("MHMcCollectionArea");
 
-    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
-    if (!fMcEvt)
+  // Look for the binnings of the histograms if they have not been already
+  // found in a previous loop.
+
+  if (!fBinsTheta)
     {
-        *fLog << err << "MMcEvt not found... abort." << endl;
-        return kFALSE;
+      fBinsTheta = (MBinning*) pList->FindObject("binsTheta");
+      if (!fBinsTheta)
+        {
+          *fLog << err << "Coarse Theta binning not found... Aborting." 
+            << endl;
+          return kFALSE;
+        }
+    }
+
+  if (!fBinsEnergy)
+    {
+      fBinsEnergy = (MBinning*) pList->FindObject("binsEnergy");
+      if (!fBinsEnergy)
+        {
+          *fLog << err << "Coarse Energy binning not found... Aborting." 
+                << endl;
+          return kFALSE;
+        }
+    }
+
+  fCollArea->SetCoarseBinnings(*fBinsEnergy, *fBinsTheta);
+
+
+  // Look for the input containers
+
+  fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
+  if (fMcEvt)
+    {
+        *fLog << inf << "MMcEvt found... I will fill MHMcCollectionArea.fHistSel..." << endl;
+        return kTRUE;
+    }
+
+  *fLog << inf << "MMcEvt not found... looking for MMcEvtBasic..." << endl;
+
+
+  fMcEvtBasic = (MMcEvtBasic*) pList->FindObject("MMcEvtBasic");
+    
+  if (fMcEvtBasic)
+    {
+      *fLog << inf << "MMcEvtBasic found... I will fill MHMcCollectionArea.fHistAll..." << endl;
+      return kTRUE;
+    }
+
+  *fLog << err << "MMcEvtBasic not found. Aborting..." << endl;
+
+  return kFALSE;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// Process. Depending on whether fMcEvt or fMcEvtBasic are available, fill 
+// MHMcCollectionArea::fHistAll, else, if fMcEvt is available, fill 
+// MHMcCollectionArea::fHistSel
+//
+Int_t MMcCollectionAreaCalc::Process()
+{
+    Double_t energy = fMcEvt? fMcEvt->GetEnergy() : fMcEvtBasic->GetEnergy();
+    Double_t impact = fMcEvt? fMcEvt->GetImpact()/100. : fMcEvtBasic->GetImpact()/100.; // in m
+    Double_t theta  = fMcEvt? fMcEvt->GetTelescopeTheta()*TMath::RadToDeg() : 
+      fMcEvtBasic->GetTelescopeTheta()*TMath::RadToDeg(); // in deg
+
+    if (fMcEvt)
+      fCollArea->FillSel(energy, impact, theta);
+    else
+      fCollArea->FillAll(energy, impact, theta);
+
+    return kTRUE;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// Postprocess. If both fHistAll and fHistSel are already filled, calculate
+// effective areas. Else it means we still have to run one more loop.
+//
+Int_t MMcCollectionAreaCalc::PostProcess()
+{
+  if ( ((TH2D*)fCollArea->GetHistAll())->GetEntries() > 0 &&
+       ((TH2D*)fCollArea->GetHistSel())->GetEntries() > 0)
+    {
+      *fLog << inf << "Calculation Collection Area..." << endl;
+      fCollArea->Calc(fSpectrum);
     }
 
     return kTRUE;
 }
-
-Bool_t MMcCollectionAreaCalc::ReInit(MParList *plist)
-{
-    fCollArea=0;
-
-    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
-    if (!runheader)
-    {
-        *fLog << err << "Error - MMcRunHeader not found... abort." << endl;
-        return kFALSE;
-    }
-
-    fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
-    *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;
-
-
-    if ( (fThetaMin >= 0 && fThetaMin != runheader->GetShowerThetaMin()) ||
-         (fThetaMax >= 0 && fThetaMax != runheader->GetShowerThetaMax()) )
-    {
-        *fLog << warn << "Warning - Read files have different Theta ranges (";
-        *fLog << "(" << fThetaMin << ", " << fThetaMax << ")   vs   (" << 
-          runheader->GetShowerThetaMin() << ", " << 
-          runheader->GetShowerThetaMax() << ")..." << endl;
-    }
-    fThetaMin = runheader->GetShowerThetaMin();
-    fThetaMax = runheader->GetShowerThetaMax();
-
-
-    if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
-        *fLog << warn << "Warning - Read files have different Corsika versions..." << endl;
-    fCorsikaVersion = runheader->GetCorsikaVersion();
-
-
-    fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
-    *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;
-
-
-    fCollArea = (MHMcCollectionArea*)plist->FindCreateObj("MHMcCollectionArea");
-    if (!fCollArea)
-        return kFALSE;
-
-    if (!fAllEvtsTriggered)
-    {
-        fMcTrig = (MMcTrig*)plist->FindObject(fObjName);
-        if (!fMcTrig)
-        {
-            *fLog << err << fObjName << " [MMcTrig] not found... anort." << endl;
-            return kFALSE;
-        }
-        return kTRUE;
-    }
-
-    MMcCorsikaRunHeader *corrunheader  = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
-    if (!corrunheader)
-    {
-        *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl;
-        return kFALSE;
-    }
-
-    //
-    // 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)) ;
-
-    TH2 &hall = *fCollArea->GetHistAll();
-
-    const Int_t nbinx = hall.GetNbinsX();
-
-    TAxis &axe = *hall.GetXaxis();
-    for (Int_t i = 1; i <= nbinx; i++)
-    {
-        const Float_t e1 = axe.GetBinLowEdge(i);
-        const Float_t e2 = axe.GetBinLowEdge(i+1);
-
-        if (e1 < emin || e2 > emax)
-            continue;
-
-        const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
-        //
-        // We fill the i-th 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.;
-        hall.Fill(energy, 1., events);
-    }
-
-    return kTRUE;
-}
-
-Int_t MMcCollectionAreaCalc::Process()
-{
-    const Double_t energy = fMcEvt->GetEnergy();
-    const Double_t impact = fMcEvt->GetImpact()/100.;
-
-    //
-    // This happens for camera files created with Camera 0.5
-    //
-    if (TMath::IsNaN(impact))
-    {
-        *fLog << err << dbginf << "ERROR - Impact=NaN (Not a number)... abort." << endl;
-        return kERROR;
-    }
-
-    if (!fAllEvtsTriggered)
-    {
-        fCollArea->FillAll(energy, impact);
-
-        if (fMcTrig->GetFirstLevel() <= 0)
-            return kTRUE;
-    }
-
-    fCollArea->FillSel(energy, impact);
-
-    return kTRUE;
-}
-
-Int_t MMcCollectionAreaCalc::PostProcess()
-{
-    if (!fCollArea)
-        return kTRUE;
-
-    //
-    //   do the calculation of the effectiv area
-    //
-    *fLog << inf << "Calculation Collection Area..." << endl;
-
-    if (!fAllEvtsTriggered)
-        fCollArea->CalcEfficiency();
-    else
-    {
-        *fLog << inf << "Total number of showers: " << fTotalNumSimulatedShowers << endl;
-        fCollArea->CalcEfficiency2();
-    }
-
-    return kTRUE;
-}
-

trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h

@@ -7 +7 @@
 
 #include <TH2.h>
+#include <TF1.h>
 
 class MParList;
 class MMcEvt;
+class MMcEvtBasic;
 class MMcTrig;
 class MHMcCollectionArea;
+class MBinning;
 
 class MMcCollectionAreaCalc : public MTask
 {
 private:
-    const MMcEvt  *fMcEvt;
-    const MMcTrig *fMcTrig;
+    const MMcEvt       *fMcEvt;
+    const MMcEvtBasic  *fMcEvtBasic;
+    const MMcTrig      *fMcTrig;
+
+    MBinning           *fBinsTheta;
+    MBinning           *fBinsEnergy;
+    // Coarse zenith angle and energy bins used in the analysis
+
+    TF1                *fSpectrum; 
+    // Tentative energy spectrum. This modifies slightly the calculation
+    // of the effective area (see MHMcCollectionArea::Calc)
+
 
     MHMcCollectionArea *fCollArea;
@@ -23 +36 @@
     TString fObjName;
 
-    UInt_t fTotalNumSimulatedShowers;
-    UInt_t fCorsikaVersion;
-    Bool_t fAllEvtsTriggered;
-    Float_t fThetaMin;
-    Float_t fThetaMax;
-
-    Bool_t ReInit(MParList *plist);
-
-    Int_t PreProcess(MParList *pList);
-    Int_t Process();
-    Int_t PostProcess();
+    Int_t  PreProcess(MParList *pList);
+    Int_t  Process();
+    Int_t  PostProcess();
 
 public:
-    MMcCollectionAreaCalc(const char *input=NULL,
-                          const char *name=NULL, const char *title=NULL);
+    MMcCollectionAreaCalc(const char *input = NULL,
+                          const char *name = NULL, const char *title = NULL);
+
+    void SetSpectrum(TF1 *f) { fSpectrum = f; }
 
     ClassDef(MMcCollectionAreaCalc, 0) // Task to calculate the collection area histogram