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

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