Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2249) +++ trunk/MagicSoft/Mars/Changelog (revision 2250) @@ -1,3 +1,16 @@ -*-*- 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 Index: trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc (revision 2249) +++ trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc (revision 2250) @@ -237,5 +237,5 @@ // 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,32 +244,9 @@ 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,5 +259,5 @@ *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!" <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); } Index: trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h (revision 2249) +++ trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h (revision 2250) @@ -35,8 +35,10 @@ 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); Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc (revision 2249) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc (revision 2250) @@ -63,4 +63,10 @@ AddToBranchList("MMcEvt.fEnergy"); AddToBranchList("MMcEvt.fImpact"); + + fCorsikaVersion = 0; + fAllEvtsTriggered = kFALSE; + fTotalNumSimulatedShowers = 0; + fTheta = -1.; + } @@ -76,16 +82,4 @@ } - 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,6 +102,6 @@ fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers(); - *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl; + if (fTheta>=0 && fTheta!=runheader->GetTelesTheta()) @@ -116,28 +110,57 @@ *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,8 +213,8 @@ { *fLog << inf << "Total number of showers: " << fTotalNumSimulatedShowers << endl; - fCollArea->CalcEfficiency(fTotalNumSimulatedShowers, fEmin, fEmax, fSlope); - } - - return kTRUE; -} - + fCollArea->CalcEfficiency2(); + } + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h (revision 2249) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h (revision 2250) @@ -5,4 +5,6 @@ #include "MTask.h" #endif + +#include class MParList; @@ -23,10 +25,6 @@ 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);