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

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