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

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