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

Last change on this file since 1874 was 1809, checked in by wittek, 22 years ago
*** empty log message ***
File size: 16.6 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 Float_t retValue;
238 if (val2x-val1x != 0.0)
239 retValue = val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
240 else
241 retValue = 0.0;
242
243 return retValue;
244}
245
246// --------------------------------------------------------------------------
247//
248// Finalize the histograms:
249// - integrate the hadroness histograms --> acceptance
250// - fill the Minimum Distance histogram (formular see class description)
251// - fill the Quality histogram (formular see class description)
252//
253Bool_t MHHadronness::Finalize()
254{
255 Int_t n = fGhness->GetNbinsX();
256
257 fGraph->Set(n);
258 fQfac->Set(n);
259
260 Stat_t sumg;
261 Stat_t sump;
262
263 sumg = fGhness->Integral(1, n);
264 sump = fPhness->Integral(1, n);
265
266 if (sumg == 0.0 || sump == 0.0)
267 {
268 *fLog << "MHHadronness::Finalize; sumg or sump is zero; sumg, sump = "
269 << sumg << ", " << sump << ". Cannot calculate hadronness"
270 << endl;
271 }
272
273
274 // Normalize photon distribution
275 Stat_t con;
276 if (sumg > 0.0)
277 for (Int_t i=1; i<=n; i++)
278 {
279 con = (fGhness->GetBinContent(i)) / sumg;
280 fGhness->SetBinContent(i, con);
281 }
282
283 // Normalize hadron distribution
284 if (sump > 0.0)
285 for (Int_t i=1; i<=n; i++)
286 {
287 con = (fPhness->GetBinContent(i)) / sump;
288 fPhness->SetBinContent(i, con);
289 }
290
291 // Calculate acceptances
292 sumg = fGhness->Integral(1, n);
293 sump = fPhness->Integral(1, n);
294
295 *fLog << "MHHadronness::Finalize; sumg, sump = " << sumg << ", "
296 << sump << endl;
297
298 Float_t max=0;
299
300 for (Int_t i=1; i<=n; i++)
301 {
302 Stat_t ip;
303 if (sump != 0.0)
304 ip = fPhness->Integral(1, i)/sump;
305 else
306 ip = 0;
307
308 Stat_t ig;
309 if (sumg != 0.0)
310 ig = fGhness->Integral(1, i)/sumg;
311 else
312 ig = 0;
313
314 fIntPhness->SetBinContent(i, ip);
315 fIntGhness->SetBinContent(i, ig);
316
317 fGraph->SetPoint(i, ip, ig);
318
319 if (ip<=0)
320 continue;
321
322 //fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1)));
323
324 const Double_t val = ig/sqrt(ip);
325 fQfac->SetPoint(i, ig, val);
326
327 if (val>max)
328 max = val;
329 }
330
331 fQfac->SetMaximum(max*1.05);
332
333 return kTRUE;
334}
335
336// --------------------------------------------------------------------------
337//
338// Search the corresponding points for the given hadron acceptance (acchad)
339// and interpolate the tow points (linear)
340//
341Double_t MHHadronness::GetGammaAcceptance(Double_t acchad) const
342{
343 const Int_t n = fGraph->GetN();
344 const Double_t *x = fGraph->GetX();
345 const Double_t *y = fGraph->GetY();
346
347 Int_t i = 0;
348 while (i<n && x[i]<acchad)
349 i++;
350
351 if (i==0 || i==n)
352 return 0;
353
354 if (i==n-1)
355 i--;
356
357 const Double_t x1 = x[i-1];
358 const Double_t y1 = y[i-1];
359
360 const Double_t x2 = x[i];
361 const Double_t y2 = y[i];
362
363 return (y2-y1)/(x2-x1) * (acchad-x2) + y2;
364}
365
366// --------------------------------------------------------------------------
367//
368// Search the corresponding points for the given gamma acceptance (accgam)
369// and interpolate the tow points (linear)
370//
371Double_t MHHadronness::GetHadronAcceptance(Double_t accgam) const
372{
373 const Int_t n = fGraph->GetN();
374 const Double_t *x = fGraph->GetX();
375 const Double_t *y = fGraph->GetY();
376
377 Int_t i = 0;
378 while (i<n && y[i]<accgam)
379 i++;
380
381 if (i==0 || i==n)
382 return 0;
383
384 if (i==n-1)
385 i--;
386
387 const Double_t x1 = y[i-1];
388 const Double_t y1 = x[i-1];
389
390 const Double_t x2 = y[i];
391 const Double_t y2 = x[i];
392
393 return (y2-y1)/(x2-x1) * (accgam-x2) + y2;
394}
395
396// --------------------------------------------------------------------------
397//
398// Print the corresponding Gammas Acceptance for a hadron acceptance of
399// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
400// acceptance and the corresponding acceptances and hadroness value is
401// printed, together with the maximum Q-factor.
402//
403void MHHadronness::Print(Option_t *) const
404{
405 *fLog << all;
406 *fLog << "Hadronness histograms:" << endl;
407 *fLog << "---------------------" << endl;
408
409 if (fGraph->GetN()==0)
410 {
411 *fLog << " <No Entries>" << endl;
412 return;
413 }
414
415 *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
416 *fLog << "acc(hadron) acc(gamma) acc(g)/acc(h)" << endl <<endl;
417
418 *fLog << " 0.005 " << Form("%6.3f", GetGammaAcceptance(0.005)) << " " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << endl;
419 *fLog << " 0.02 " << Form("%6.3f", GetGammaAcceptance(0.02)) << " " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << endl;
420 *fLog << " 0.05 " << Form("%6.3f", GetGammaAcceptance(0.05)) << " " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << endl;
421 *fLog << " 0.1 " << Form("%6.3f", GetGammaAcceptance(0.1 )) << " " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << endl;
422 *fLog << " 0.2 " << Form("%6.3f", GetGammaAcceptance(0.2 )) << " " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << endl;
423 *fLog << " 0.3 " << Form("%6.3f", GetGammaAcceptance(0.3 )) << " " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << endl;
424 *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << " 0.1 " << endl;
425 *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << " 0.2 " << endl;
426 *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << " 0.3 " << endl;
427 *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << " 0.4 " << endl;
428 *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << " 0.5 " << endl;
429 *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << " 0.6 " << endl;
430 *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << " 0.7 " << endl;
431 *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << " 0.8 " << endl;
432 *fLog << endl;
433
434 /*
435 const Int_t h = fMinDist->GetMaximumBin();
436 *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " @ H=" << fMinDist->GetBinCenter(h) << endl;
437 *fLog << " Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, ";
438 *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl;
439 */
440
441 *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
442 *fLog << " Acc Hadrons = " << Form("%5.1f", GetHadronAcceptance(0.5)*100) << "%" << endl;
443 *fLog << endl;
444}
445
446// --------------------------------------------------------------------------
447//
448// Draw clone of all histograms. (For the Meaning see class description)
449//
450TObject *MHHadronness::DrawClone(Option_t *opt) const
451{
452 if (fGraph->GetN()==0)
453 return NULL;
454
455 TCanvas &c = *MakeDefCanvas("Hadronness", fTitle);
456 c.Divide(2, 2);
457
458 gROOT->SetSelectedPad(NULL);
459
460 c.cd(1);
461 //gStyle->SetOptStat(10);
462 Getghness()->DrawCopy();
463 Getphness()->SetLineColor(kRed);
464 Getphness()->DrawCopy("same");
465 c.cd(2);
466 //gStyle->SetOptStat(0);
467 Getighness()->DrawCopy();
468 Getiphness()->SetLineColor(kRed);
469 Getiphness()->DrawCopy("same");
470 /*
471 fMinDist->SetLineColor(kBlue);
472 fMinDist->DrawCopy("same");
473 */
474 c.cd(3);
475 //GetQfac()->DrawCopy();
476 TGraph &g2 = (TGraph&)*fQfac->DrawClone("A*");
477 g2.SetBit(kCanDelete);
478 gPad->Modified();
479 gPad->Update();
480 if (g2.GetHistogram())
481 {
482 g2.GetXaxis()->SetRangeUser(0, 1);
483 g2.GetXaxis()->SetTitle("Acceptance Gammas");
484 g2.GetYaxis()->SetTitle("Quality");
485 g2.SetMarkerStyle(kFullDotSmall);
486 g2.Draw("P");
487
488 gPad->Modified();
489 gPad->Update();
490 }
491 c.cd(4);
492 gPad->Modified();
493 gPad->Update();
494 TGraph &g = (TGraph&)*fGraph->DrawClone("AC");
495 g.SetBit(kCanDelete);
496 gPad->Modified();
497 gPad->Update();
498 if (g.GetHistogram())
499 {
500 g.GetXaxis()->SetRangeUser(0, 1);
501 g.GetXaxis()->SetTitle("Acceptance Hadrons");
502 g.GetYaxis()->SetTitle("Acceptance Gammas");
503 g.SetMarkerStyle(kFullDotSmall);
504 g.Draw("P");
505
506 gPad->Modified();
507 gPad->Update();
508 }
509 /*
510 const Int_t h = fMinDist->GetMaximumBin();
511 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
512 fIntGhness->GetBinContent(h), kStar);
513 m->SetMarkerColor(kBlue);
514 m->SetBit(kCanDelete);
515 m->Draw();
516 */
517 return &c;
518}
519
520// --------------------------------------------------------------------------
521//
522// Draw all histograms. (For the Meaning see class description)
523//
524void MHHadronness::Draw(Option_t *)
525{
526 if (fGraph->GetN()==0)
527 return;
528
529 if (!gPad)
530 MakeDefCanvas("Hadronness", fTitle);
531
532 gPad->Divide(2, 2);
533
534 gPad->cd(1);
535 //gStyle->SetOptStat(10);
536 Getghness()->Draw();
537 Getphness()->SetLineColor(kRed);
538 Getphness()->Draw("same");
539 gPad->cd(2);
540 //gStyle->SetOptStat(0);
541 Getighness()->Draw();
542 Getiphness()->SetLineColor(kRed);
543 Getiphness()->Draw("same");
544 /*
545 fMinDist->SetLineColor(kBlue);
546 fMinDist->Draw("same");
547 */
548 gPad->cd(3);
549 //GetQfac()->Draw();
550 fQfac->Draw("A*");
551 gPad->Modified();
552 gPad->Update();
553 if (fQfac->GetHistogram())
554 {
555 fQfac->GetXaxis()->SetRangeUser(0, 1);
556 fQfac->GetXaxis()->SetTitle("Acceptance Gammas");
557 fQfac->GetYaxis()->SetTitle("Quality");
558 fQfac->SetMarkerStyle(kFullDotSmall);
559 fQfac->Draw("P");
560
561 gPad->Modified();
562 gPad->Update();
563 }
564 gPad->cd(4);
565 gPad->Modified();
566 gPad->Update();
567 fGraph->Draw("AC");
568 gPad->Modified();
569 gPad->Update();
570 if (fGraph->GetHistogram())
571 {
572 fGraph->GetXaxis()->SetRangeUser(0, 1);
573 fGraph->GetXaxis()->SetTitle("Acceptance Hadrons");
574 fGraph->GetYaxis()->SetTitle("Acceptance Gammas");
575 fGraph->SetMarkerStyle(kFullDotSmall);
576 fGraph->Draw("P");
577 gPad->Modified();
578 gPad->Update();
579 }
580 /*
581 const Int_t h = fMinDist->GetMaximumBin();
582 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
583 fIntGhness->GetBinContent(h), kStar);
584 m->SetMarkerColor(kBlue);
585 m->SetBit(kCanDelete);
586 m->Draw();
587 */
588}
Note: See TracBrowser for help on using the repository browser.