source: branches/Mars_MC/mhistmc/MHMcCollectionArea.cc@ 17054

Last change on this file since 17054 was 6539, checked in by moralejo, 20 years ago
*** empty log message ***
File size: 11.6 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 ! Author(s): Abelardo Moralejo 2/2005 <mailto:moralejo@pd.infn.it>
21 !
22 ! Copyright: MAGIC Software Development, 2000-2005
23 !
24 !
25 \* ======================================================================== */
26
27//////////////////////////////////////////////////////////////////////////////
28// //
29// MHMcCollectionArea //
30// //
31//////////////////////////////////////////////////////////////////////////////
32
33#include "MHMcCollectionArea.h"
34
35#include <TH2.h>
36#include <TH3.h>
37#include <TCanvas.h>
38#include <THStack.h>
39#include <TLegend.h>
40#include <TArrayD.h>
41
42#include "MH.h"
43#include "MBinning.h"
44
45#include "MLog.h"
46#include "MLogManip.h"
47
48
49ClassImp(MHMcCollectionArea);
50
51using namespace std;
52
53////////////////////////////////////////////////////////////////////////////////
54//
55// Constructor. Creates the three necessary histograms:
56// - selected showers (input)
57// - all showers (input)
58// - collection area (result)
59//
60MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title):
61 fImpactBins(50), fImpactMax(500.), fMinEvents(10)
62{
63 fName = name ? name : "MHMcCollectionArea";
64 fTitle = title ? title : "Collection Area vs. Theta vs. Energy";
65
66 //
67 // Initialize the histogram for the distribution
68 // Theta vs impact parameter vs E (z, y, x)
69 //
70 // As default we set the energy range from 2 Gev to 20000 GeV (in log 4
71 // orders of magnitude) and for each order we take 25 subdivisions -->
72 // 100 xbins. We set the radius range from 0 m to 500 m with 10 m bin -->
73 // 50 ybins. We make bins equally spaced in cos(theta)
74 //
75 // The coarse binning (of fHistColCoarse) is not set by default, the
76 // PreProcess of mmc/MMcCollectionAreaCalc will do it with the binnings
77 // found in the parameter list.
78 //
79
80 MBinning binsx;
81 MBinning binsy;
82 MBinning binsz;
83
84 Int_t nbins = 32;
85 TArrayD edges(nbins+1);
86
87 edges[0] = 0;
88
89 for(int i = 0; i < nbins; i++)
90 {
91 Double_t x = 1 - i*0.01; // x = cos(theta)
92 edges[i+1] = acos(x-0.005)*kRad2Deg;
93 }
94
95 binsx.SetEdgesLog(100, 2., 20000); // Energy [GeV]
96 binsy.SetEdges (fImpactBins, 0, fImpactMax); // Impact parameter [m]
97 binsz.SetEdges (edges); // Theta [deg]
98
99 fHistAll = new TH3D();
100 fHistSel = new TH3D();
101 fHistCol = new TH2D();
102 fHistColCoarse = new TH2D();
103
104 MH::SetBinning(fHistAll, &binsx, &binsy, &binsz);
105 MH::SetBinning(fHistSel, &binsx, &binsy, &binsz);
106
107 fHistColCoarse->SetName("CollectionArea");
108 fHistCol->SetName("CollAreaFineBins");
109 fHistAll->SetName("AllEvents");
110 fHistSel->SetName("SelectedEvents");
111
112 fHistAll->Sumw2();
113 fHistSel->Sumw2();
114
115 fHistColCoarse->SetTitle(fTitle);
116 fHistCol->SetTitle(fTitle);
117 fHistAll->SetTitle("All showers - Theta vs Radius vs Energy");
118 fHistSel->SetTitle("Selected showers - Theta vs Radius vs Energy");
119
120 fHistAll->SetDirectory(NULL);
121 fHistSel->SetDirectory(NULL);
122 fHistCol->SetDirectory(NULL);
123 fHistColCoarse->SetDirectory(NULL);
124
125 fHistAll->UseCurrentStyle();
126 fHistSel->UseCurrentStyle();
127 fHistCol->UseCurrentStyle();
128 fHistColCoarse->UseCurrentStyle();
129
130 fHistAll->SetXTitle("E [GeV]");
131 fHistAll->SetYTitle("r [m]");
132 fHistAll->SetZTitle("\\theta [\\circ]");
133
134 fHistSel->SetXTitle("E [GeV]");
135 fHistSel->SetYTitle("r [m]");
136 fHistSel->SetZTitle("\\theta [\\circ]");
137
138 fHistCol->SetXTitle("E [GeV]");
139 fHistCol->SetYTitle("\\theta [\\circ]");
140 fHistCol->SetZTitle("A [m^{2}]");
141
142 fHistColCoarse->SetXTitle("E [GeV]");
143 fHistColCoarse->SetYTitle("\\theta [\\circ]");
144 fHistColCoarse->SetZTitle("A [m^{2}]");
145}
146
147// --------------------------------------------------------------------------
148//
149// Delete the three histograms
150//
151MHMcCollectionArea::~MHMcCollectionArea()
152{
153 delete fHistAll;
154 delete fHistSel;
155 delete fHistCol;
156}
157
158// --------------------------------------------------------------------------
159//
160// Set the (fine) binnings of histograms fHistAll, fHistSel used in the
161// calculations. We do not need to change impact parameter binning.
162//
163void MHMcCollectionArea::SetBinnings(const MBinning &binsEnergy,
164 const MBinning &binsTheta)
165{
166 MBinning binsImpact;
167 binsImpact.SetEdges(fImpactBins, 0., fImpactMax);
168
169 MH::SetBinning(fHistAll, &binsEnergy, &binsImpact, &binsTheta);
170 MH::SetBinning(fHistSel, &binsEnergy, &binsImpact, &binsTheta);
171
172 fHistAll->Sumw2();
173 fHistSel->Sumw2();
174}
175
176
177// --------------------------------------------------------------------------
178//
179// Set the binnings of the histogram fHistColCoarse, the effective areas
180// in the coarse bins used in the analysis.
181//
182//
183void MHMcCollectionArea::SetCoarseBinnings(const MBinning &binsEnergy,
184 const MBinning &binsTheta)
185{
186 MH::SetBinning(fHistColCoarse, &binsEnergy, &binsTheta);
187}
188
189// --------------------------------------------------------------------------
190//
191// Fill data into the histogram which contains all showers
192//
193void MHMcCollectionArea::FillAll(Double_t energy, Double_t radius, Double_t theta)
194{
195 fHistAll->Fill(energy, radius, theta);
196}
197
198// --------------------------------------------------------------------------
199//
200// Fill data into the histogram which contains the selected showers
201//
202void MHMcCollectionArea::FillSel(Double_t energy, Double_t radius, Double_t theta)
203{
204 fHistSel->Fill(energy, radius, theta);
205}
206
207// --------------------------------------------------------------------------
208//
209// Draw
210//
211void MHMcCollectionArea::Draw(Option_t* option)
212{
213 //
214 // Lego plot
215 //
216 TCanvas *c1 = new TCanvas();
217 c1->SetLogx();
218 c1->SetLogz();
219 c1->SetGridx();
220 c1->SetGridy();
221
222 fHistCol->Draw("lego2");
223
224 //
225 // Averagye Aeff
226 //
227 TCanvas *c2 = new TCanvas();
228 c2->SetLogx();
229 c2->SetLogy();
230 c2->SetGridx();
231 c2->SetGridy();
232
233 TH1D* harea = fHistCol->ProjectionX();
234 harea->Draw("e1");
235
236 //
237 // Plot the Aeff for the different theta
238 //
239 TCanvas *c3 = new TCanvas();
240 c3->SetLogx();
241 c3->SetLogy();
242 c3->SetGridx();
243 c3->SetGridy();
244
245 TLegend * leg = new TLegend(0.73,0.65,0.89,0.89);
246
247 TAxis* yaxis = fHistCol->GetYaxis();
248 const Int_t nbiny = fHistCol->GetYaxis()->GetNbins();
249
250 THStack* hs = new THStack("aa","aa");
251
252 hs->Add(harea,"e1");
253 leg->AddEntry(harea,"All","l");
254
255 for(Int_t iy=1; iy<=nbiny; iy++)
256 {
257
258 TH1D* h1= fHistCol->ProjectionX(Form("%d",iy),iy,iy);
259
260 if(h1->GetEntries()==0)
261 continue;
262
263 cout <<h1->GetEntries() << endl;
264
265 leg->AddEntry(h1,Form("\\theta = %.0f",yaxis->GetBinCenter(iy)),"l");
266 h1->SetLineColor(iy);
267 hs->Add(h1,"e1");
268 }
269 hs->SetMinimum(1);
270
271 hs->Draw("nostack");
272 leg->Draw();
273
274}
275
276// --------------------------------------------------------------------------
277//
278// Calculate the collection area and set the 'ReadyToSave' flag
279// We first calculate the area in fine energy bins, and then do a
280// weighted mean to obtain the area in coarse bins. The weights in
281// the coarse bins are intended to account for the effect of the
282// energy spectrum in the effective area itself. The weights
283// are taken from the tentative differential spectrum dN_gam/dE given
284// through the function "spectrum". If no such function is supplied,
285// then no weights are applied (and hence the spectrum will be as a
286// flat spectrum in dN_gam/dE). Of course we have a "generated" MC
287// spectrum, but within each fine bin the differences in spectrum
288// should not change the result (if bins are fine enough). With no
289// supplied tentative spectrum, each fine bin is weighted equally in
290// calculating the area in the coarse bin, and so it is like having a
291// flat spectrum.
292//
293// You can run this Calc procedure on an already existing
294// MHMcCollectionArea object, as long as it is filled.
295//
296void MHMcCollectionArea::Calc(TF1 *spectrum)
297{
298 // Search last impact parameter bin containing events
299 // FIXME: this should be done independently for each theta angle.
300 //
301 TH1D &himpact = *(TH1D*)fHistAll->Project3D("y");
302
303 Int_t impbin;
304 for (impbin = himpact.GetNbinsX(); impbin > 0; impbin--)
305 if (himpact.GetBinContent(impbin)>0)
306 break;
307
308 Float_t max_radius = himpact.GetBinLowEdge(impbin);
309
310 Float_t total_area = TMath::Pi()*max_radius*max_radius;
311
312 for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
313 for (Int_t iz = 1; iz <= fHistAll->GetNbinsZ(); iz++)
314 {
315 fHistAll->SetBinContent(ix, impbin, iz, 0.);
316 fHistSel->SetBinContent(ix, impbin, iz, 0.);
317 fHistAll->SetBinError(ix, impbin, iz, 0.);
318 fHistSel->SetBinError(ix, impbin, iz, 0.);
319 }
320
321 TH2D &histsel = *(TH2D*)fHistSel->Project3D("zx,e");
322 TH2D &histall = *(TH2D*)fHistAll->Project3D("zx,e");
323 // "e" option means that errors are computed!
324
325
326 TAxis &xaxis = *histsel.GetXaxis();
327 TAxis &yaxis = *histsel.GetYaxis();
328 MH::SetBinning(fHistCol, &xaxis, &yaxis);
329
330 cout << "Total considered MC area = pi * " << max_radius
331 << "^2 square meters" << endl;
332
333 fHistCol->Sumw2();
334 fHistCol->Divide(&histsel, &histall, total_area, 1., "b");
335
336 //
337 // Now get the effective area in the selected coarse bins. Weight
338 // the values in the small bins according the supplied tentative
339 // spectrum, if it has been supplied as argument of Calc.
340 //
341
342 for (Int_t ibin = 1; ibin <= fHistColCoarse->GetNbinsX(); ibin++)
343 for (Int_t jbin = 1; jbin <= fHistColCoarse->GetNbinsY(); jbin++)
344 {
345 Float_t maxenergy = fHistColCoarse->GetXaxis()->GetBinUpEdge(ibin);
346 Float_t minenergy = fHistColCoarse->GetXaxis()->GetBinLowEdge(ibin);
347
348 Float_t maxtheta = fHistColCoarse->GetYaxis()->GetBinUpEdge(jbin);
349 Float_t mintheta = fHistColCoarse->GetYaxis()->GetBinLowEdge(jbin);
350
351 // Fine bins ranges covered by the coarse bin ibin, jbin:
352 Int_t ibin2max = fHistCol->GetXaxis()->FindBin(maxenergy);
353 Int_t ibin2min = fHistCol->GetXaxis()->FindBin(minenergy);
354
355 Int_t jbin2max = fHistCol->GetYaxis()->FindBin(maxtheta);
356 Int_t jbin2min = fHistCol->GetYaxis()->FindBin(mintheta);
357
358 Float_t area = 0.;
359 Float_t errarea = 0.;
360 Float_t norm = 0;
361
362 for (Int_t ibin2 = ibin2min; ibin2 <= ibin2max; ibin2++)
363 {
364 Float_t weight = spectrum? spectrum->
365 Eval(fHistCol->GetXaxis()->GetBinCenter(ibin2)) : 1.;
366
367 for (Int_t jbin2 = jbin2min; jbin2 <= jbin2max; jbin2++)
368 {
369 // Skip bins with too few produced MC events
370 if (histall.GetBinContent(ibin2,jbin2) < fMinEvents)
371 continue;
372
373 area += weight * fHistCol->GetBinContent(ibin2,jbin2);
374 norm += weight;
375 errarea += pow(weight * fHistCol->
376 GetBinError(ibin2,jbin2), 2.);
377 }
378 }
379 if (norm > 0.)
380 {
381 area /= norm;
382 errarea = sqrt(errarea)/norm;
383 }
384
385 fHistColCoarse->SetBinContent(ibin, jbin, area);
386 fHistColCoarse->SetBinError(ibin, jbin, errarea);
387 }
388
389 SetReadyToSave();
390}
391
392
Note: See TracBrowser for help on using the repository browser.