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

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