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

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