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

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