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 | !
|
---|
21 | ! Copyright: MAGIC Software Development, 2000-2001
|
---|
22 | !
|
---|
23 | ! Modified 4/7/2002 Abelardo Moralejo: now the dimension of fTrigger is
|
---|
24 | ! set dinamically, to allow an arbitrary large number of trigger
|
---|
25 | ! conditions to be processed.
|
---|
26 | !
|
---|
27 | !
|
---|
28 | \* ======================================================================== */
|
---|
29 |
|
---|
30 | #include "MMcTriggerRateCalc.h"
|
---|
31 |
|
---|
32 | #include "MLog.h"
|
---|
33 | #include "MLogManip.h"
|
---|
34 |
|
---|
35 | #include "MParList.h"
|
---|
36 | #include "MHMcRate.h"
|
---|
37 |
|
---|
38 | #include "MMcEvt.hxx"
|
---|
39 | #include "MMcTrig.hxx"
|
---|
40 | #include "MMcRunHeader.hxx"
|
---|
41 | #include "MMcTrigHeader.hxx"
|
---|
42 | #include "MMcCorsikaRunHeader.h"
|
---|
43 |
|
---|
44 | #include "MH.h"
|
---|
45 | #include <TCanvas.h>
|
---|
46 |
|
---|
47 | ClassImp(MMcTriggerRateCalc);
|
---|
48 |
|
---|
49 | void MMcTriggerRateCalc::Init(int dim, float *trigbg, float simbg,
|
---|
50 | const char *name, const char *title)
|
---|
51 | {
|
---|
52 | fName = name ? name : "MMcTriggerRateCalc";
|
---|
53 | fTitle = title ? title : "Task to calc the trigger rate ";
|
---|
54 |
|
---|
55 | fMcTrig = NULL;
|
---|
56 | fMcRate = NULL;
|
---|
57 |
|
---|
58 | fTrigNSB = trigbg;
|
---|
59 | fSimNSB = simbg;
|
---|
60 |
|
---|
61 | fPartId = -1;
|
---|
62 |
|
---|
63 | fShowers = 0;
|
---|
64 | fAnalShow = 0;
|
---|
65 |
|
---|
66 | fFirst = dim>0 ? 1 : -dim;
|
---|
67 | fLast = dim>0 ? dim : -dim;
|
---|
68 |
|
---|
69 | fNum = fLast-fFirst+1;
|
---|
70 | fTrigger = new float[fNum];
|
---|
71 |
|
---|
72 | for (UInt_t i=0;i<fNum;i++)
|
---|
73 | fTrigger[i]=0;
|
---|
74 |
|
---|
75 | AddToBranchList("MMcEvt.fPartId");
|
---|
76 | AddToBranchList("MMcEvt.fImpact");
|
---|
77 | AddToBranchList("MMcEvt.fEnergy");
|
---|
78 | AddToBranchList("MMcEvt.fPhi");
|
---|
79 | AddToBranchList("MMcEvt.fTheta");
|
---|
80 | AddToBranchList("MMcEvt.fPhotElfromShower");
|
---|
81 | AddToBranchList("MMcTrig", "fNumFirstLevel", fFirst, fLast);
|
---|
82 |
|
---|
83 | }
|
---|
84 |
|
---|
85 | // ReInit: read run header to get some info later:
|
---|
86 |
|
---|
87 | Bool_t MMcTriggerRateCalc::ReInit(MParList *pList)
|
---|
88 | {
|
---|
89 | fMcRunHeader = (MMcRunHeader*) pList->FindObject("MMcRunHeader");
|
---|
90 | if (!fMcRunHeader)
|
---|
91 | {
|
---|
92 | *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl;
|
---|
93 | return kFALSE;
|
---|
94 | }
|
---|
95 |
|
---|
96 | fMcCorRunHeader = (MMcCorsikaRunHeader*) pList->FindObject("MMcCorsikaRunHeader");
|
---|
97 | if (!fMcCorRunHeader)
|
---|
98 | {
|
---|
99 | *fLog << err << dbginf << "Error - MMcCorsikaRunHeader not found... exit." << endl;
|
---|
100 | return kFALSE;
|
---|
101 | }
|
---|
102 |
|
---|
103 | fShowers += fMcRunHeader->GetNumSimulatedShowers();
|
---|
104 |
|
---|
105 | if (fMcRunHeader->GetAllEvtsTriggered())
|
---|
106 | {
|
---|
107 | *fLog << endl << all << endl <<
|
---|
108 | "WARNING: the input data file contains only the" << endl <<
|
---|
109 | "events that triggered. I will assume the standard value" << endl <<
|
---|
110 | "for maximum impact parameter (400 m)" <<endl;
|
---|
111 |
|
---|
112 |
|
---|
113 | if (fTrigNSB[0] > 0)
|
---|
114 | *fLog << endl << all <<
|
---|
115 | "WARNING: NSB rate can be overestimated by up to 5%." << endl <<
|
---|
116 | "For a precise estimate of the total rate including NSB" << endl <<
|
---|
117 | "accidental triggers I need a file containing all event headers."
|
---|
118 | << endl;
|
---|
119 | else
|
---|
120 | *fLog << endl << all <<
|
---|
121 | "WARNING: calculating only shower rate (no NSB accidental triggers)" << endl;
|
---|
122 | }
|
---|
123 |
|
---|
124 | *fLog << endl << all <<
|
---|
125 | "WARNING: I will assume the standard maximum off axis angle" << endl <<
|
---|
126 | "(5 degrees) for the original showers" << endl;
|
---|
127 |
|
---|
128 |
|
---|
129 | for (UInt_t i=0; i<fNum; i++)
|
---|
130 | {
|
---|
131 | if (fMcRunHeader->GetAllEvtsTriggered())
|
---|
132 | {
|
---|
133 | GetRate(i)->SetImpactMin(0.);
|
---|
134 | GetRate(i)->SetImpactMax(40000.); // in cm
|
---|
135 | }
|
---|
136 | GetRate(i)->SetSolidAngle(2.390941702e-2); // sr
|
---|
137 |
|
---|
138 | //
|
---|
139 | // Set limits of energy range, coverting from GeV to TeV:
|
---|
140 | //
|
---|
141 | GetRate(i)->SetEnergyMin(fMcCorRunHeader->GetELowLim()/1000.);
|
---|
142 | GetRate(i)->SetEnergyMax(fMcCorRunHeader->GetEUppLim()/1000.);
|
---|
143 |
|
---|
144 | TString th("MMcTrigHeader;");
|
---|
145 | th += i+1;
|
---|
146 |
|
---|
147 | MMcTrigHeader* trighead = (MMcTrigHeader*) pList->FindObject(th);
|
---|
148 | GetRate(i)->SetMeanThreshold(trighead->GetMeanThreshold());
|
---|
149 | GetRate(i)->SetMultiplicity(trighead->GetMultiplicity());
|
---|
150 |
|
---|
151 | }
|
---|
152 |
|
---|
153 | return kTRUE;
|
---|
154 | }
|
---|
155 |
|
---|
156 | // --------------------------------------------------------------------------
|
---|
157 | //
|
---|
158 | // overloaded constructor I
|
---|
159 | //
|
---|
160 | // dim: fDimension
|
---|
161 | // *trigbg: number of NSB triggers for
|
---|
162 | // a given trigger condition.
|
---|
163 | // simbg: Number of simulated events for the NSB
|
---|
164 | // rate: rate of incident showers
|
---|
165 | //
|
---|
166 |
|
---|
167 | MMcTriggerRateCalc::MMcTriggerRateCalc(float rate, int dim,
|
---|
168 | float *trigbg, float simbg,
|
---|
169 | const char *name, const char *title)
|
---|
170 | {
|
---|
171 | Init(dim, trigbg, simbg, name, title);
|
---|
172 | }
|
---|
173 |
|
---|
174 |
|
---|
175 | // --------------------------------------------------------------------------
|
---|
176 | //
|
---|
177 | // overloaded constructor II
|
---|
178 | //
|
---|
179 | // dim: fDimension
|
---|
180 | // *trigbg: number of NSB triggers for
|
---|
181 | // a given trigger condition.
|
---|
182 | // simbg: Number of simulated events for the NSB
|
---|
183 | //
|
---|
184 | MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, float *trigbg,float simbg,
|
---|
185 | const char *name, const char *title)
|
---|
186 | {
|
---|
187 | Init(dim, trigbg, simbg, name, title);
|
---|
188 | }
|
---|
189 |
|
---|
190 | MMcTriggerRateCalc::~MMcTriggerRateCalc()
|
---|
191 | {
|
---|
192 | if (fMcTrig)
|
---|
193 | delete fMcTrig;
|
---|
194 |
|
---|
195 | if (fMcRate)
|
---|
196 | delete fMcRate;
|
---|
197 | }
|
---|
198 |
|
---|
199 |
|
---|
200 | // --------------------------------------------------------------------------
|
---|
201 | //
|
---|
202 | // The PreProcess connects the raw data with this task. It checks if the
|
---|
203 | // input containers exist, if not a kFalse flag is returned. It also checks
|
---|
204 | // if the output contaniers exist, if not they are created.
|
---|
205 | // This task can read either Montecarlo files with multiple trigger
|
---|
206 | // options, either Montecarlo files with a single trigger option.
|
---|
207 | //
|
---|
208 | Bool_t MMcTriggerRateCalc::PreProcess (MParList *pList)
|
---|
209 | {
|
---|
210 | // connect the raw data with this task
|
---|
211 |
|
---|
212 | fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
|
---|
213 | if (!fMcEvt)
|
---|
214 | {
|
---|
215 | *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
|
---|
216 | return kFALSE;
|
---|
217 | }
|
---|
218 |
|
---|
219 | UInt_t num;
|
---|
220 |
|
---|
221 | fMcTrig = new TObjArray(pList->FindObjectList("MMcTrig", fFirst, fLast));
|
---|
222 | num = fMcTrig->GetEntriesFast();
|
---|
223 | if (num != fNum)
|
---|
224 | {
|
---|
225 | *fLog << err << dbginf << fNum << " MMcTrig objects requested, ";
|
---|
226 | *fLog << num << " are available... aborting." << endl;
|
---|
227 | return kFALSE;
|
---|
228 | }
|
---|
229 |
|
---|
230 | fMcRate = new TObjArray(pList->FindObjectList("MHMcRate", fFirst, fLast));
|
---|
231 | num = fMcRate->GetEntriesFast();
|
---|
232 | if (num != fNum)
|
---|
233 | {
|
---|
234 | *fLog << err << dbginf << fNum << " MHMcRate objects requested, ";
|
---|
235 | *fLog << num << " are available... aborting." << endl;
|
---|
236 | return kFALSE;
|
---|
237 | }
|
---|
238 |
|
---|
239 | for (UInt_t i=0;i<fNum;i++)
|
---|
240 | {
|
---|
241 | MHMcRate &rate = *GetRate(i);
|
---|
242 |
|
---|
243 | if (fTrigNSB)
|
---|
244 | rate.SetBackground(fTrigNSB[i], fSimNSB);
|
---|
245 | else
|
---|
246 | rate.SetBackground(0., fSimNSB);
|
---|
247 | }
|
---|
248 |
|
---|
249 | return kTRUE;
|
---|
250 | }
|
---|
251 |
|
---|
252 | // --------------------------------------------------------------------------
|
---|
253 | //
|
---|
254 | // The Process-function counts the number of simulated showers, the
|
---|
255 | // number of analised showers and the number of triggers. It also updates
|
---|
256 | // the limits for theta, phi, energy and impact parameter in the
|
---|
257 | // MHMcRate container.
|
---|
258 | //
|
---|
259 | Bool_t MMcTriggerRateCalc::Process()
|
---|
260 | {
|
---|
261 | //
|
---|
262 | // Counting analysed showers (except in the case in which the file
|
---|
263 | // contains only events that triggered, since then we do not know
|
---|
264 | // how many showers were analysed).
|
---|
265 | //
|
---|
266 |
|
---|
267 | if ( ! fMcRunHeader->GetAllEvtsTriggered() &&
|
---|
268 | fMcEvt->GetPhotElfromShower() )
|
---|
269 | fAnalShow++;
|
---|
270 |
|
---|
271 | //
|
---|
272 | // Set PartId and check it is the same for all events
|
---|
273 | //
|
---|
274 | if (fPartId < 0)
|
---|
275 | fPartId = fMcEvt->GetPartId();
|
---|
276 | else if (fPartId != fMcEvt->GetPartId())
|
---|
277 | {
|
---|
278 | *fLog << err << dbginf << "Incident particle type seems to change";
|
---|
279 | *fLog << "from " << fPartId << " to " << fMcEvt->GetPartId() << endl;
|
---|
280 | *fLog << "in input data files... aborting." << endl;
|
---|
281 | return kFALSE;
|
---|
282 | }
|
---|
283 |
|
---|
284 | //
|
---|
285 | // Getting angles, energy and impact parameter to set boundaries
|
---|
286 | //
|
---|
287 | const Float_t theta = fMcEvt->GetTheta();
|
---|
288 | const Float_t phi = fMcEvt->GetPhi();
|
---|
289 | const Float_t param = fMcEvt->GetImpact();
|
---|
290 | const Float_t ener = fMcEvt->GetEnergy()/1000.0;
|
---|
291 |
|
---|
292 | //
|
---|
293 | // Counting number of triggers
|
---|
294 | //
|
---|
295 | for (UInt_t i=0; i<fNum; i++)
|
---|
296 | {
|
---|
297 | if (GetTrig(i)->GetFirstLevel())
|
---|
298 | fTrigger[i] ++;
|
---|
299 |
|
---|
300 | if ( ! fMcRunHeader->GetAllEvtsTriggered())
|
---|
301 | GetRate(i)->UpdateBoundaries(ener, theta, phi, param);
|
---|
302 | }
|
---|
303 |
|
---|
304 | return kTRUE;
|
---|
305 | }
|
---|
306 |
|
---|
307 | // --------------------------------------------------------------------------
|
---|
308 | //
|
---|
309 | // The PostProcess-function calculates and shows the trigger rate
|
---|
310 | //
|
---|
311 | Bool_t MMcTriggerRateCalc::PostProcess()
|
---|
312 | {
|
---|
313 | for (UInt_t i=0; i<fNum; i++)
|
---|
314 | {
|
---|
315 | MHMcRate &rate = *GetRate(i);
|
---|
316 |
|
---|
317 | rate.SetParticle(fPartId);
|
---|
318 | switch (fPartId)
|
---|
319 | {
|
---|
320 | case kPROTON:
|
---|
321 | if (fMcCorRunHeader->GetSlopeSpec() != -2.75)
|
---|
322 | {
|
---|
323 | *fLog << err << dbginf <<
|
---|
324 | "Spectrum slope as read from input file ("<<
|
---|
325 | fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) for protons" << endl << "... aborting." << endl;
|
---|
326 | return kFALSE;
|
---|
327 | }
|
---|
328 | rate.SetFlux(0.1091, 2.75);
|
---|
329 | break;
|
---|
330 | case kHELIUM:
|
---|
331 | if (fMcCorRunHeader->GetSlopeSpec() != -2.62)
|
---|
332 | {
|
---|
333 | *fLog << err << dbginf <<
|
---|
334 | "Spectrum slope as read from input file ("<<
|
---|
335 | fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) for Helium" << endl << "... aborting." << endl;
|
---|
336 | return kFALSE;
|
---|
337 | }
|
---|
338 | rate.SetFlux(0.0660, 2.62);
|
---|
339 | break;
|
---|
340 | default:
|
---|
341 | *fLog << err << dbginf << "Unknown incident flux parameters for ";
|
---|
342 | *fLog << fPartId<< " particle Id ... aborting." << endl;
|
---|
343 | return kFALSE;
|
---|
344 | }
|
---|
345 | }
|
---|
346 |
|
---|
347 | //
|
---|
348 | // Computing trigger rate
|
---|
349 | //
|
---|
350 | for (UInt_t i=0; i<fNum; i++)
|
---|
351 | GetRate(i)->CalcRate(fTrigger[i], fAnalShow, fShowers);
|
---|
352 |
|
---|
353 | return kTRUE;
|
---|
354 | }
|
---|
355 |
|
---|
356 |
|
---|
357 | // Draw rate as a funtion of discriminator threshold.
|
---|
358 |
|
---|
359 | void MMcTriggerRateCalc::Draw()
|
---|
360 | {
|
---|
361 | TCanvas *c = MH::MakeDefCanvas("Rate");
|
---|
362 |
|
---|
363 | Float_t xmin = GetRate(0)->GetMeanThreshold()-0.55;
|
---|
364 | Float_t xmax = GetRate(fNum-1)->GetMeanThreshold()+0.55;
|
---|
365 | Int_t nbins = 10*(xmax-xmin);
|
---|
366 |
|
---|
367 | fHist[1] = new TH1F("Rate2","Trigger rate, mult. 2", nbins, xmin, xmax);
|
---|
368 | fHist[2] = new TH1F("Rate3","Trigger rate, mult. 3", nbins, xmin, xmax);
|
---|
369 | fHist[3] = new TH1F("Rate4","Trigger rate, mult. 4", nbins, xmin, xmax);
|
---|
370 | fHist[4] = new TH1F("Rate5","Trigger rate, mult. 5", nbins, xmin, xmax);
|
---|
371 |
|
---|
372 | for (UInt_t i=0; i<fNum; i++)
|
---|
373 | {
|
---|
374 | Short_t mult = GetRate(i)->GetMultiplicity();
|
---|
375 |
|
---|
376 | fHist[mult-1]->SetBinContent(fHist[mult-1]->FindBin(GetRate(i)->GetMeanThreshold()), GetRate(i)->GetTriggerRate());
|
---|
377 |
|
---|
378 | fHist[mult-1]->SetBinError(fHist[mult-1]->FindBin(GetRate(i)->GetMeanThreshold()), GetRate(i)->GetTriggerRateError());
|
---|
379 | }
|
---|
380 |
|
---|
381 | for (Int_t i = 1; i <=4; i++)
|
---|
382 | {
|
---|
383 | fHist[i]->SetLineWidth(2);
|
---|
384 | fHist[i]->SetMarkerStyle(20);
|
---|
385 | fHist[i]->SetMarkerSize(.5);
|
---|
386 | fHist[i]->SetMaximum(1.2*GetRate(0)->GetTriggerRate());
|
---|
387 | fHist[i]->SetMinimum(0.5*GetRate(fNum-1)->GetTriggerRate());
|
---|
388 | }
|
---|
389 |
|
---|
390 | fHist[2]->SetLineColor(1);
|
---|
391 | fHist[2]->SetMarkerColor(1);
|
---|
392 | fHist[2]->Draw();
|
---|
393 |
|
---|
394 | fHist[3]->SetLineColor(3);
|
---|
395 | fHist[3]->SetMarkerColor(3);
|
---|
396 | fHist[3]->Draw("same");
|
---|
397 |
|
---|
398 | fHist[4]->SetLineColor(4);
|
---|
399 | fHist[4]->Draw("same");
|
---|
400 | fHist[4]->SetMarkerColor(4);
|
---|
401 |
|
---|
402 | c->SetLogy();
|
---|
403 | c->SetGridy();
|
---|
404 | c->SetGridx();
|
---|
405 |
|
---|
406 | return;
|
---|
407 | }
|
---|
408 |
|
---|