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

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