source: trunk/Mars/mhistmc/MHMcRate.cc@ 18282

Last change on this file since 18282 was 2173, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 8.2 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 12/2000 <mailto:tbretz@uni-sw.gwdg.de>
19! Author(s): Harald Kornmayer 1/2001
20! Author(s): Abelardo Moralejo 2/2003
21!
22! Copyright: MAGIC Software Development, 2000-2001
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// Explanations on the rate calculation can be found in
30// chapter 7 of the following diploma thesis:
31// http://www.pd.infn.it/magic/tesi2.ps.gz (in Italian)
32//
33////////////////////////////////////////////////////////////////////////////
34
35#include "MHMcRate.h"
36
37#include "MLog.h"
38#include "MLogManip.h"
39
40ClassImp(MHMcRate);
41
42using namespace std;
43
44void MHMcRate::Init(const char *name, const char *title)
45{
46 fName = name ? name : "MMcTriggerRate";
47 fTitle = title ? title : "Task to calc the collection area ";
48
49 fPartId=0; // Type of particle
50
51 fEnergyMax=0.0; // Maximum Energy (TeV)
52 fEnergyMin=1000000.0; // Minimum Energy (TeV)
53
54 fSolidAngle = -1.; // Solid angle within which incident directions
55 // are distributed
56 fThetaMax=0.0; // Maximum theta angle of run
57 fThetaMin=370.0; // Minimum theta angle of run
58 fPhiMax=0.0; // Maximum phi angle of run
59 fPhiMin=370.0; // Minimum phi angle of run
60
61 fImpactMax=0.0; // Maximum impact parameter
62 fImpactMin=100000.0; // Minimum impact parameter
63
64 fBackTrig=-1.0; // Number of triggers from background
65 fBackSim=-1.0; // Number of simulated showers for the background
66
67 fTriggerRate= -1.0; // Trigger rate in Hz
68 fTriggerRateError= -1.0; // Estimated error for the trigger rate in Hz
69
70 fMultiplicity = -1; // Multiplicity of the trigger condition
71 fMeanThreshold = -1.; // Discr. threshold of the trigger condition
72 fTriggerCondNum = 0; // Trigger condition number within camera file
73
74}
75
76// --------------------------------------------------------------------------
77//
78// default constructor
79// fills all member data with initial values
80//
81MHMcRate::MHMcRate(const char *name, const char *title)
82{
83 Init(name, title);
84
85 fSpecIndex=0.0; // dn/dE = k * e^{- fSpecIndex}
86 fFlux0=-1.0; // dn/dE = fFlux0 * E^{-a}
87
88 fShowerRate= -1.0; // Showers rate in Hz
89 fShowerRateError=0.0; // Estimated error of shower rate in Hz
90}
91
92// --------------------------------------------------------------------------
93//
94// overloaded constructor I
95// fills all member data with initial values and sets the rate of
96// incident showers to ShowRate
97//
98MHMcRate::MHMcRate(Float_t showrate,
99 const char *name, const char *title)
100{
101 Init(name, title);
102
103 fSpecIndex=0.0; // dn/dE = k * e^{- fSpecIndex}
104 fFlux0=-1.0; // dn/dE = fFlux0 * E^{-a}
105
106 fShowerRate= showrate; // Showers rate in Hz
107 fShowerRateError=sqrt(showrate); // Estimated error of shower rate in Hz
108}
109
110// --------------------------------------------------------------------------
111//
112// overloaded constructor I
113// fills all member data with initial values and sets the
114// spectral index and the initial flux to SpecIndex and Flux0
115//
116MHMcRate::MHMcRate(Float_t specindex, Float_t flux0,
117 const char *name, const char *title)
118{
119 Init(name, title);
120
121 fSpecIndex=specindex; // dn/dE = k * e^{- fSpecIndex}
122 fFlux0=flux0; // dn/dE = fFlux0 * E^{-a}
123
124 fShowerRate= -1.0;
125 fShowerRateError=0.0;
126}
127
128// --------------------------------------------------------------------------
129//
130// set the particle that produces the showers in the athmosphere
131//
132void MHMcRate:: SetParticle(UShort_t part)
133{
134 fPartId=part;
135}
136
137// --------------------------------------------------------------------------
138//
139// Set the information about trigger due only to the Night Sky Background:
140//
141void MHMcRate::SetBackground (Float_t showers, Float_t triggers)
142{
143 fBackTrig=showers; // Number of triggers from background
144 fBackSim=triggers; // Number of simulated showers for the background
145}
146
147// --------------------------------------------------------------------------
148//
149// set the parameters to compute the incident rate
150//
151void MHMcRate:: SetFlux(Float_t flux0, Float_t specindx)
152{
153 fFlux0=flux0;
154 fSpecIndex=specindx;
155
156}
157
158// --------------------------------------------------------------------------
159//
160// set the incident rate
161//
162void MHMcRate:: SetIncidentRate(Float_t showerrate)
163{
164 fShowerRate=showerrate;
165}
166
167// --------------------------------------------------------------------------
168//
169// update the limits for energy, theta, phi and impact parameter
170//
171void MHMcRate::UpdateBoundaries(Float_t energy, Float_t theta,
172 Float_t phi, Float_t impact)
173{
174 // It updates the limit values
175
176 if (fThetaMax<theta) fThetaMax=theta;
177 if (fThetaMin>theta) fThetaMin=theta;
178
179 if (fPhiMax<phi) fPhiMax=phi;
180 if (fPhiMin>phi) fPhiMin=phi;
181
182 if (fImpactMax<impact) fImpactMax=impact;
183 if (fImpactMin>impact) fImpactMin=impact;
184
185 if (fEnergyMax<energy) fEnergyMax=energy;
186 if (fEnergyMin>energy) fEnergyMin=energy;
187}
188
189// --------------------------------------------------------------------------
190//
191// compute the trigger rate and set the ReadyToSave bit
192//
193void MHMcRate::CalcRate(Float_t trig, Float_t anal, Float_t simu)
194{
195 // It computes the trigger rate
196
197 // First one computes the rate of incident showers.
198 const Double_t specidx = 1.0-fSpecIndex;
199
200 const Double_t epowmax = pow((double)fEnergyMax, specidx);
201 const Double_t epowmin = pow((double)fEnergyMin, specidx);
202
203 if (fShowerRate <= 0)
204 fShowerRate = fFlux0/specidx*(epowmax-epowmin);
205
206 if (fSolidAngle < 0.)
207 fSolidAngle = (fPhiMax-fPhiMin)*(cos(fThetaMin)-cos(fThetaMax));
208
209 if (fPartId!=1)
210 fShowerRate *= fSolidAngle;
211
212 fShowerRate *= TMath::Pi()*(fImpactMax/100.0*fImpactMax/100.0 -
213 fImpactMin/100.0*fImpactMin/100.0);
214
215 fShowerRateError = sqrt(fShowerRate);
216
217 // The simulated trigger time in the camera program is 160 ns:
218 // 9/10/2002, AM: Fixed error below in calculation of "anal2"
219 // ( added factor fShowerRate/simu )
220
221 const Double_t anal2 = 1.0-fShowerRate*(anal/simu)*160.0e-9;
222 const Double_t back2 = fBackSim*160.0e-9;
223
224 // Then the trigger rate and its error is evaluated
225 if(fBackTrig<0){
226 fTriggerRateError = sqrt((trig*fShowerRate*fShowerRate/(simu*simu)) +
227 (anal2*anal2*1/(fBackSim*back2*back2)));
228 fBackTrig=0;
229 }
230 else
231 fTriggerRateError = sqrt((trig*fShowerRate*fShowerRate/(simu*simu)) +
232 (anal2*anal2*fBackTrig/(back2*back2)));
233
234 fTriggerRate = trig*fShowerRate/simu + anal2*fBackTrig/back2;
235
236 SetReadyToSave();
237}
238
239// --------------------------------------------------------------------------
240//
241// print the trigger rate
242//
243void MHMcRate::Print(Option_t *) const
244{
245 *fLog << all << "Incident rate " << fShowerRate << " Hz " << endl;
246 *fLog << "Multiplicity: " << fMultiplicity << ", Discr. threshold: " << fMeanThreshold << endl;
247 *fLog << "Trigger Rate " << fTriggerRate << " +- " << fTriggerRateError << " Hz" << endl;
248}
249
250// --------------------------------------------------------------------------
251//
252// draw the trigger rate
253//
254void MHMcRate::Draw(Option_t *)
255{
256 *fLog << all << dbginf << " - MHMcRate::Draw: To be iplemented" << endl;
257}
258
Note: See TracBrowser for help on using the repository browser.