source: branches/MarsMoreSimulationTruth/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc

Last change on this file was 6486, checked in by rico, 20 years ago
*** empty log message ***
File size: 12.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! Author(s): Javier Rico 02/2005 <mailto:jrico@ifae.es>
18!
19! Copyright: MAGIC Software Development, 2000-2005
20!
21!
22\* ======================================================================== */
23
24//////////////////////////////////////////////////////////////////////////////
25//
26// Computes the Effective areas and coefficients for unfolding for a given
27// spectrum that can be parametrized by a function
28//
29//////////////////////////////////////////////////////////////////////////////
30
31#include <fstream>
32#include <math.h>
33
34#include "MEffAreaAndCoeffCalc.h"
35
36#include "TF1.h"
37#include "MHillas.h"
38#include "MMcEvt.hxx"
39#include "TH2F.h"
40#include "TFile.h"
41
42#include "MLog.h"
43#include "MLogManip.h"
44
45ClassImp(MEffAreaAndCoeffCalc);
46
47using namespace std;
48
49const Int_t ntbins=1; // default number of theta bins
50const Double_t tbin[ntbins+1] = {0,90}; // default theta bins bounds
51const Int_t nebins = 10; // default number of energy bins
52const Float_t emin = 10.; // default minimum energy value
53const Float_t emax = 10000.; // default maximum energy value
54const Int_t nsubbins = 20; // default number of subbins per bin
55const Char_t* deff = "4.e9*pow(x,-2.6+1)"; // default spectrum function
56
57// -------------------------------------------------------------------------
58//
59// Constructor
60//
61MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
62 : fSpec(NULL), fEmin(emin), fEmax(emax), fEbins(nebins), fEsubbins(nsubbins),
63 fWeight(NULL), fCoeff(NULL), fEffA(NULL)//, fFile(NULL)
64{
65 // set the default function
66 SetFunction(deff);
67
68 // create the TChains
69 fCini = new TChain("OriginalMC");
70 fCcut = new TChain("Events");
71
72 // define some useful aliases
73 fCini->SetAlias("logenergy","log10(MMcEvtBasic.fEnergy)");
74 fCini->SetAlias("theta","MMcEvtBasic.fTelescopeTheta*180./3.14159");
75
76 fCcut->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
77 fCcut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159");
78
79 fCcut->SetBranchAddress("MHillas.",&fHillas);
80 fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
81
82 // initial value of the zenith angle binning
83 SetThetaBinning(ntbins,tbin);
84
85 // borra
86 // fFile = new TFile("coeftest.root","RECREATE");
87}
88
89// -------------------------------------------------------------------------
90//
91// Destructor
92//
93MEffAreaAndCoeffCalc::~MEffAreaAndCoeffCalc()
94{
95 if(fSpec)
96 delete fSpec;
97
98 if(fTbin)
99 delete [] fTbin;
100
101 if(fWeight)
102 delete [] fWeight;
103
104 if(fCoeff)
105 delete fCoeff;
106
107 if(fEffA)
108 delete fEffA;
109
110
111 delete fCini;
112 delete fCcut;
113}
114
115// -------------------------------------------------------------------------
116//
117// Set the binning in zenith angle
118//
119void MEffAreaAndCoeffCalc::SetThetaBinning(Int_t n, const Double_t* binlist)
120{
121 fNTbins=n;
122 if(fTbin) delete [] fTbin;
123 fTbin = new Double_t[n+1];
124 for(Int_t i=0;i<n+1;i++)
125 fTbin[i] = binlist[i];
126}
127
128// -------------------------------------------------------------------------
129//
130// Set the function by expression and minimum and maximum energies
131//
132void MEffAreaAndCoeffCalc::SetFunction(const Char_t* chfunc, Float_t emin, Float_t emax)
133{
134 if(fSpec)
135 delete fSpec;
136 if(emin<=0 || emax<=0 || emax<emin)
137 {
138 emin = fEmin;
139 emax = fEmax;
140 }
141 else
142 {
143 fEmin = emin;
144 fEmax = emax;
145 }
146 fSpec = new TF1("fspec",chfunc,emin,emax);
147}
148
149// -------------------------------------------------------------------------
150//
151// Set the function by function pointer
152//
153void MEffAreaAndCoeffCalc::SetFunction(TF1* newf)
154{
155 if(fSpec)
156 delete fSpec;
157 fSpec = newf;
158}
159
160// -------------------------------------------------------------------------
161//
162// fill the histogram containing the original sample energy spectrum
163//
164void MEffAreaAndCoeffCalc::FillOriginalSpectrum()
165{
166}
167
168// -------------------------------------------------------------------------
169//
170// compute the weights for a particular input spectrum
171//
172void MEffAreaAndCoeffCalc::ComputeWeights()
173{
174 // OJO!! maybe this should be hard-coded somewhere else
175 const Float_t abslogmin=0.;
176 const Float_t abslogmax=5.;
177
178 const Float_t logemin = TMath::Log10(fEmin);
179 const Float_t logemax = TMath::Log10(fEmax);
180 const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
181 const Float_t desub = de/fEsubbins; // subbin size (in log)
182 const Int_t esbins = fEbins*fEsubbins; // total number of subbins
183
184 // compute the binning for weights histogram
185 const Int_t nmindist = (logemin>abslogmin)? Int_t((logemin-abslogmin)/desub) : 0;
186 const Int_t nmaxdist = (logemax<abslogmax)? Int_t((abslogmax-logemax)/desub) : 0;
187
188 fLogEWmin = logemin-desub*nmindist;
189 fLogEWmax = logemax+desub*nmaxdist;
190 fEWbins = nmindist+esbins+nmaxdist;
191
192 // reset the weights array
193 if(fWeight)
194 delete [] fWeight;
195 fWeight = new Double_t[fEWbins];
196
197
198 TH1F* horigs = new TH1F("horigs","Original energy spectrum",fEWbins,fLogEWmin,fLogEWmax);
199 fCini->Draw("logenergy>>horigs","","goff");
200 // borra
201 // horigs->Write();
202
203 // borra
204 // TH1F hw("hw","hw",fEWbins,fLogEWmin,fLogEWmax);
205 for(Int_t i=0;i<fEWbins;i++)
206 {
207 const Float_t denom = horigs->GetBinContent(i+1); // number of events in initial spectrum
208 const Float_t ew = TMath::Power(10,fLogEWmin+(i+0.5)*desub); // real energy
209 const Float_t numer = fSpec->Eval(ew); // number of events for the required spectrum
210 if(denom>10)
211 fWeight[i]=numer/denom;
212 else
213 {
214 // cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 " << endl;
215 fWeight[i]=-1;
216 }
217 // borra
218 // hw.SetBinContent(i+1,fWeight[i]);
219 }
220 // borra
221 // hw.Write();
222
223 delete horigs;
224}
225
226// --------------------------------------------------------------
227//
228// compute the coefficients used for the (iterative) unfolding
229//
230void MEffAreaAndCoeffCalc::ComputeCoefficients()
231{
232 if(!fWeight)
233 {
234 cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: No weights computed! nothing done" << endl;
235 return;
236 }
237
238 const Float_t logemin = TMath::Log10(fEmin);
239 const Float_t logemax = TMath::Log10(fEmax);
240 const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
241 const Float_t desub = de/fEsubbins; // subbin size (in log)
242 const Int_t nentries = Int_t(fCcut->GetEntries());
243
244 // declare needed histos
245 TH2F* hest = new TH2F("hest","Estimated energy",fEbins,logemin,logemax,fNTbins,fTbin);
246 TH2F* hmc = new TH2F("hmc","MC energy",fEbins,logemin,logemax,fNTbins,fTbin);
247
248 // borra
249 // TH1F* hest1 = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax);
250 // TH1F* hmc1 = new TH1F("hmc1","MC energy",fEbins,logemin,logemax);
251 // TH1F* coef1 = new TH1F("coef1","coefficients",fEbins,logemin,logemax);
252
253 // reset the coefficients
254 if(fCoeff)
255 delete fCoeff;
256 fCoeff = new TH2F("fcoeff","Coefficients for unfolding",fEbins,logemin,logemax,fNTbins,fTbin);
257
258 // fill the histograms of estimated and measured energy for weighted events
259 for(Int_t i=0;i<nentries;i++)
260 {
261 fCcut->GetEntry(i);
262
263 // mc and estimated energy
264 // OJO!! Estimated energy will be taken directly from some input container
265 Float_t emc = fMcEvt->GetEnergy();
266 Float_t estim = fHillas->GetSize()/0.18/15.;
267
268 if((emc<fEmin && estim<fEmin) || (emc>fEmax && estim>fEmax) ||
269 (emc<fEmin && estim>fEmax) || (emc>fEmax && estim<fEmin))
270 continue;
271
272 Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees)
273
274 // compute the bin where the weight is taken from
275 Int_t wbin = Int_t((TMath::Log10(emc)-fLogEWmin)/desub);
276
277 // fill the histograms
278 if(wbin<fEWbins)
279 {
280 if(fWeight[wbin]>0)
281 {
282 hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]);
283 hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
284 // borra
285// if(theta<fTbin[1])
286// {
287// hest1->Fill(TMath::Log10(estim),fWeight[wbin]);
288// hmc1->Fill(TMath::Log10(emc),fWeight[wbin]);
289// }
290 }
291 else
292 cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has no computed weight (lack of MC statistics), skipping" << endl;
293 }
294 else
295 cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
296 }
297
298 // divide the previous histos to get the coefficients for unfolding
299 for(Int_t j=0;j<fNTbins;j++)
300 for(Int_t i=0;i<fEbins;i++)
301 {
302 const Float_t num = hmc->GetBinContent(i+1,j+1);
303 const Float_t den = hest->GetBinContent(i+1,j+1);
304 //borra
305// const Float_t num1 = hmc1->GetBinContent(i+1);
306// const Float_t den1 = hest1->GetBinContent(i+1);
307 if(den)
308 {
309 fCoeff->SetBinContent(i+1,j+1,num/den);
310 //borra
311// if(j==0 && den1)
312// coef1->SetBinContent(i+1,num1/den1);
313 }
314 else
315 {
316 cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: energy bin " << i << ", thetabin " << j << " has no survivors, setting coefficient to 0" << endl;
317 fCoeff->SetBinContent(i+1,j+1,0.);
318 }
319 }
320
321
322 //borra
323// hmc1->Write();
324// hest1->Write();
325// coef1->Write();
326
327 delete hmc;
328 delete hest;
329}
330
331// --------------------------------------------------------------
332//
333// compute the coefficients used for the (iterative) unfolding
334//
335void MEffAreaAndCoeffCalc::ComputeEffectiveAreas()
336{
337 if(!fWeight)
338 {
339 cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: No weights computed! nothing done" << endl;
340 return;
341 }
342
343 //OJO!! radius should be read from somewhere!
344 const Float_t radius = 30000.; // generation radius (in cm)
345 const Float_t Atot = 3.14159*radius*radius; //total area (in cm^2)
346 const Float_t logemin = TMath::Log10(fEmin);
347 const Float_t logemax = TMath::Log10(fEmax);
348 const Int_t esbins = fEbins*fEsubbins; // total number of subbins
349 const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
350 const Float_t desub = de/fEsubbins; // subbin size (in log)
351
352 // reset the effective areas
353 if(fEffA)
354 delete fEffA;
355 fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin);
356 //borra
357// TH1F eff("eff","Effective area",fEbins,logemin,logemax);
358
359 // fill the spectrum of the survivors
360 TH2F* hpass= new TH2F("hpass","Survivors",esbins,logemin,logemax,fNTbins,fTbin);
361 TH2F* horig= new TH2F("horig","Original events",esbins,logemin,logemax,fNTbins,fTbin);
362
363 fCcut->Draw("theta:logenergy>>hpass","","goff");
364 fCini->Draw("theta:logenergy>>horig","","goff");
365
366 // compute the effective areas
367 for(Int_t j=0;j<fNTbins;j++)
368 for(Int_t i=0;i<fEbins;i++)
369 {
370 Float_t effarea =0;
371 Float_t wtot = 0;
372 for(Int_t k=0;k<fEsubbins;k++)
373 {
374 Float_t numerator = hpass->GetBinContent(i*fEsubbins+k+1,j+1);
375 Float_t denominator = horig->GetBinContent(i*fEsubbins+k+1,j+1);
376
377 if(denominator<=0)
378 cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy subbin " << i*fEsubbins+k <<", theta bin " << j << " contains no events" << endl;
379 else
380 {
381 const Float_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
382 const Float_t ww = fSpec->Eval(ew);
383 effarea+=Atot*numerator/denominator*ww;
384 wtot += ww;
385 }
386 }
387 if(!wtot)
388 {
389 cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy bin " << i <<", theta bin " << j << " contains no events setting effective area to 0" << endl;
390 fEffA->SetBinContent(i+1,j+1,0);
391 }
392 else
393 {
394 fEffA->SetBinContent(i+1,j+1,effarea/wtot);
395 // borra
396// if(j==0)
397// {
398// // cout << "*****" << i << ": " << effarea/wtot << endl;
399// eff.SetBinContent(i+1,effarea/wtot);
400// }
401 }
402 }
403
404 // borra
405// eff.Write();
406
407 delete hpass;
408 delete horig;
409}
410
411// --------------------------------------------------------------
412//
413// Call the internal functions to compute all the factors
414//
415void MEffAreaAndCoeffCalc::ComputeAllFactors()
416{
417 ComputeWeights();
418 ComputeCoefficients();
419 ComputeEffectiveAreas();
420}
Note: See TracBrowser for help on using the repository browser.