/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Thomas Bretz 12/2000 ! Author(s): Harald Kornmayer 1/2001 ! Author(s): Abelardo Moralejo 2/2003 ! ! Copyright: MAGIC Software Development, 2000-2001 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // 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) // //////////////////////////////////////////////////////////////////////////// #include "MHMcRate.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHMcRate); void MHMcRate::Init(const char *name, const char *title) { fName = name ? name : "MMcTriggerRate"; fTitle = title ? title : "Task to calc the collection area "; fPartId=0; // Type of particle 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 fPhiMax=0.0; // Maximum phi angle of run fPhiMin=370.0; // Minimum phi angle of run fImpactMax=0.0; // Maximum impact parameter fImpactMin=100000.0; // Minimum impact parameter fBackTrig=-1.0; // Number of triggers from background fBackSim=-1.0; // Number of simulated showers for the background fTriggerRate= -1.0; // Trigger rate in Hz fTriggerRateError= -1.0; // Estimated error for the trigger rate in Hz fMultiplicity = -1; // Multiplicity of the trigger condition fMeanThreshold = -1.; // Discr. threshold of the trigger condition fTriggerCondNum = 0; // Trigger condition number within camera file } // -------------------------------------------------------------------------- // // default constructor // fills all member data with initial values // MHMcRate::MHMcRate(const char *name, const char *title) { Init(name, title); fSpecIndex=0.0; // dn/dE = k * e^{- fSpecIndex} fFlux0=-1.0; // dn/dE = fFlux0 * E^{-a} fShowerRate= -1.0; // Showers rate in Hz fShowerRateError=0.0; // Estimated error of shower rate in Hz } // -------------------------------------------------------------------------- // // overloaded constructor I // fills all member data with initial values and sets the rate of // incident showers to ShowRate // MHMcRate::MHMcRate(Float_t showrate, const char *name, const char *title) { Init(name, title); fSpecIndex=0.0; // dn/dE = k * e^{- fSpecIndex} fFlux0=-1.0; // dn/dE = fFlux0 * E^{-a} fShowerRate= showrate; // Showers rate in Hz fShowerRateError=sqrt(showrate); // Estimated error of shower rate in Hz } // -------------------------------------------------------------------------- // // overloaded constructor I // fills all member data with initial values and sets the // spectral index and the initial flux to SpecIndex and Flux0 // MHMcRate::MHMcRate(Float_t specindex, Float_t flux0, const char *name, const char *title) { Init(name, title); fSpecIndex=specindex; // dn/dE = k * e^{- fSpecIndex} fFlux0=flux0; // dn/dE = fFlux0 * E^{-a} fShowerRate= -1.0; fShowerRateError=0.0; } // -------------------------------------------------------------------------- // // set the particle that produces the showers in the athmosphere // void MHMcRate:: SetParticle(UShort_t part) { fPartId=part; } // -------------------------------------------------------------------------- // // Set the information about trigger due only to the Night Sky Background: // void MHMcRate::SetBackground (Float_t showers, Float_t triggers) { fBackTrig=showers; // Number of triggers from background fBackSim=triggers; // Number of simulated showers for the background } // -------------------------------------------------------------------------- // // set the parameters to compute the incident rate // void MHMcRate:: SetFlux(Float_t flux0, Float_t specindx) { fFlux0=flux0; fSpecIndex=specindx; } // -------------------------------------------------------------------------- // // set the incident rate // void MHMcRate:: SetIncidentRate(Float_t showerrate) { fShowerRate=showerrate; } // -------------------------------------------------------------------------- // // update the limits for energy, theta, phi and impact parameter // void MHMcRate::UpdateBoundaries(Float_t energy, Float_t theta, Float_t phi, Float_t impact) { // It updates the limit values if (fThetaMaxtheta) fThetaMin=theta; if (fPhiMaxphi) fPhiMin=phi; if (fImpactMaximpact) fImpactMin=impact; if (fEnergyMaxenergy) fEnergyMin=energy; } // -------------------------------------------------------------------------- // // compute the trigger rate and set the ReadyToSave bit // void MHMcRate::CalcRate(Float_t trig, Float_t anal, Float_t simu) { // It computes the trigger rate // First one computes the rate of incident showers. const Double_t specidx = 1.0-fSpecIndex; const Double_t epowmax = pow(fEnergyMax, specidx); const Double_t epowmin = pow(fEnergyMin, specidx); if (fShowerRate <= 0) fShowerRate = fFlux0/specidx*(epowmax-epowmin); if (fSolidAngle < 0.) fSolidAngle = (fPhiMax-fPhiMin)*(cos(fThetaMin)-cos(fThetaMax)); if (fPartId!=1) fShowerRate *= fSolidAngle; fShowerRate *= TMath::Pi()*(fImpactMax/100.0*fImpactMax/100.0 - fImpactMin/100.0*fImpactMin/100.0); fShowerRateError = sqrt(fShowerRate); // The simulated trigger time in the camera program is 160 ns: // 9/10/2002, AM: Fixed error below in calculation of "anal2" // ( added factor fShowerRate/simu ) const Double_t anal2 = 1.0-fShowerRate*(anal/simu)*160.0e-9; const Double_t back2 = fBackSim*160.0e-9; // Then the trigger rate and its error is evaluated if(fBackTrig<0){ fTriggerRateError = sqrt((trig*fShowerRate*fShowerRate/(simu*simu)) + (anal2*anal2*1/(fBackSim*back2*back2))); fBackTrig=0; } else fTriggerRateError = sqrt((trig*fShowerRate*fShowerRate/(simu*simu)) + (anal2*anal2*fBackTrig/(back2*back2))); fTriggerRate = trig*fShowerRate/simu + anal2*fBackTrig/back2; SetReadyToSave(); } // -------------------------------------------------------------------------- // // print the trigger rate // void MHMcRate::Print(Option_t *) const { *fLog << all << "Incident rate " << fShowerRate << " Hz " << endl; *fLog << "Multiplicity: " << fMultiplicity << ", Discr. threshold: " << fMeanThreshold << endl; *fLog << "Trigger Rate " << fTriggerRate << " +- " << fTriggerRateError << " Hz" << endl; } // -------------------------------------------------------------------------- // // draw the trigger rate // void MHMcRate::Draw(Option_t *) { *fLog << all << dbginf << " - MHMcRate::Draw: To be iplemented" << endl; } q