Index: trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc (revision 1715) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc (revision 1783) @@ -38,9 +38,10 @@ #include "MMcEvt.hxx" #include "MMcTrig.hxx" +#include "MMcRunHeader.hxx" +#include "MMcCorsikaRunHeader.h" ClassImp(MMcTriggerRateCalc); -void MMcTriggerRateCalc::Init(int dim, int part, float *trigbg, - float simbg, +void MMcTriggerRateCalc::Init(int dim, float *trigbg, float simbg, const char *name, const char *title) { @@ -51,8 +52,11 @@ fMcRate = NULL; + fTrigNSB = trigbg; + fSimNSB = simbg; + + fPartId = -1; + fShowers = 0; - fAnalShow = simbg; - - fPartId=part; + fAnalShow = 0; fFirst = dim>0 ? 1 : -dim; @@ -60,10 +64,10 @@ fNum = fLast-fFirst+1; - fTrigger = new float[fNum]; for (UInt_t i=0;iFindObject("MMcRunHeader"); + if (!fMcRunHeader) + { + *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl; + return kFALSE; + } + + fMcCorRunHeader = (MMcCorsikaRunHeader*) pList->FindObject("MMcCorsikaRunHeader"); + if (!fMcCorRunHeader) + { + *fLog << err << dbginf << "Error - MMcCorsikaRunHeader not found... exit." << endl; + return kFALSE; + } + + if (fMcRunHeader->GetAllEvtsTriggered()) + fShowers += fMcRunHeader->GetNumSimulatedShowers(); + + if (fMcRunHeader->GetAllEvtsTriggered()) + { + *fLog << endl << all << endl << + "WARNING: the input data file contains only the" << endl << + "events that triggered. I will assume the standard value" << endl << + "for maximum impact parameter (400 m)" < 0) + *fLog << endl << all << + "WARNING: NSB rate can be overestimated by up to 5%." << endl << + "For a precise estimate of the total rate including NSB" << endl << + "accidental triggers I need a file containing all event headers." + << endl; + else + *fLog << endl << all << + "WARNING: calculating only shower rate (no NSB accidental triggers)" << endl; + } + + *fLog << endl << all << + "WARNING: I will assume the standard maximum off axis angle" << endl << + "(5 degrees) for the original showers" << endl; + + for (UInt_t i=0; iGetAllEvtsTriggered()) + { + GetRate(i)->SetImpactMin(0.); + GetRate(i)->SetImpactMax(40000.); // in cm + } + GetRate(i)->SetSolidAngle(2.390941702e-2); // sr + + // + // Set limits of energy range, coverting from GeV to TeV: + // + GetRate(i)->SetEnergyMin(fMcCorRunHeader->GetELowLim()/1000.); + GetRate(i)->SetEnergyMax(fMcCorRunHeader->GetEUppLim()/1000.); + + } + + return kTRUE; } @@ -79,16 +148,15 @@ // // dim: fDimension -// part: fPartId -// *trigbg: number of shower from bacground that triggers +// *trigbg: number of NSB triggers for // a given trigger condition. -// simbg: Number of simulated showers for the background +// simbg: Number of simulated events for the NSB // rate: rate of incident showers // -MMcTriggerRateCalc::MMcTriggerRateCalc(float rate, int dim, int part, +MMcTriggerRateCalc::MMcTriggerRateCalc(float rate, int dim, float *trigbg, float simbg, const char *name, const char *title) { - Init(dim, part, trigbg, simbg, name, title); + Init(dim, trigbg, simbg, name, title); } @@ -99,14 +167,12 @@ // // dim: fDimension -// part: fPartId -// *trigbg: number of shower from bacground that triggers +// *trigbg: number of NSB triggers for // a given trigger condition. -// simbg: Number of simulated showers for the background -// -MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, int part, float *trigbg, - float simbg, +// simbg: Number of simulated events for the NSB +// +MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, float *trigbg,float simbg, const char *name, const char *title) { - Init(dim, part, trigbg, simbg, name, title); + Init(dim, trigbg, simbg, name, title); } @@ -160,4 +226,78 @@ } + for (UInt_t i=0;iGetAllEvtsTriggered() && + fMcEvt->GetPhotElfromShower() ) + fAnalShow++; + + // + // Set PartId and check it is the same for all events + // + if (fPartId < 0) + fPartId = fMcEvt->GetPartId(); + else if (fPartId != fMcEvt->GetPartId()) + { + *fLog << err << dbginf << "Incident particle type seems to change"; + *fLog << "from " << fPartId << " to " << fMcEvt->GetPartId() << endl; + *fLog << "in input data files... aborting." << endl; + return kFALSE; + } + + // + // Getting angles, energy and impact parameter to set boundaries + // + const Float_t theta = fMcEvt->GetTheta(); + const Float_t phi = fMcEvt->GetPhi(); + const Float_t param = fMcEvt->GetImpact(); + const Float_t ener = fMcEvt->GetEnergy()/1000.0; + + // + // Counting number of triggers + // + for (UInt_t i=0; iGetFirstLevel()) + fTrigger[i] ++; + + if ( ! fMcRunHeader->GetAllEvtsTriggered()) + GetRate(i)->UpdateBoundaries(ener, theta, phi, param); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// The PostProcess-function calculates and shows the trigger rate +// +Bool_t MMcTriggerRateCalc::PostProcess() +{ for (UInt_t i=0; iGetSlopeSpec() != -2.75) + { + *fLog << err << dbginf << + "Spectrum slope as read from input file ("<< + fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) for protons" << endl << "... aborting." << endl; + return kFALSE; + } + rate.SetFlux(0.1091, 2.75); + break; case kHELIUM: - rate.SetFlux(0.0660, 2.62); - break; + if (fMcCorRunHeader->GetSlopeSpec() != -2.62) + { + *fLog << err << dbginf << + "Spectrum slope as read from input file ("<< + fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) for Helium" << endl << "... aborting." << endl; + return kFALSE; + } + rate.SetFlux(0.0660, 2.62); + break; default: - *fLog << err << dbginf << "Unknown incident flux parameters for "; - *fLog << fPartId<< " particle Id ... aborting." << endl; - return kFALSE; + *fLog << err << dbginf << "Unknown incident flux parameters for "; + *fLog << fPartId<< " particle Id ... aborting." << endl; + return kFALSE; } - rate.SetBackground(fTrigger[i], fAnalShow); - - fTrigger[i]=0; - } - - fAnalShow=0.0; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// The Process-function counts the number of simulated showers, the -// number of analised showers and the number of triggers. It also updates -// the limits for theta, phi, energy and impact parameter in the -// MHMcRate container. -// -Bool_t MMcTriggerRateCalc::Process() -{ - // - // Counting analysed and simulated showers - // - fShowers++; - if (fMcEvt->GetPhotElfromShower()) - fAnalShow++; - - // - // Getting angles, energy and impact parameter to set boundaries - // - const Float_t theta = fMcEvt->GetTheta(); - const Float_t phi = fMcEvt->GetPhi(); - const Float_t param = fMcEvt->GetImpact(); - const Float_t ener = fMcEvt->GetEnergy()/1000.0; - - // - // Counting number of triggers - // - for (UInt_t i=0; iGetFirstLevel()) - fTrigger[i] ++; - - GetRate(i)->UpdateBoundaries(ener, theta, phi, param); - } - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// The PostProcess-function calculates and shows the trigger rate -// -Bool_t MMcTriggerRateCalc::PostProcess() -{ + } + // // Computing trigger rate Index: trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.h (revision 1715) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.h (revision 1783) @@ -11,4 +11,6 @@ class MParList; class MMcEvt; +class MMcRunHeader; +class MMcCorsikaRunHeader; class MMcTrig; class MHMcRate; @@ -22,10 +24,15 @@ TObjArray *fMcTrig; + MMcRunHeader *fMcRunHeader; + MMcCorsikaRunHeader *fMcCorRunHeader; + UInt_t fNum; // decoded dimension UInt_t fFirst; UInt_t fLast; - Float_t* fTrigger; // Number of triggered showers + Float_t* fTrigNSB; // Number of triggers due to NSB alone + Float_t fSimNSB; // Number of simulated NSB-only events + Float_t* fTrigger; // Number of triggered showers Float_t fShowers; // Number of simulated showers Float_t fAnalShow; // Number of analysed showers @@ -33,5 +40,5 @@ Int_t fPartId; // Incident particle that generates showers - void Init(int dim, int part, float *trigbg, + void Init(int dim, float *trigbg, float simbg, const char *name, const char *title); @@ -40,13 +47,13 @@ public: - MMcTriggerRateCalc(int dim=0, int part=14, float *trigbg=NULL, - float simbg=100000, + MMcTriggerRateCalc(int dim=0, float *trigbg=NULL, float simbg=100000, const char *name=NULL, const char *title=NULL); - MMcTriggerRateCalc(float rate, int dim, int part, float *trigbg, - float simbg, + MMcTriggerRateCalc(float rate, int dim, float *trigbg, float simbg, const char *name=NULL, const char *title=NULL); ~MMcTriggerRateCalc(); + + Bool_t ReInit(MParList *plist); Bool_t PreProcess(MParList *pList);