Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc (revision 5340) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc (revision 6491) @@ -17,5 +17,6 @@ ! ! Author(s): Thomas Bretz 12/2000 -! Author(s): Harald Kornmayer 1/2001 +! Author(s): Harald Kornmayer 1/2001 +! Author(s): Abelardo Moralejo 2/2005 ! ! Copyright: MAGIC Software Development, 2000-2002 @@ -26,9 +27,5 @@ ////////////////////////////////////////////////////////////////////////////// // -// 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,5 +39,6 @@ #include "MMcEvt.hxx" -#include "MMcTrig.hxx" +#include "MMcEvtBasic.h" + #include "MMcRunHeader.hxx" #include "MMcCorsikaRunHeader.h" @@ -52,181 +50,122 @@ 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; -} - Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h (revision 5340) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h (revision 6491) @@ -7,15 +7,28 @@ #include +#include 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,19 +36,13 @@ 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