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

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