source: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc@ 9344

Last change on this file since 9344 was 9301, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 13.3 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-2008
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MHCollectionArea
29//
30// Class Version 2:
31// ----------------
32// + added //! to fMcEvt which was missing before
33// - removed obsolete data member fEnergy
34//
35//////////////////////////////////////////////////////////////////////////////
36#include "MHCollectionArea.h"
37
38#include <TMath.h>
39
40#include <TLatex.h>
41#include <TCanvas.h>
42#include <TPaveStats.h>
43
44#include "MLog.h"
45#include "MLogManip.h"
46
47#include "MString.h"
48#include "MBinning.h"
49
50#include "MMcEvt.hxx"
51#include "MMcRunHeader.hxx"
52#include "MMcCorsikaRunHeader.h"
53
54#include "MParList.h"
55#include "MParameters.h"
56
57ClassImp(MHCollectionArea);
58
59using namespace std;
60
61// --------------------------------------------------------------------------
62//
63// Creates the three necessary histograms:
64// - selected showers (input)
65// - all showers (input)
66// - collection area (result)
67//
68MHCollectionArea::MHCollectionArea(const char *name, const char *title)
69 : fMcEvt(0), fMcAreaRadius(-1), fIsExtern(kFALSE)
70{
71 // initialize the histogram for the distribution r vs E
72 //
73 // we set the energy range from 2 Gev to 20000 GeV (in log 4 orders
74 // of magnitude) and for each order we take 25 subdivision --> 100 xbins
75 //
76 // we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins
77 //
78 fName = name ? name : "MHCollectionArea";
79 fTitle = title ? title : "Collection Area vs. Energy/Theta";
80
81 fHistSel.SetName("SelEvts");
82 fHistSel.SetTitle("Number of Events after cuts");
83 fHistSel.SetXTitle("\\Theta [deg]");
84 fHistSel.SetYTitle("E [GeV]");
85 fHistSel.SetDirectory(NULL);
86 fHistSel.UseCurrentStyle();
87 fHistSel.SetLineColor(kBlue);
88
89 fHistAll.SetName("AllEvts");
90 fHistAll.SetTitle("Number of events produced");
91 fHistAll.SetXTitle("\\Theta [deg]");
92 fHistAll.SetYTitle("E_{mc} [GeV]");
93 fHistAll.SetDirectory(NULL);
94 fHistAll.UseCurrentStyle();
95
96 fHEnergy.SetName("CollEnergy");
97 fHEnergy.SetTitle("Collection Area");
98 fHEnergy.SetXTitle("E [GeV]");
99 fHEnergy.SetYTitle("A_{eff} [m^{2}]");
100 fHEnergy.SetDirectory(NULL);
101 fHEnergy.UseCurrentStyle();
102
103 MBinning binsa, binse, binst;
104 binse.SetEdgesLog(21, 6.3, 100000);
105 binst.SetEdgesASin(67, -0.005, 0.665);
106
107 binse.Apply(fHEnergy);
108
109 MH::SetBinning(&fHistSel, &binst, &binse);
110 MH::SetBinning(&fHistAll, &binst, &binse);
111
112 // For some unknown reasons this must be called after
113 // the binning has been initialized at least once
114 fHistSel.Sumw2();
115 fHistAll.Sumw2();
116 fHEnergy.Sumw2();
117}
118
119// --------------------------------------------------------------------------
120//
121// Return the Area defined by fMcAreaRadius
122//
123Double_t MHCollectionArea::GetCollectionAreaAbs() const
124{
125 return TMath::Pi()*fMcAreaRadius*fMcAreaRadius;
126}
127
128// --------------------------------------------------------------------------
129//
130// Calculate the Efficiency (collection area) and set the 'ReadyToSave'
131// flag
132//
133void MHCollectionArea::CalcEfficiency()
134{
135 TH1D *hsel = fHistSel.ProjectionY("Spy", -1, -1, "E");;
136 TH1D *hall = fHistAll.ProjectionY("Apy", -1, -1, "E");
137
138 //
139 // Impact parameter range.
140 //
141 const Float_t totalarea = GetCollectionAreaAbs();//TMath::Pi() * (r2*r2 - r1*r1);
142
143 // "b" option: calculate binomial errors
144 // Do not use totalarea inside the binomial error calculation:
145 // it is not a weight.
146 fHEnergy.Divide(hsel, hall, 1, 1, "b");
147#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
148 MH::SetBinomialErrors(fHEnergy, *hsel, *hall);
149#endif
150
151 fHEnergy.Scale(totalarea);
152
153 delete hsel;
154 delete hall;
155}
156
157
158Bool_t MHCollectionArea::SetupFill(const MParList *pl)
159{
160 fHistSel.Reset();
161 if (!fIsExtern)
162 fHistAll.Reset();
163
164 fMcEvt = (MMcEvt*)pl->FindObject("MMcEvt");
165 if (!fMcEvt)
166 {
167 *fLog << err << "MMcEvt not found... abort." << endl;
168 return kFALSE;
169 }
170 /*
171 fEnergy = (MParameterD*)pl->FindObject("MEnergyEst", "MParameterD");
172 if (!fEnergy)
173 {
174 *fLog << err << "MEnergyEst [MParameterD] not found... abort." << endl;
175 return kFALSE;
176 }
177 */
178 MBinning binst, binse;
179 binst.SetEdges(fHistAll, 'x');
180 binse.SetEdges(fHistAll, 'y');
181
182 //if (!fIsExtern)
183 {
184 MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning");
185 if (bins)
186 binst.SetEdges(*bins);
187
188 bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
189 if (bins)
190 binse.SetEdges(*bins);
191 }
192
193 binse.Apply(fHEnergy);
194
195 MH::SetBinning(&fHistSel, &binst, &binse);
196 MH::SetBinning(&fHistAll, &binst, &binse);
197
198 fMcAreaRadius = -1;
199 fCorsikaVersion = 0;
200
201 return kTRUE;
202}
203
204Bool_t MHCollectionArea::ReInit(MParList *plist)
205{
206 MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
207 if (!runheader)
208 {
209 *fLog << err << "MMcRunHeader not found... abort." << endl;
210 return kFALSE;
211 }
212
213 if (runheader->GetImpactMax()>fMcAreaRadius*100)
214 {
215 fMcAreaRadius = 0.01*runheader->GetImpactMax(); // cm->m
216 *fLog << inf << "Maximum simulated impact: " << fMcAreaRadius << "m" << endl;
217 }
218
219 if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
220 {
221 *fLog << warn;
222 *fLog << "Warning - Read files have different Corsika versions..." << endl;
223 *fLog << " Last file=" << fCorsikaVersion << " New file=" << runheader->GetCorsikaVersion() << endl;
224 }
225 fCorsikaVersion = runheader->GetCorsikaVersion();
226
227 if (fIsExtern)
228 return kTRUE;
229
230 fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
231 *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;
232
233 fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
234 *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;
235
236 MMcCorsikaRunHeader *crh = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
237 if (!crh)
238 {
239 *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl;
240 return kFALSE;
241 }
242
243 //
244 // Calculate approximately the original number of events in each
245 // energy bin:
246 //
247 const Float_t emin = crh->GetELowLim();
248 const Float_t emax = crh->GetEUppLim();
249 const Float_t expo = 1 + crh->GetSlopeSpec();
250 const Float_t k = runheader->GetNumSimulatedShowers() /
251 (pow(emax,expo) - pow(emin,expo));
252
253 const Int_t nbiny = fHistAll.GetNbinsY();
254
255 TAxis &axe = *fHistAll.GetYaxis();
256 for (Int_t i = 1; i <= nbiny; i++)
257 {
258 const Float_t e1 = axe.GetBinLowEdge(i);
259 const Float_t e2 = axe.GetBinLowEdge(i+1);
260
261 if (e1 < emin || e2 > emax)
262 continue;
263
264 const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
265 //
266 // We fill the i-th energy bin, with the total number of events
267 // Second argument of Fill would be impact parameter of each
268 // event, but we don't really need it for the collection area,
269 // so we just put a dummy value (1.)
270 //
271
272 const Float_t energy = (e1+e2)/2.;
273 fHistAll.Fill(20, energy, events);
274 // you have MMcRunHeader.fShowerThetaMin and MMcRunHeader.fShowerThetaMax
275 }
276
277 return kTRUE;
278}
279
280void MHCollectionArea::Paint(Option_t *option)
281{
282 if (TString(option)=="paint3")
283 {
284 /*
285 TH1 *h = dynamic_cast<TH1*>(gPad->FindObject("Efficiency"));
286 if (h)
287 {
288 const TString txt = Form("N/N_{0}=%.2f",
289 GetCollectionAreaEff(),
290 GetCollectionAreaAbs(), fMcAreaRadius);
291
292 TLatex text(0.31, 0.95, txt);
293 text.SetBit(TLatex::kTextNDC);
294 text.SetTextSize(0.04);
295 text.Paint();*/
296 return;
297 }
298 if (TString(option)=="paint4")
299 {
300 //const TString txt = Form("A_{eff}=%.0fm^{2} A_{abs}=%.0fm^{2} r=%.0fm",
301 // GetCollectionAreaEff(),
302 // GetCollectionAreaAbs(), fMcAreaRadius);
303 const TString txt = MString::Format("r_{max}=%.0fm --> A_{max}=%.0fm^{2}",
304 fMcAreaRadius, GetCollectionAreaAbs());
305
306 TLatex text(0.31, 0.95, txt);
307 text.SetBit(TLatex::kTextNDC);
308 text.SetTextSize(0.04);
309 text.Paint();
310 return;
311 }
312
313 TVirtualPad *pad;
314
315 TPaveStats *st=0;
316 for (int x=0; x<4; x++)
317 {
318 pad=gPad->GetPad(x+1);
319 if (!pad || !(st = (TPaveStats*)pad->GetPrimitive("stats")))
320 continue;
321
322 if (st->GetOptStat()==11)
323 continue;
324
325 const Double_t y1 = st->GetY1NDC();
326 const Double_t y2 = st->GetY2NDC();
327 const Double_t x1 = st->GetX1NDC();
328 const Double_t x2 = st->GetX2NDC();
329
330 st->SetY1NDC((y2-y1)/3+y1);
331 st->SetX1NDC((x2-x1)/3+x1);
332 st->SetOptStat(11);
333 }
334
335 pad = gPad;
336
337 TH1 *h1=0, *h2=0;
338
339 pad->cd(1);
340 if (gPad->FindObject("ProjSelX"))
341 fHistSel.ProjectionX("ProjSelX", -1, -1, "E");
342
343 pad->cd(2);
344 if (gPad->FindObject("ProjAllY"))
345 h1=fHistAll.ProjectionY("ProjAllY", -1, -1, "E");
346 if (gPad->FindObject("ProjSelY"))
347 h2=fHistSel.ProjectionY("ProjSelY", -1, -1, "E");
348
349 if (h1 && h1->GetMaximum()>0)
350 {
351 gPad->SetLogx();
352 gPad->SetLogy();
353 }
354
355 pad->cd(3);
356 TH1 *h=dynamic_cast<TH1*>(gPad->FindObject("Efficiency"));
357 if (h1 && h2 && h)
358 {
359 h->Divide(h2, h1, 1, 1, "b");
360#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
361 MH::SetBinomialErrors(*h, *h2, *h1);
362#endif
363 h->SetMinimum(0);
364 }
365
366 pad->cd(4);
367 CalcEfficiency();
368 if (fHEnergy.GetMaximum()>0)
369 {
370 gPad->SetLogx();
371 gPad->SetLogy();
372 }
373}
374
375void MHCollectionArea::Draw(Option_t *option)
376{
377 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
378
379 // Do the projection before painting the histograms into
380 // the individual pads
381 AppendPad();
382
383 pad->SetBorderMode(0);
384 pad->Divide(2,2);
385
386 TH1 *h=0, *h1=0, *h2=0;
387
388 if (fHistSel.GetNbinsX()>1)
389 {
390 pad->cd(1);
391 gPad->SetBorderMode(0);
392 gPad->SetGridx();
393 gPad->SetGridy();
394 /*
395 h = fHistAll.ProjectionX("ProjAllX", -1, -1, "E");
396 h->SetXTitle("\\Theta [\\circ]");
397 h->SetDirectory(NULL);
398 h->SetLineColor(kGreen);
399 h->SetBit(kCanDelete);
400 h->Draw();
401 */
402 h = fHistSel.ProjectionX("ProjSelX", -1, -1, "E");
403 h->SetXTitle("\\Theta [\\circ]");
404 h->SetDirectory(NULL);
405 h->SetLineColor(kRed);
406 h->SetBit(kCanDelete);
407 h->Draw("hist"/*"same"*/);
408 }
409 else
410 delete pad->GetPad(1);
411
412 if (fHistSel.GetNbinsY()>1)
413 {
414 pad->cd(2);
415 gPad->SetBorderMode(0);
416 gPad->SetGridx();
417 gPad->SetGridy();
418
419 h1 = fHistAll.ProjectionY("ProjAllY", -1, -1, "E");
420 h1->SetDirectory(NULL);
421 h1->SetLineColor(kGreen);
422 h1->SetXTitle("E [GeV]");
423 h1->SetBit(kCanDelete);
424 h1->Draw();
425
426 h2 = fHistSel.ProjectionY("ProjSelY", -1, -1, "E");
427 h2->SetDirectory(NULL);
428 h2->SetLineColor(kRed);
429 h2->SetBit(kCanDelete);
430 h2->Draw("same");
431 }
432 else
433 delete pad->GetPad(2);
434
435 if (h1 && h2)
436 {
437 pad->cd(3);
438 gPad->SetBorderMode(0);
439 gPad->SetGridx();
440 gPad->SetGridy();
441 gPad->SetLogx();
442 h = h2->DrawCopy();
443 h->Divide(h2, h1, 1, 1, "b");
444#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
445 MH::SetBinomialErrors(*h, *h2, *h1);
446#endif
447 h->SetNameTitle("Efficiency", "Combined cut and trigger efficiency");
448 h->SetDirectory(NULL);
449 AppendPad("paint3");
450 }
451 else
452 delete pad->GetPad(4);
453
454 if (fHEnergy.GetNbinsX()>1)
455 {
456 pad->cd(4);
457 gPad->SetBorderMode(0);
458 gPad->SetGridx();
459 gPad->SetGridy();
460 fHEnergy.Draw();
461 AppendPad("paint4");
462 }
463 else
464 delete pad->GetPad(4);
465}
466
467Int_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight)
468{
469 const Double_t energy = fMcEvt->GetEnergy();
470 const Double_t theta = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
471
472 fHistSel.Fill(theta, energy, weight);
473
474 return kTRUE;
475}
476
477Bool_t MHCollectionArea::Finalize()
478{
479 *fLog << all << "Maximum simulated impact found: " << fMcAreaRadius << "m" << endl;
480
481 CalcEfficiency();
482
483 return kTRUE;
484}
Note: See TracBrowser for help on using the repository browser.