source: branches/removing_cpp11_features/mmontecarlo/MMcTriggerRateCalc.cc@ 18558

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