source: trunk/MagicSoft/Mars/mhist/MHHadronness.cc@ 1683

Last change on this file since 1683 was 1668, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 15.5 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): Thomas Bretz, 5/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Abelardo Moralejo <mailto:moralejo@pd.infn.it>
20!
21! Copyright: MAGIC Software Development, 2000-2002
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MHHadronness
29//
30// This is histogram is a way to evaluate the quality of a gamma/hadron
31// seperation method. It is filled from a MHadronness container, which
32// stores a hadroness for the current event. The Value must be in the
33// range [0,1]. To fill the histograms correctly the information
34// whether it is a gamma or hadron (not a gamma) must be available from
35// a MMcEvt container.
36//
37// In the constructor you can change the number of used bns for the
38// evaluation.
39//
40// The meaning of the histograms (Draw, DrawClone) are the following:
41// * Upper Left Corner:
42// - black: histogram of all hadronesses for gammas
43// - red: histogram of all hadronesses for non gammas
44// * Upper Right Corner:
45// - black: acceptance of gammas (Ag) vs. the hadroness
46// - red: acceptance of non gammas (Ah) vs. the hadroness
47// * Bottom Left Corner:
48// Naive quality factor: Ag/sqrt(Ah)
49// * Bottom Right Corner:
50// - black: Acceprtance Gammas vs. Acceptance Hadrons
51// - blue cross: minimum of distance to (0, 1)
52//
53////////////////////////////////////////////////////////////////////////////
54#include "MHHadronness.h"
55
56/*
57 // - blue: 2D distance of (acceptance_hadrons, acceptances_gammas)
58 // to optimum (0, 1): 1-sqrt(Ag*Ag + (1-Ah)*(1-Ah))
59 */
60
61#include <TPad.h>
62#include <TGraph.h>
63#include <TStyle.h>
64#include <TCanvas.h>
65#include <TMarker.h>
66
67#include "MLog.h"
68#include "MLogManip.h"
69
70#include "MParList.h"
71#include "MBinning.h"
72#include "MHadronness.h"
73
74#include "MMcEvt.hxx"
75
76ClassImp(MHHadronness);
77
78// --------------------------------------------------------------------------
79//
80// Setup histograms, nbins is the number of bins used for the evaluation.
81// The default is 100 bins.
82//
83MHHadronness::MHHadronness(Int_t nbins, const char *name, const char *title)
84{
85 //
86 // set the name and title of this object
87 //
88 fName = name ? name : "MHHadronness";
89 fTitle = title ? title : "Gamma/Hadron Separation Quality Histograms";
90
91 fGraph = new TGraph;
92 fGraph->SetTitle("Acceptance Gammas vs. Hadrons");
93 fGraph->SetMaximum(1);
94
95 fGhness = new TH1D("Ghness", "Hadronness", nbins, 0, 1);
96 fPhness = new TH1D("Phness", "Hadronness", nbins, 0, 1);
97 fGhness->SetXTitle("Hadronness");
98 fPhness->SetXTitle("Hadronness");
99 fGhness->SetYTitle("Counts");
100 fPhness->SetYTitle("Counts");
101
102 fIntGhness = new TH1D("AccGammas", "Acceptance", nbins, 0, 1);
103 fIntPhness = new TH1D("AccHadrons", "Acceptance", nbins, 0, 1);
104 fIntGhness->SetXTitle("Hadronness");
105 fIntPhness->SetXTitle("Hadronness");
106 fIntGhness->SetYTitle("Acceptance");
107 fIntPhness->SetYTitle("Acceptance");
108
109 /*
110 fQfac = new TH1D("Qfac", "Naive Quality factor", nbins, 0, 1);
111 fQfac->SetXTitle("Hadronness");
112 fQfac->SetYTitle("Quality");
113 */
114 fQfac = new TGraph;
115 fQfac->SetTitle(" Naive Quality factor ");
116
117 /*
118 fMinDist = new TH1D("MinDist", "Minimum Dist to (0, 1)", nbins, 0, 1);
119 fMinDist->SetXTitle("Hadronness");
120 fMinDist->SetYTitle("Distance");
121 */
122
123 //fQfac->SetDirectory(NULL);
124 fGhness->SetDirectory(NULL);
125 fPhness->SetDirectory(NULL);
126 //fMinDist->SetDirectory(NULL);
127 fIntGhness->SetDirectory(NULL);
128 fIntPhness->SetDirectory(NULL);
129}
130
131// --------------------------------------------------------------------------
132//
133// Delete the histograms.
134//
135MHHadronness::~MHHadronness()
136{
137 delete fGhness;
138 delete fIntGhness;
139 delete fPhness;
140 delete fIntPhness;
141 delete fQfac;
142 // delete fMinDist;
143 delete fGraph;
144}
145
146// --------------------------------------------------------------------------
147//
148// Setup Filling of the histograms. It needs:
149// MMcEvt and MHadronness
150//
151Bool_t MHHadronness::SetupFill(const MParList *plist)
152{
153 fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
154 if (!fMcEvt)
155 {
156 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
157 return kFALSE;
158 }
159
160 fHadronness = (MHadronness*)plist->FindObject("MHadronness");
161 if (!fHadronness)
162 {
163 *fLog << err << dbginf << "MHadronness not found... aborting." << endl;
164 return kFALSE;
165 }
166
167 /*
168 MBinning* bins = (MBinning*)plist->FindObject("BinningHadronness");
169 if (!bins)
170 {
171 *fLog << err << dbginf << "BinningHadronness [MBinning] not found... aborting." << endl;
172 return kFALSE;
173 }
174
175 SetBinning(&fHist, binsalpha, binsenergy, binstheta);
176
177 fHist.Sumw2();
178 */
179
180 return kTRUE;
181}
182
183// --------------------------------------------------------------------------
184//
185// Fill the Hadronness from a MHadronness container into the corresponding
186// histogram dependant on the particle id.
187//
188// Sometimes a distance is calculated as NaN (not a number). Such events
189// are skipped at the moment.
190//
191Bool_t MHHadronness::Fill(const MParContainer *par)
192{
193 // Preliminary Workaround: FIXME!
194
195 const Double_t h = fHadronness->GetHadronness();
196
197 if (TMath::IsNaN(h))
198 return kCONTINUE;
199
200 if (fMcEvt->GetPartId()==kGAMMA)
201 fGhness->Fill(h);
202 else
203 fPhness->Fill(h);
204
205 return kTRUE;
206}
207
208Float_t MHHadronness::GetQ05() const
209{
210 Int_t n = fQfac->GetN();
211
212 Double_t val1x=0;
213 Double_t val2x=1;
214
215 Double_t val1y=0;
216 Double_t val2y=0;
217
218 for (Int_t i=1; i<=n; i++)
219 {
220 Double_t x, y;
221
222 fQfac->GetPoint(i, x, y);
223
224 if (x<0.5 && x>val1x)
225 {
226 val1x = x;
227 val1y = y;
228 }
229
230 if (x>0.5 && x<val2x)
231 {
232 val2x = x;
233 val2y = y;
234 }
235 }
236
237 return val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
238}
239
240// --------------------------------------------------------------------------
241//
242// Finalize the histograms:
243// - integrate the hadroness histograms --> acceptance
244// - fill the Minimum Distance histogram (formular see class description)
245// - fill the Quality histogram (formular see class description)
246//
247Bool_t MHHadronness::Finalize()
248{
249 Int_t n = fGhness->GetNbinsX();
250
251 fGraph->Set(n);
252 fQfac->Set(n);
253
254 const Stat_t sumg = fGhness->Integral(1, n+1);
255 const Stat_t sump = fPhness->Integral(1, n+1);
256
257 Float_t max=0;
258
259 for (Int_t i=1; i<=n; i++)
260 {
261 const Stat_t ip = fPhness->Integral(1, i)/sump;
262 const Stat_t ig = fGhness->Integral(1, i)/sumg;
263
264 fIntPhness->SetBinContent(i, ip);
265 fIntGhness->SetBinContent(i, ig);
266
267 fGraph->SetPoint(i, ip, ig);
268
269 if (ip<=0)
270 continue;
271
272 //fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1)));
273
274 const Double_t val = ig/sqrt(ip);
275 fQfac->SetPoint(i, ig, val);
276
277 if (val>max)
278 max = val;
279 }
280
281 fQfac->SetMaximum(max*1.05);
282
283 return kTRUE;
284}
285
286// --------------------------------------------------------------------------
287//
288// Search the corresponding points for the given hadron acceptance (acchad)
289// and interpolate the tow points (linear)
290//
291Double_t MHHadronness::GetGammaAcceptance(Double_t acchad) const
292{
293 const Int_t n = fGraph->GetN();
294 const Double_t *x = fGraph->GetX();
295 const Double_t *y = fGraph->GetY();
296
297 Int_t i = 0;
298 while (i<n && x[i]<acchad)
299 i++;
300
301 if (i==0 || i==n)
302 return 0;
303
304 if (i==n-1)
305 i--;
306
307 const Double_t x1 = x[i-1];
308 const Double_t y1 = y[i-1];
309
310 const Double_t x2 = x[i];
311 const Double_t y2 = y[i];
312
313 return (y2-y1)/(x2-x1) * (acchad-x2) + y2;
314}
315
316// --------------------------------------------------------------------------
317//
318// Search the corresponding points for the given gamma acceptance (accgam)
319// and interpolate the tow points (linear)
320//
321Double_t MHHadronness::GetHadronAcceptance(Double_t accgam) const
322{
323 const Int_t n = fGraph->GetN();
324 const Double_t *x = fGraph->GetX();
325 const Double_t *y = fGraph->GetY();
326
327 Int_t i = 0;
328 while (i<n && y[i]<accgam)
329 i++;
330
331 if (i==0 || i==n)
332 return 0;
333
334 if (i==n-1)
335 i--;
336
337 const Double_t x1 = y[i-1];
338 const Double_t y1 = x[i-1];
339
340 const Double_t x2 = y[i];
341 const Double_t y2 = x[i];
342
343 return (y2-y1)/(x2-x1) * (accgam-x2) + y2;
344}
345
346// --------------------------------------------------------------------------
347//
348// Print the corresponding Gammas Acceptance for a hadron acceptance of
349// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
350// acceptance and the corresponding acceptances and hadroness value is
351// printed, together with the maximum Q-factor.
352//
353void MHHadronness::Print(Option_t *) const
354{
355 *fLog << all;
356 *fLog << "Hadronness histograms:" << endl;
357 *fLog << "---------------------" << endl;
358
359 if (fGraph->GetN()==0)
360 {
361 *fLog << " <No Entries>" << endl;
362 return;
363 }
364
365 *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
366 *fLog << "acc(hadron) acc(gamma) acc(g)/acc(h)" << endl <<endl;
367
368 *fLog << " 0.005 " << Form("%6.3f", GetGammaAcceptance(0.005)) << " " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << endl;
369 *fLog << " 0.02 " << Form("%6.3f", GetGammaAcceptance(0.02)) << " " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << endl;
370 *fLog << " 0.05 " << Form("%6.3f", GetGammaAcceptance(0.05)) << " " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << endl;
371 *fLog << " 0.1 " << Form("%6.3f", GetGammaAcceptance(0.1 )) << " " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << endl;
372 *fLog << " 0.2 " << Form("%6.3f", GetGammaAcceptance(0.2 )) << " " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << endl;
373 *fLog << " 0.3 " << Form("%6.3f", GetGammaAcceptance(0.3 )) << " " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << endl;
374 *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << " 0.1 " << endl;
375 *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << " 0.2 " << endl;
376 *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << " 0.3 " << endl;
377 *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << " 0.4 " << endl;
378 *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << " 0.5 " << endl;
379 *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << " 0.6 " << endl;
380 *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << " 0.7 " << endl;
381 *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << " 0.8 " << endl;
382 *fLog << endl;
383
384 /*
385 const Int_t h = fMinDist->GetMaximumBin();
386 *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " @ H=" << fMinDist->GetBinCenter(h) << endl;
387 *fLog << " Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, ";
388 *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl;
389 */
390
391 *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
392 *fLog << " Acc Hadrons = " << Form("%5.1f", GetHadronAcceptance(0.5)*100) << "%" << endl;
393 *fLog << endl;
394}
395
396// --------------------------------------------------------------------------
397//
398// Draw clone of all histograms. (For the Meaning see class description)
399//
400TObject *MHHadronness::DrawClone(Option_t *opt) const
401{
402 if (fGraph->GetN()==0)
403 return NULL;
404
405 TCanvas &c = *MakeDefCanvas("Hadronness", fTitle);
406 c.Divide(2, 2);
407
408 gROOT->SetSelectedPad(NULL);
409
410 c.cd(1);
411 gStyle->SetOptStat(10);
412 Getghness()->DrawCopy();
413 Getphness()->SetLineColor(kRed);
414 Getphness()->DrawCopy("same");
415 c.cd(2);
416 gStyle->SetOptStat(0);
417 Getighness()->DrawCopy();
418 Getiphness()->SetLineColor(kRed);
419 Getiphness()->DrawCopy("same");
420 /*
421 fMinDist->SetLineColor(kBlue);
422 fMinDist->DrawCopy("same");
423 */
424 c.cd(3);
425 //GetQfac()->DrawCopy();
426 TGraph &g2 = (TGraph&)*fQfac->DrawClone("A*");
427 g2.SetBit(kCanDelete);
428 gPad->Modified();
429 gPad->Update();
430 if (g2.GetHistogram())
431 {
432 g2.GetXaxis()->SetRangeUser(0, 1);
433 g2.GetXaxis()->SetTitle("Acceptance Gammas");
434 g2.GetYaxis()->SetTitle("Quality");
435 g2.SetMarkerStyle(kFullDotSmall);
436 g2.Draw("P");
437
438 gPad->Modified();
439 gPad->Update();
440 }
441 c.cd(4);
442 gPad->Modified();
443 gPad->Update();
444 TGraph &g = (TGraph&)*fGraph->DrawClone("AC");
445 g.SetBit(kCanDelete);
446 gPad->Modified();
447 gPad->Update();
448 if (g.GetHistogram())
449 {
450 g.GetXaxis()->SetRangeUser(0, 1);
451 g.GetXaxis()->SetTitle("Acceptance Hadrons");
452 g.GetYaxis()->SetTitle("Acceptance Gammas");
453 g.SetMarkerStyle(kFullDotSmall);
454 g.Draw("P");
455
456 gPad->Modified();
457 gPad->Update();
458 }
459 /*
460 const Int_t h = fMinDist->GetMaximumBin();
461 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
462 fIntGhness->GetBinContent(h), kStar);
463 m->SetMarkerColor(kBlue);
464 m->SetBit(kCanDelete);
465 m->Draw();
466 */
467 return &c;
468}
469
470// --------------------------------------------------------------------------
471//
472// Draw all histograms. (For the Meaning see class description)
473//
474void MHHadronness::Draw(Option_t *)
475{
476 if (fGraph->GetN()==0)
477 return;
478
479 if (!gPad)
480 MakeDefCanvas("Hadronness", fTitle);
481
482 gPad->Divide(2, 2);
483
484 gPad->cd(1);
485 gStyle->SetOptStat(10);
486 Getghness()->Draw();
487 Getphness()->SetLineColor(kRed);
488 Getphness()->Draw("same");
489 gPad->cd(2);
490 gStyle->SetOptStat(0);
491 Getighness()->Draw();
492 Getiphness()->SetLineColor(kRed);
493 Getiphness()->Draw("same");
494 /*
495 fMinDist->SetLineColor(kBlue);
496 fMinDist->Draw("same");
497 */
498 gPad->cd(3);
499 //GetQfac()->Draw();
500 fQfac->Draw("A*");
501 gPad->Modified();
502 gPad->Update();
503 if (fQfac->GetHistogram())
504 {
505 fQfac->GetXaxis()->SetRangeUser(0, 1);
506 fQfac->GetXaxis()->SetTitle("Acceptance Gammas");
507 fQfac->GetYaxis()->SetTitle("Quality");
508 fQfac->SetMarkerStyle(kFullDotSmall);
509 fQfac->Draw("P");
510
511 gPad->Modified();
512 gPad->Update();
513 }
514 gPad->cd(4);
515 gPad->Modified();
516 gPad->Update();
517 fGraph->Draw("AC");
518 gPad->Modified();
519 gPad->Update();
520 if (fGraph->GetHistogram())
521 {
522 fGraph->GetXaxis()->SetRangeUser(0, 1);
523 fGraph->GetXaxis()->SetTitle("Acceptance Hadrons");
524 fGraph->GetYaxis()->SetTitle("Acceptance Gammas");
525 fGraph->SetMarkerStyle(kFullDotSmall);
526 fGraph->Draw("P");
527 gPad->Modified();
528 gPad->Update();
529 }
530 /*
531 const Int_t h = fMinDist->GetMaximumBin();
532 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
533 fIntGhness->GetBinContent(h), kStar);
534 m->SetMarkerColor(kBlue);
535 m->SetBit(kCanDelete);
536 m->Draw();
537 */
538}
Note: See TracBrowser for help on using the repository browser.