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

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