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

Last change on this file since 1800 was 1800, checked in by moralejo, 22 years ago
*** empty log message ***
File size: 11.9 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 if (GetRate(i)->GetTriggerCondNum() > 0)
146 {
147 th += ";";
148 th += GetRate(i)->GetTriggerCondNum();
149 }
150
151 MMcTrigHeader* trighead = (MMcTrigHeader*) pList->FindObject(th);
152 GetRate(i)->SetMeanThreshold(trighead->GetMeanThreshold());
153 GetRate(i)->SetMultiplicity(trighead->GetMultiplicity());
154
155 }
156
157 return kTRUE;
158}
159
160// --------------------------------------------------------------------------
161//
162// overloaded constructor I
163//
164// dim: fDimension
165// *trigbg: number of NSB triggers for
166// a given trigger condition.
167// simbg: Number of simulated events for the NSB
168// rate: rate of incident showers
169//
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 (fMcCorRunHeader->GetSlopeSpec() != -2.75)
326 {
327 *fLog << err << dbginf <<
328 "Spectrum slope as read from input file ("<<
329 fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected 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 (fMcCorRunHeader->GetSlopeSpec() != -2.62)
336 {
337 *fLog << err << dbginf <<
338 "Spectrum slope as read from input file ("<<
339 fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) 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// Draw rate as a funtion of discriminator threshold.
362
363void MMcTriggerRateCalc::Draw()
364{
365 TCanvas *c = MH::MakeDefCanvas("Rate");
366
367 Float_t xmin = GetRate(0)->GetMeanThreshold()-0.55;
368 Float_t xmax = GetRate(fNum-1)->GetMeanThreshold()+0.55;
369 Int_t nbins = 10*(xmax-xmin);
370
371 fHist[1] = new TH1F("Rate2","Trigger rate, mult. 2", nbins, xmin, xmax);
372 fHist[2] = new TH1F("Rate3","Trigger rate, mult. 3", nbins, xmin, xmax);
373 fHist[3] = new TH1F("Rate4","Trigger rate, mult. 4", nbins, xmin, xmax);
374 fHist[4] = new TH1F("Rate5","Trigger rate, mult. 5", nbins, xmin, xmax);
375
376 for (UInt_t i=0; i<fNum; i++)
377 {
378 Short_t mult = GetRate(i)->GetMultiplicity();
379
380 fHist[mult-1]->SetBinContent(fHist[mult-1]->FindBin(GetRate(i)->GetMeanThreshold()), GetRate(i)->GetTriggerRate());
381
382 fHist[mult-1]->SetBinError(fHist[mult-1]->FindBin(GetRate(i)->GetMeanThreshold()), GetRate(i)->GetTriggerRateError());
383 }
384
385 for (Int_t i = 1; i <=4; i++)
386 {
387 fHist[i]->SetLineWidth(2);
388 fHist[i]->SetMarkerStyle(20);
389 fHist[i]->SetMarkerSize(.5);
390 }
391
392 c->DrawFrame (xmin, 0.5*GetRate(fNum-1)->GetTriggerRate(), xmax, 1.2*GetRate(0)->GetTriggerRate(), "Trigger rate for multiplicity = 3, 4, 5");
393
394 c->SetLogy();
395 c->SetGridy();
396 c->SetGridx();
397
398 fHist[2]->SetLineColor(1);
399 fHist[2]->SetMarkerColor(1);
400 fHist[2]->SetMinimum(0.5*GetRate(fNum-1)->GetTriggerRate());
401 fHist[2]->GetXaxis()->SetTitle("Discr. threshold (mV)");
402 fHist[2]->GetYaxis()->SetTitle("Trigger rate (Hz)");
403 fHist[2]->GetYaxis()->SetTitleOffset(1.2);
404 fHist[2]->Draw("axis");
405 fHist[2]->Draw("same");
406
407 fHist[3]->SetLineColor(3);
408 fHist[3]->SetMarkerColor(3);
409 fHist[3]->Draw("same");
410
411 fHist[4]->SetLineColor(4);
412 fHist[4]->Draw("same");
413 fHist[4]->SetMarkerColor(4);
414
415 return;
416}
417
Note: See TracBrowser for help on using the repository browser.