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

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