source: trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc@ 8936

Last change on this file since 8936 was 8935, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 16.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): Wolfgang Wittek 4/2002 <mailto:wittek@mppmu.mpg.de>
19! Author(s): Abelardo Moralejo 5/2003 <mailto:moralejo@pd.infn.it>
20! Author(s): Thomas Bretz 1/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2008
23!
24!
25\* ======================================================================== */
26
27//////////////////////////////////////////////////////////////////////////////
28//
29// MHEnergyEst
30//
31// calculates the migration matrix E-est vs. E-true
32// for different bins in Theta
33//
34// Class Version 2:
35// - fHResolution
36// + fHResolutionEst
37// + fHResolutionMC
38//
39//////////////////////////////////////////////////////////////////////////////
40#include "MHEnergyEst.h"
41
42#include <TF1.h>
43#include <TLine.h>
44#include <TMath.h>
45#include <TCanvas.h>
46#include <TStyle.h>
47#include <TProfile.h>
48#include <TPaveStats.h>
49
50#include "MString.h"
51
52#include "MMcEvt.hxx"
53
54#include "MBinning.h"
55#include "MParList.h"
56#include "MParameters.h"
57#include "MHMatrix.h"
58
59#include "MLog.h"
60#include "MLogManip.h"
61
62ClassImp(MHEnergyEst);
63
64using namespace std;
65
66// --------------------------------------------------------------------------
67//
68// Default Constructor. It sets name and title of the histogram.
69//
70MHEnergyEst::MHEnergyEst(const char *name, const char *title)
71 : fMcEvt(0), fEnergy(0), fResult(0), fMatrix(0)
72{
73 //
74 // set the name and title of this object
75 //
76 fName = name ? name : "MHEnergyEst";
77 fTitle = title ? title : "Histogram for the result of the energy reconstruction";
78
79 //fNameEnergy = "MEnergyEst";
80 //fNameResult = "MinimizationValue";
81
82 fHEnergy.SetDirectory(NULL);
83 fHEnergy.SetName("EnergyEst");
84 fHEnergy.SetTitle("Histogram in E_{est}, E_{mc} and \\Theta");
85 fHEnergy.SetXTitle("E_{est} [GeV]");
86 fHEnergy.SetYTitle("E_{mc} [GeV]");
87 fHEnergy.SetZTitle("\\Theta [\\circ]");
88
89 fHResolutionEst.SetDirectory(NULL);
90 fHResolutionEst.SetName("ResEnergyEst");
91 fHResolutionEst.SetTitle("Histogram in \\Delta E/E_{est} vs E_{est}");
92 fHResolutionEst.SetXTitle("E_{est} [GeV]");
93 fHResolutionEst.SetYTitle("1-E_{mc}/E_{est}");
94
95 fHResolutionMC.SetDirectory(NULL);
96 fHResolutionMC.SetName("ResEnergyMC");
97 fHResolutionMC.SetTitle("Histogram in \\Delta E/E_{mc} vs E_{mc}");
98 fHResolutionMC.SetXTitle("E_{mc} [GeV]");
99 fHResolutionMC.SetYTitle("E_{est}/E_{mc}-1");
100
101 fHImpact.SetDirectory(NULL);
102 fHImpact.SetName("Impact");
103 fHImpact.SetTitle("\\Delta E/E vs Impact parameter");
104 fHImpact.SetXTitle("Impact parameter [m]");
105 fHImpact.SetYTitle("E_{est}/E_{mc}-1");
106
107 MBinning binsi, binse, binst, binsr;
108 binse.SetEdgesLog(21, 6.3, 100000);
109 binst.SetEdgesASin(51, -0.005, 0.505);
110 binsr.SetEdges(75, -1.75, 1.75);
111
112 // Use the binning in impact to do efficiency studies
113 binsi.SetEdges(1, 0, 1000);
114
115 SetBinning(&fHEnergy, &binse, &binse, &binst);
116 SetBinning(&fHImpact, &binsi, &binsr);
117
118 SetBinning(&fHResolutionEst, &binse, &binsr);
119 SetBinning(&fHResolutionMC, &binse, &binsr);
120
121 // For some unknown reasons this must be called after
122 // the binning has been initialized at least once
123 fHEnergy.Sumw2();
124 fHImpact.Sumw2();
125 fHResolutionEst.Sumw2();
126 fHResolutionMC.Sumw2();
127}
128
129// --------------------------------------------------------------------------
130//
131// Set the binnings and prepare the filling of the histograms
132//
133Bool_t MHEnergyEst::SetupFill(const MParList *plist)
134{
135 if (!fMatrix)
136 {
137 fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
138 if (!fMcEvt)
139 {
140 *fLog << err << "MMcEvt not found... aborting." << endl;
141 return kFALSE;
142 }
143 }
144
145 fEnergy = (MParameterD*)plist->FindObject("MEnergyEst", "MParameterD");
146 if (!fEnergy)
147 {
148 *fLog << err << "MEnergyEst [MParameterD] not found... aborting." << endl;
149 return kFALSE;
150 }
151
152 fResult = (MParameterD*)const_cast<MParList*>(plist)->FindCreateObj("MParameterD", "MinimizationValue");
153 if (!fResult)
154 return kFALSE;
155
156 MBinning binst, binse, binsi, binsr;
157 binse.SetEdges(fHEnergy, 'x');
158 binst.SetEdges(fHEnergy, 'z');
159 binsi.SetEdges(fHImpact, 'x');
160 binsr.SetEdges(fHImpact, 'y');
161
162 binst.SetEdges(*plist, "BinningTheta");
163 binse.SetEdges(*plist, "BinningEnergyEst");
164 binsi.SetEdges(*plist, "BinningImpact");
165 binsr.SetEdges(*plist, "BinningEnergyRes");
166
167 SetBinning(&fHEnergy, &binse, &binse, &binst);
168 SetBinning(&fHImpact, &binsi, &binsr);
169
170 SetBinning(&fHResolutionEst, &binse, &binsr);
171 SetBinning(&fHResolutionMC, &binse, &binsr);
172
173 fChisq = 0;
174 fBias = 0;
175
176 fHEnergy.Reset();
177 fHImpact.Reset();
178
179 fHResolutionEst.Reset();
180 fHResolutionMC.Reset();
181
182 return kTRUE;
183}
184
185// --------------------------------------------------------------------------
186//
187// Fill the histogram
188//
189Bool_t MHEnergyEst::Fill(const MParContainer *par, const Stat_t w)
190{
191 const Double_t eest = fEnergy->GetVal();
192 const Double_t etru = fMatrix ? GetVal(0) : fMcEvt->GetEnergy();
193 const Double_t imp = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100;
194 const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
195
196 const Double_t resEst = (eest-etru)/eest;
197 const Double_t resMC = (eest-etru)/etru;
198
199 fHEnergy.Fill(eest, etru, theta, w);
200 fHImpact.Fill(imp, resEst, w);
201
202 fHResolutionEst.Fill(eest, resEst, w);
203 fHResolutionMC.Fill(etru, resMC, w);
204
205 // For the fit we use a different quantity
206 //const Double_t res = TMath::Log10(eest/etru);
207 const Double_t res = eest-etru;
208
209 fChisq += res*res;
210 fBias += res;
211
212 return kTRUE;
213}
214
215// --------------------------------------------------------------------------
216//
217// Divide chisq and bias by number of executions
218// Print result
219//
220Bool_t MHEnergyEst::Finalize()
221{
222 fChisq /= GetNumExecutions();
223 fBias /= GetNumExecutions();
224
225 fResult->SetVal(fChisq);
226
227 *fLog << all << endl;
228 Print();
229
230 return kTRUE;
231}
232
233// --------------------------------------------------------------------------
234//
235// Print result
236//
237void MHEnergyEst::Print(Option_t *o) const
238{
239 const Double_t res = TMath::Sqrt(fChisq-fBias*fBias);
240 if (!TMath::Finite(res))
241 {
242 *fLog << all << "MHEnergyEst::Print: ERROR - Resolution is not finite (eg. NaN)." << endl;
243 return;
244 }
245
246 TH1D *h = (TH1D*)fHResolutionEst.ProjectionY("Dummy", -1, -1, "s");
247 h->Fit("gaus", "Q0", "", -1.0, 0.25);
248
249 TF1 *f = h->GetFunction("gaus");
250
251 *fLog << all << "F=" << fChisq << endl;
252 *fLog << "Results from Histogram:" << endl;
253 *fLog << " Mean of Delta E/E: " << Form("%+4.2f%%", 100*h->GetMean()) << endl;
254 *fLog << " RMS of Delta E/E: " << Form("%4.2f%%", 100*h->GetRMS()) << endl;
255 *fLog << "Results from Histogram-Fit:" << endl;
256 *fLog << " Bias of Delta E/E: " << Form("%+4.2f%%", 100*f->GetParameter(1)) << endl;
257 *fLog << " Sigma of Delta E/E: " << Form("%4.2f%%", 100*f->GetParameter(2)) << endl;
258 *fLog << " Res of Delta E/E: " << Form("%4.2f%%", 100*TMath::Hypot(f->GetParameter(1), f->GetParameter(2))) << endl;
259
260 delete h;
261}
262
263// --------------------------------------------------------------------------
264//
265// Return Correction Coefficients (weights)
266// hist = E_mc/E_est
267//
268void MHEnergyEst::GetWeights(TH1D &hist) const
269{
270 // Project into EnergyEst_ey
271 // the "e" ensures that errors are calculated
272 TH1D *h1 = (TH1D*)fHEnergy.Project3D("rtlprmft_ex"); // E_Est
273 TH1D *h2 = (TH1D*)fHEnergy.Project3D("rtlprmft_ey"); // E_mc
274
275 h2->Copy(hist);
276 hist.Divide(h1);
277
278 delete h1;
279 delete h2;
280
281 hist.SetNameTitle("EnergyRatio", "Ratio between estimated and monte carlo energy");
282 hist.SetXTitle("E [GeV]");
283 hist.SetYTitle("N_{mc}/N_{est} [1]");
284 hist.SetDirectory(0);
285}
286
287void MHEnergyEst::Paint(Option_t *opt)
288{
289 TVirtualPad *pad = gPad;
290
291 pad->cd(1);
292
293 TH1 *hx=0;
294 TH1 *hy=0;
295
296 if (pad->GetPad(1))
297 {
298 pad->GetPad(1)->cd(1);
299
300 if (gPad->FindObject("EnergyEst_ex"))
301 hx = fHEnergy.Project3D("ex");
302
303 if (gPad->FindObject("EnergyEst_ey"))
304 hy = fHEnergy.Project3D("ey");
305
306 if (hx && hy)
307 {
308 hx->SetLineColor(kBlue);
309 hx->SetMarkerColor(kBlue);
310 hy->SetMaximum();
311 hy->SetMaximum(TMath::Max(hx->GetMaximum(), hy->GetMaximum())*1.2);
312 if (hy->GetMaximum()>0)
313 gPad->SetLogy();
314 }
315
316 if (pad->GetPad(1)->GetPad(2))
317 {
318 pad->GetPad(1)->GetPad(2)->cd(1);
319 if ((hx=(TH1D*)gPad->FindObject("EnergyEst_ez")))
320 fHEnergy.Project3D("ez");
321
322 pad->GetPad(1)->GetPad(2)->cd(2);
323 hx = (TH1D*)fHResolutionEst.ProjectionY("Resolution", -1, -1, "e");
324 TPaveStats *stats = dynamic_cast<TPaveStats*>(hx->FindObject("stats"));
325 if (stats)
326 {
327 stats->SetBit(BIT(17)); // TStyle.cxx: kTakeStyle=BIT(17)
328 stats->SetX1NDC(0.63);
329 stats->SetY1NDC(0.68);
330 }
331
332 hx->Fit("gaus", "Q", "", -1.0, 0.25);
333 hx->GetFunction("gaus")->SetLineColor(kBlue);
334 hx->GetFunction("gaus")->SetLineWidth(2);
335 gPad=NULL;
336 gStyle->SetOptFit(101);
337 }
338 }
339
340 if (pad->GetPad(2))
341 {
342 pad->GetPad(2)->cd(1);
343 if (gPad->FindObject("EnergyEst_yx"))
344 {
345 TH2D *hyx = static_cast<TH2D*>(fHEnergy.Project3D("yx"));
346 UpdateProf(*hyx, kTRUE);
347 }
348
349 TLine *l = (TLine*)gPad->FindObject("TLine");
350 if (l)
351 {
352 const Float_t min = TMath::Max(fHEnergy.GetXaxis()->GetXmin(), fHEnergy.GetYaxis()->GetXmin());
353 const Float_t max = TMath::Min(fHEnergy.GetXaxis()->GetXmax(), fHEnergy.GetYaxis()->GetXmax());
354
355 l->SetX1(min);
356 l->SetX2(max);
357 l->SetY1(min);
358 l->SetY2(max);
359 }
360
361 pad->GetPad(2)->cd(2);
362 UpdateProf(fHResolutionEst, kFALSE);
363
364 pad->GetPad(2)->cd(3);
365 UpdateProf(fHResolutionMC, kFALSE);
366 }
367}
368
369void MHEnergyEst::UpdateProf(TH2 &h, Bool_t logy)
370{
371 const TString pname = MString::Format("Prof%s", h.GetName());
372
373 if (!gPad->FindObject(pname))
374 return;
375
376 TH1D *hx = h.ProfileX(pname, -1, -1, "s");
377 hx->SetLineColor(kBlue);
378 hx->SetMarkerColor(kBlue);
379
380 if (logy && hx->GetMaximum()>0)
381 gPad->SetLogy();
382}
383
384TH1 *MHEnergyEst::MakeProj(const char *how)
385{
386 TH1 *p = fHEnergy.Project3D(how);
387 p->SetDirectory(NULL);
388 p->SetBit(kCanDelete);
389 p->SetBit(TH1::kNoStats);
390 p->SetMarkerStyle(kFullDotMedium);
391 p->SetLineColor(kBlue);
392
393 return p;
394}
395
396TH1 *MHEnergyEst::MakeProf(TH2 &h)
397{
398 TH1 *p = h.ProfileX(Form("Prof%s", h.GetName()), -1, -1, "s");
399 p->SetDirectory(NULL);
400 p->SetBit(kCanDelete);
401 p->SetLineWidth(2);
402 p->SetLineColor(kBlue);
403 p->SetFillStyle(4000);
404 p->SetStats(kFALSE);
405
406 return p;
407}
408
409// --------------------------------------------------------------------------
410//
411// Draw the histogram
412//
413void MHEnergyEst::Draw(Option_t *opt)
414{
415 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
416
417 // Do the projection before painting the histograms into
418 // the individual pads
419 AppendPad("");
420
421 pad->SetBorderMode(0);
422 pad->Divide(2, 1, 1e-10, 1e-10);
423
424 TH1 *h;
425
426 // ----------------------------------------
427
428 pad->cd(1);
429 gPad->SetBorderMode(0);
430
431 gPad->Divide(1, 2, 1e-10, 1e-10);
432
433 TVirtualPad *pad2 = gPad;
434
435 // ----------------------------------------
436
437 pad2->cd(1);
438 gPad->SetBorderMode(0);
439
440 gPad->SetGridx();
441 gPad->SetGridy();
442 gPad->SetLogx();
443
444 h = MakeProj("ey");
445 h->SetTitle("Energy disribution: Monte Carlo E_{mc} (black), Estimated E_{est} (blue)");
446 h->SetXTitle("E [GeV]"); // E_mc
447 h->SetYTitle("Counts");
448 h->Draw();
449
450 h = MakeProj("ex");
451 h->SetLineColor(kBlue);
452 h->SetMarkerColor(kBlue);
453 h->Draw("same");
454
455 // ----------------------------------------
456
457 pad2->cd(2);
458 gPad->SetBorderMode(0);
459
460 TVirtualPad *pad3 = gPad;
461 pad3->Divide(2, 1, 1e-10, 1e-10);
462 pad3->cd(1);
463 gPad->SetBorderMode(0);
464 gPad->SetGridx();
465 gPad->SetGridy();
466
467 h = MakeProj("ez");
468 h->SetTitle("Zenith Angle Distribution");
469 h->GetXaxis()->SetMoreLogLabels();
470 h->GetXaxis()->SetNoExponent();
471 h->Draw();
472
473 // ----------------------------------------
474
475 pad3->cd(2);
476 gPad->SetBorderMode(0);
477 gPad->SetGridx();
478 gPad->SetGridy();
479
480 h = fHResolutionEst.ProjectionY("_py");
481 h->SetTitle("Distribution of \\Delta E/E_{est}");
482 h->SetDirectory(NULL);
483 h->SetBit(kCanDelete);
484 h->GetXaxis()->SetRangeUser(-1.3, 1.3);
485 h->Draw();
486 // ----------------------------------------
487
488 pad->cd(2);
489 gPad->SetBorderMode(0);
490
491 gPad->Divide(1, 3, 1e-10, 1e-10);
492 pad2 = gPad;
493
494 // ----------------------------------------
495
496 pad2->cd(1);
497 gPad->SetBorderMode(0);
498 gPad->SetLogy();
499 gPad->SetLogx();
500 gPad->SetGridx();
501 gPad->SetGridy();
502
503 // Results in crashes....
504 //gROOT->GetListOfCleanups()->Add(gPad); // WHY?
505
506 TH2D *h2 = (TH2D*)fHEnergy.Project3D("yx");
507 h2->SetDirectory(NULL);
508 h2->SetBit(kCanDelete);
509 h2->SetFillColor(kBlue);
510 h2->SetTitle("Estimated energy E_{est} vs Monte Carlo energy E_{mc}");
511
512 TH1 *h1 = MakeProf(*h2);
513
514 h2->Draw("");
515 h1->Draw("E0 same");
516
517 TLine line;
518 line.DrawLine(0,0,1,1);
519
520 line.SetLineColor(kBlue);
521 line.SetLineWidth(2);
522 line.SetLineStyle(kDashed);
523
524 // ----------------------------------------
525
526 pad2->cd(2);
527 gPad->SetBorderMode(0);
528 gPad->SetLogx();
529 gPad->SetGridx();
530 gPad->SetGridy();
531 fHResolutionEst.Draw();
532 MakeProf(fHResolutionEst)->Draw("E0 same");
533
534 fHResolutionEst.GetXaxis()->SetMoreLogLabels();
535 fHResolutionEst.GetXaxis()->SetNoExponent();
536
537 line.DrawLine(fHResolutionEst.GetXaxis()->GetXmin(), 0, fHResolutionEst.GetXaxis()->GetXmax(), 0);
538
539 // ----------------------------------------
540
541 pad2->cd(3);
542 gPad->SetBorderMode(0);
543 gPad->SetLogx();
544 gPad->SetGridx();
545 gPad->SetGridy();
546 fHResolutionMC.Draw();
547 MakeProf(fHResolutionMC)->Draw("E0 same");
548
549 fHResolutionMC.GetXaxis()->SetMoreLogLabels();
550 fHResolutionMC.GetXaxis()->SetNoExponent();
551
552 line.DrawLine(fHResolutionMC.GetXaxis()->GetXmin(), 0, fHResolutionMC.GetXaxis()->GetXmax(), 0);
553}
554
555// --------------------------------------------------------------------------
556//
557// Returns the mapped value from the Matrix
558//
559Double_t MHEnergyEst::GetVal(Int_t i) const
560{
561 return (*fMatrix)[fMap[i]];
562}
563
564// --------------------------------------------------------------------------
565//
566// You can use this function if you want to use a MHMatrix instead of the
567// given containers. This function adds all necessary columns to the
568// given matrix. Afterward you should fill the matrix with the corresponding
569// data (eg from a file by using MHMatrix::Fill). If you now loop
570// through the matrix (eg using MMatrixLoop) MEnergyEstParam2::Process
571// will take the values from the matrix instead of the containers.
572//
573void MHEnergyEst::InitMapping(MHMatrix *mat)
574{
575 if (fMatrix)
576 return;
577
578 fMatrix = mat;
579
580 fMap[0] = fMatrix->AddColumn("MMcEvt.fEnergy");
581 fMap[1] = fMatrix->AddColumn("MMcEvt.fImpact/100");
582 fMap[2] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta*kRad2Deg");
583}
584
585void MHEnergyEst::StopMapping()
586{
587 fMatrix = NULL;
588}
Note: See TracBrowser for help on using the repository browser.