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

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