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

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