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

Last change on this file since 1790 was 1790, checked in by moralejo, 22 years ago
*** empty log message ***
File size: 11.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@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
47ClassImp(MMcTriggerRateCalc);
48
49void 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
87Bool_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
167MMcTriggerRateCalc::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//
184MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, float *trigbg,float simbg,
185 const char *name, const char *title)
186{
187 Init(dim, trigbg, simbg, name, title);
188}
189
190MMcTriggerRateCalc::~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//
208Bool_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//
259Bool_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//
311Bool_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
359void 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
Note: See TracBrowser for help on using the repository browser.