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

Last change on this file since 7008 was 6976, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 11.9 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 <TCanvas.h>
34#include <TPaveStats.h>
35
36#include "MLog.h"
37#include "MLogManip.h"
38
39#include "MH.h"
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 [GeV]");
84 fHistAll.SetDirectory(NULL);
85 fHistAll.UseCurrentStyle();
86
87 fHEnergy.SetName("CollEnergy");
88 fHEnergy.SetTitle("Collection Area vs Energy");
89 fHEnergy.SetXTitle("E [GeV]");
90 fHEnergy.SetYTitle("A [m^{2}]");
91 fHEnergy.SetDirectory(NULL);
92 fHEnergy.UseCurrentStyle();
93
94 MBinning binsa, binse, binst;
95 binse.SetEdgesLog(15, 10, 1000000);
96 binst.SetEdgesCos(50, 0, 60);
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 r1 = 0;
122 const Float_t r2 = fMcAreaRadius;
123
124 //*fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of " << fMCAreaRadius << " meters for the MC events. Check that this is the true one!" <<endl<<endl;
125 const Float_t total_area = TMath::Pi() * (r2*r2 - r1*r1);
126
127 for (Int_t ix=1; ix<=nbinx; ix++)
128 {
129 const Float_t Na = hall->GetBinContent(ix);
130
131 if (Na <= 0)
132 continue;
133
134 const Float_t Ns = hsel->GetBinContent(ix);
135
136 // Since Na is an estimate of the total number of showers generated
137 // in the energy bin, it may happen that Ns (triggered showers) is
138 // larger than Na. In that case, the bin is skipped:
139
140 if (Na < Ns)
141 continue;
142
143 const Double_t eff = Ns/Na;
144
145 const Double_t efferr = TMath::Sqrt((1.-eff)*Ns)/Na;
146
147 const Float_t col_area = eff * total_area;
148 const Float_t col_area_error = efferr * total_area;
149
150 fHEnergy.SetBinContent(ix, col_area);
151 fHEnergy.SetBinError(ix, col_area_error);
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("BinningEnergy", "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 TVirtualPad *pad;
278
279 TPaveStats *st=0;
280 for (int x=0; x<4; x++)
281 {
282 pad=gPad->GetPad(x+1);
283 if (!pad || !(st = (TPaveStats*)pad->GetPrimitive("stats")))
284 continue;
285
286 if (st->GetOptStat()==11)
287 continue;
288
289 const Double_t y1 = st->GetY1NDC();
290 const Double_t y2 = st->GetY2NDC();
291 const Double_t x1 = st->GetX1NDC();
292 const Double_t x2 = st->GetX2NDC();
293
294 st->SetY1NDC((y2-y1)/3+y1);
295 st->SetX1NDC((x2-x1)/3+x1);
296 st->SetOptStat(11);
297 }
298
299 pad = gPad;
300
301 TH1 *h1=0, *h2=0;
302
303 pad->cd(1);
304 if (gPad->FindObject("ProjSelX"))
305 fHistSel.ProjectionX("ProjSelX", -1, 9999, "E");
306
307 pad->cd(2);
308 if (gPad->FindObject("ProjAllY"))
309 h1=fHistAll.ProjectionY("ProjAllY", -1, 9999, "E");
310 if (gPad->FindObject("ProjSelY"))
311 h2=fHistSel.ProjectionY("ProjSelY", -1, 9999, "E");
312
313 if (h1 && h1->GetMaximum()>0)
314 {
315 gPad->SetLogx();
316 gPad->SetLogy();
317 }
318
319 pad->cd(3);
320 TH1 *h=dynamic_cast<TH1*>(gPad->FindObject("Efficiency"));
321 if (h1 && h2 && h)
322 {
323 h->Divide(h2, h1);
324 h->SetMinimum(0);
325 }
326
327 pad->cd(4);
328 CalcEfficiency();
329 if (fHEnergy.GetMaximum()>0)
330 {
331 gPad->SetLogx();
332 gPad->SetLogy();
333 }
334}
335
336void MHCollectionArea::Draw(Option_t *option)
337{
338 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
339
340 // Do the projection before painting the histograms into
341 // the individual pads
342 AppendPad();
343
344 pad->SetBorderMode(0);
345 pad->Divide(2,2);
346
347 TH1 *h=0, *h1=0, *h2=0;
348
349 if (fHistSel.GetNbinsX()>1)
350 {
351 pad->cd(1);
352 gPad->SetBorderMode(0);
353 gPad->SetGridx();
354 gPad->SetGridy();
355 /*
356 h = fHistAll.ProjectionX("ProjAllX", -1, 9999, "E");
357 h->SetXTitle("\\Theta [\\circ]");
358 h->SetDirectory(NULL);
359 h->SetLineColor(kGreen);
360 h->SetBit(kCanDelete);
361 h->Draw();
362 */
363 h = fHistSel.ProjectionX("ProjSelX", -1, 9999, "E");
364 h->SetXTitle("\\Theta [\\circ]");
365 h->SetDirectory(NULL);
366 h->SetLineColor(kRed);
367 h->SetBit(kCanDelete);
368 h->Draw("hist"/*"same"*/);
369 }
370 else
371 delete pad->GetPad(1);
372
373 if (fHistSel.GetNbinsY()>1)
374 {
375 pad->cd(2);
376 gPad->SetBorderMode(0);
377 gPad->SetGridx();
378 gPad->SetGridy();
379
380 h1 = fHistAll.ProjectionY("ProjAllY", -1, 9999, "E");
381 h1->SetDirectory(NULL);
382 h1->SetLineColor(kGreen);
383 h1->SetXTitle("E [GeV]");
384 h1->SetBit(kCanDelete);
385 h1->Draw();
386
387 h2 = fHistSel.ProjectionY("ProjSelY", -1, 9999, "E");
388 h2->SetDirectory(NULL);
389 h2->SetLineColor(kRed);
390 h2->SetBit(kCanDelete);
391 h2->Draw("same");
392 }
393 else
394 delete pad->GetPad(2);
395
396 if (h1 && h2)
397 {
398 pad->cd(3);
399 gPad->SetBorderMode(0);
400 gPad->SetGridx();
401 gPad->SetGridy();
402 gPad->SetLogx();
403 h = h2->DrawCopy();
404 h->Divide(h1);
405 h->SetNameTitle("Efficiency", "Combined cut and trigger efficiency");
406 h->SetDirectory(NULL);
407 }
408 else
409 delete pad->GetPad(4);
410
411 if (fHEnergy.GetNbinsX()>1)
412 {
413 pad->cd(4);
414 gPad->SetBorderMode(0);
415 gPad->SetGridx();
416 gPad->SetGridy();
417 fHEnergy.Draw();
418 }
419 else
420 delete pad->GetPad(4);
421}
422
423Bool_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight)
424{
425 const Double_t energy = fMcEvt->GetEnergy()/*fEnergy->GetVal()*/;
426 const Double_t theta = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
427
428 //*fLog << energy << " " << theta << endl;
429
430 fHistSel.Fill(theta, energy);
431
432 return kTRUE;
433}
Note: See TracBrowser for help on using the repository browser.