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

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