Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 1782) +++ trunk/MagicSoft/Mars/Changelog (revision 1783) @@ -1,4 +1,20 @@ -*-*- END -*-*- + 2003/02/21: Abelardo Moralejo + + * mmontecarlo/MMcTriggerRateCalc.[h,cc] + - adapted to new camera files, added warnings. + - added ReInit procedure to read relevant info from from the run headers + + * mhist/MHMcRate.[h,cc] + - adapted accordingly. Added Set functions for several members. + + * mmc/MMcCorsikaRunHeader.h + - added Get functions for fELowLim, fEUppLim and fSlopeSpec. + + * mmain/MMontecarlo.cc, macros/trigrate.C + - adapted to changes above, changed MReadTree to MReadMarsFile to + be able to read the run headers. + 2003/02/21: Antonio Stamerra Index: trunk/MagicSoft/Mars/macros/trigrate.C =================================================================== --- trunk/MagicSoft/Mars/macros/trigrate.C (revision 1782) +++ trunk/MagicSoft/Mars/macros/trigrate.C (revision 1783) @@ -124,5 +124,6 @@ // analised trigger conditions should be set (BgR[]) // - MReadTree reader("Events", filename); + + MReadMarsFile reader("Events", filename); tasklist.AddToList(&reader); @@ -141,5 +142,6 @@ cout << "Number of Trigger conditions: " << num << endl; - MMcTriggerRateCalc rate(dim, kPROTON, BgR, numnsbevents); + MMcTriggerRateCalc rate(dim, BgR, numnsbevents); + tasklist.AddToList(&rate); Index: trunk/MagicSoft/Mars/mhist/MHMcRate.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMcRate.cc (revision 1782) +++ trunk/MagicSoft/Mars/mhist/MHMcRate.cc (revision 1783) @@ -18,4 +18,9 @@ ! Author(s): Thomas Bretz 12/2000 ! Author(s): Harald Kornmayer 1/2001 +! Author(s): Abelardo Moralejo 2/2003 +! +! Explanations on the rate calculation can be found in +! chapter 7 of the following diploma thesis: +! http://www.pd.infn.it/magic/tesi2.ps.gz (in Italian) ! ! Copyright: MAGIC Software Development, 2000-2001 @@ -38,7 +43,9 @@ fPartId=0; // Type of particle - fEnergyMax=0.0; // Maximum Energy in GeV - fEnergyMin=1000000.0; // Minimum Energy in GeV - + fEnergyMax=0.0; // Maximum Energy (TeV) + fEnergyMin=1000000.0; // Minimum Energy (TeV) + + fSolidAngle = -1.; // Solid angle within which incident directions + // are distributed fThetaMax=0.0; // Maximum theta angle of run fThetaMin=370.0; // Minimum theta angle of run @@ -119,5 +126,5 @@ // -------------------------------------------------------------------------- // -// Set the information about trigger due only to the background conditions +// Set the information about trigger due only to the Night Sky Background: // void MHMcRate::SetBackground (Float_t showers, Float_t triggers) @@ -184,12 +191,14 @@ if (fShowerRate <= 0) - fShowerRate = fFlux0/specidx*(epowmin-epowmax); + fShowerRate = fFlux0/specidx*(epowmax-epowmin); + + if (fSolidAngle < 0.) + fSolidAngle = (fPhiMax-fPhiMin)*(cos(fThetaMin)-cos(fThetaMax)); if (fPartId!=1) - fShowerRate *= (fPhiMax-fPhiMin)*(cos(fThetaMax)-cos(fThetaMin)); - - const Double_t impactdiff = fImpactMax-fImpactMin; - - fShowerRate *= TMath::Pi()*(impactdiff/100.0*impactdiff/100.0); + fShowerRate *= fSolidAngle; + + fShowerRate *= TMath::Pi()*(fImpactMax/100.0*fImpactMax/100.0 - + fImpactMin/100.0*fImpactMin/100.0); fShowerRateError = sqrt(fShowerRate); Index: trunk/MagicSoft/Mars/mhist/MHMcRate.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMcRate.h (revision 1782) +++ trunk/MagicSoft/Mars/mhist/MHMcRate.h (revision 1783) @@ -12,6 +12,6 @@ UShort_t fPartId; // Type of particle - Float_t fEnergyMax; // Maximum Energy in GeV - Float_t fEnergyMin; // Minimum Energy in GeV + Float_t fEnergyMax; // Maximum Energy in TeV + Float_t fEnergyMin; // Minimum Energy in TeV Float_t fThetaMax; // Maximum theta angle of run @@ -20,6 +20,9 @@ Float_t fPhiMin; // Minimum phi angle of run - Float_t fImpactMax; // Maximum impact parameter - Float_t fImpactMin; // Minimum impact parameter + Float_t fSolidAngle; // Solid angle within which incident directions + // are distributed (sr) + + Float_t fImpactMax; // Maximum impact parameter (cm) + Float_t fImpactMin; // Minimum impact parameter (cm) Float_t fBackTrig; // Number of triggers from background @@ -50,4 +53,16 @@ void SetIncidentRate(Float_t showerrate); + void SetImpactMax(Float_t Impact) {fImpactMax=Impact;} + void SetImpactMin(Float_t Impact) {fImpactMin=Impact;} + + void SetThetaMax(Float_t Theta) {fThetaMax=Theta;} + void SetThetaMin(Float_t Theta) {fThetaMin=Theta;} + void SetPhiMax(Float_t Phi) {fPhiMax=Phi;} + void SetPhiMin(Float_t Phi) {fPhiMin=Phi;} + + void SetSolidAngle(Float_t Solid) {fSolidAngle=Solid;} + void SetEnergyMax(Float_t Energy) {fEnergyMax=Energy;} + void SetEnergyMin(Float_t Energy) {fEnergyMin=Energy;} + void UpdateBoundaries(Float_t energy, Float_t theta, Float_t phi, Float_t impact); Index: trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc =================================================================== --- trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc (revision 1782) +++ trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc (revision 1783) @@ -349,10 +349,12 @@ // analised trigger conditions should be set (BgR[]) // - MReadTree read("Events", fInputFile); + MReadMarsFile read("Events", fInputFile); tlist.AddToList(&read); - Float_t BgR[10]={660, 4, 0, 0, 0, 0, 0, 0, 0, 0}; - - MMcTriggerRateCalc crate(dim, 14, BgR, 100000); + // We calculate only shower rate (not including NSB-only triggers) + Float_t* BgR = new Float_t[dim]; + memset(BgR, 0, num*sizeof(Float_t)); + + MMcTriggerRateCalc crate(dim, BgR, 100000); tlist.AddToList(&crate); Index: trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc (revision 1782) +++ 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 1782) +++ 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); Index: trunk/MagicSoft/include-Classes/MMcFormat/MMcCorsikaRunHeader.h =================================================================== --- trunk/MagicSoft/include-Classes/MMcFormat/MMcCorsikaRunHeader.h (revision 1782) +++ trunk/MagicSoft/include-Classes/MMcFormat/MMcCorsikaRunHeader.h (revision 1783) @@ -79,4 +79,8 @@ ); + Float_t GetELowLim() const { return fELowLim; } + Float_t GetEUppLim() const { return fEUppLim; } + Float_t GetSlopeSpec() const { return fSlopeSpec; } + ClassDef(MMcCorsikaRunHeader, 2) // storage container for corsika setup information };