source: trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc@ 2180

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