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

Last change on this file since 1962 was 1952, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 15.2 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-2003
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// As a default MHHadronness searches for the container "MHadronness".
54// This can be overwritten by a different pointer specified in the
55// Fill function (No type check is done!) Use a different name in
56// MFillH.
57//
58////////////////////////////////////////////////////////////////////////////
59#include "MHHadronness.h"
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 fIntGhness->SetMaximum(1.1);
109 fIntPhness->SetMaximum(1.1);
110
111 fQfac = new TGraph;
112 fQfac->SetTitle(" Naive Quality factor ");
113
114 fGhness->SetDirectory(NULL);
115 fPhness->SetDirectory(NULL);
116 fIntGhness->SetDirectory(NULL);
117 fIntPhness->SetDirectory(NULL);
118}
119
120// --------------------------------------------------------------------------
121//
122// Delete the histograms.
123//
124MHHadronness::~MHHadronness()
125{
126 delete fGhness;
127 delete fIntGhness;
128 delete fPhness;
129 delete fIntPhness;
130 delete fQfac;
131 delete fGraph;
132}
133
134// --------------------------------------------------------------------------
135//
136// Setup Filling of the histograms. It needs:
137// MMcEvt and MHadronness
138//
139Bool_t MHHadronness::SetupFill(const MParList *plist)
140{
141 fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
142 if (!fMcEvt)
143 {
144 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
145 return kFALSE;
146 }
147
148 fHadronness = (MHadronness*)plist->FindObject("MHadronness");
149
150 /*
151 MBinning* bins = (MBinning*)plist->FindObject("BinningHadronness");
152 if (!bins)
153 {
154 *fLog << err << dbginf << "BinningHadronness [MBinning] not found... aborting." << endl;
155 return kFALSE;
156 }
157
158 SetBinning(&fHist, binsalpha, binsenergy, binstheta);
159
160 fHist.Sumw2();
161 */
162
163 return kTRUE;
164}
165
166// --------------------------------------------------------------------------
167//
168// Fill the Hadronness from a MHadronness container into the corresponding
169// histogram dependant on the particle id.
170//
171// Every particle Id different than kGAMMA is considered a hadron.
172//
173// If you call Fill with a pointer (eg. from MFillH) this container is
174// used as a hadronness (Warning: No type check is done!) otherwise
175// the default is used ("MHadronness")
176//
177// Sometimes a distance is calculated as NaN (not a number). Such events
178// are skipped at the moment.
179//
180Bool_t MHHadronness::Fill(const MParContainer *par)
181{
182 // Preliminary Workaround: FIXME!
183 if (!par && !fHadronness)
184 {
185 *fLog << err << "MHHadronness::Fill: No MHadronness container specified!" << endl;
186 return kFALSE;
187 }
188
189 const MHadronness &had = par ? *(MHadronness*)par : *fHadronness;
190
191 const Double_t h = had.GetHadronness();
192
193 if (TMath::IsNaN(h))
194 return kCONTINUE;
195
196 if (fMcEvt->GetPartId()==kGAMMA)
197 fGhness->Fill(h);
198 else
199 fPhness->Fill(h);
200
201 return kTRUE;
202}
203
204Float_t MHHadronness::GetQ05() const
205{
206 Int_t n = fQfac->GetN();
207
208 Double_t val1x=0;
209 Double_t val2x=1;
210
211 Double_t val1y=0;
212 Double_t val2y=0;
213
214 for (Int_t i=1; i<=n; i++)
215 {
216 Double_t x, y;
217
218 fQfac->GetPoint(i, x, y);
219
220 if (x<0.5 && x>val1x)
221 {
222 val1x = x;
223 val1y = y;
224 }
225
226 if (x>0.5 && x<val2x)
227 {
228 val2x = x;
229 val2y = y;
230 }
231 }
232
233 return val2x-val1x == 0 ? 0 : val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
234}
235
236// --------------------------------------------------------------------------
237//
238// Finalize the histograms:
239// - integrate the hadroness histograms --> acceptance
240// - fill the Minimum Distance histogram (formular see class description)
241// - fill the Quality histogram (formular see class description)
242//
243Bool_t MHHadronness::Finalize()
244{
245 Int_t n = fGhness->GetNbinsX();
246
247 fGraph->Set(n);
248 fQfac->Set(n);
249
250 const Stat_t sumg = fGhness->Integral();
251 const Stat_t sump = fPhness->Integral();
252
253 *fLog << inf << "Sum Hadronness: gammas=" << sumg << " hadrons=" << sump << endl;
254
255 if (sumg==0 )
256 *fLog << warn << "Cannot calculate hadronness for 'gammas'." << endl;
257 if (sump==0)
258 *fLog << warn << "Cannot calculate hadronness for 'hadrons'." << endl;
259
260 // Normalize photon distribution
261 if (sumg>0)
262 fGhness->Scale(1./sumg);
263
264 // Normalize hadron distribution
265 if (sump>0)
266 fPhness->Scale(1./sump);
267
268 // Calculate acceptances
269 Float_t max=0;
270
271 for (Int_t i=1; i<=n; i++)
272 {
273 const Stat_t ip = fPhness->Integral(1, i);
274 const Stat_t ig = fGhness->Integral(1, i);
275
276 fIntPhness->SetBinContent(i, ip);
277 fIntGhness->SetBinContent(i, ig);
278
279 fGraph->SetPoint(i, ip, ig);
280
281 if (ip<=0)
282 continue;
283
284 const Double_t val = ig/sqrt(ip);
285 fQfac->SetPoint(i, ig, val);
286
287 if (val>max)
288 max = val;
289 }
290
291 fQfac->SetMaximum(max*1.05);
292
293 return kTRUE;
294}
295
296// --------------------------------------------------------------------------
297//
298// Search the corresponding points for the given hadron acceptance (acchad)
299// and interpolate the tow points (linear)
300//
301Double_t MHHadronness::GetGammaAcceptance(Double_t acchad) const
302{
303 const Int_t n = fGraph->GetN();
304 const Double_t *x = fGraph->GetX();
305 const Double_t *y = fGraph->GetY();
306
307 Int_t i = 0;
308 while (i<n && x[i]<acchad)
309 i++;
310
311 if (i==0 || i==n)
312 return 0;
313
314 if (i==n-1)
315 i--;
316
317 const Double_t x1 = x[i-1];
318 const Double_t y1 = y[i-1];
319
320 const Double_t x2 = x[i];
321 const Double_t y2 = y[i];
322
323 return (y2-y1)/(x2-x1) * (acchad-x2) + y2;
324}
325
326// --------------------------------------------------------------------------
327//
328// Search the corresponding points for the given gamma acceptance (accgam)
329// and interpolate the tow points (linear)
330//
331Double_t MHHadronness::GetHadronAcceptance(Double_t accgam) const
332{
333 const Int_t n = fGraph->GetN();
334 const Double_t *x = fGraph->GetX();
335 const Double_t *y = fGraph->GetY();
336
337 Int_t i = 0;
338 while (i<n && y[i]<accgam)
339 i++;
340
341 if (i==0 || i==n)
342 return 0;
343
344 if (i==n-1)
345 i--;
346
347 const Double_t x1 = y[i-1];
348 const Double_t y1 = x[i-1];
349
350 const Double_t x2 = y[i];
351 const Double_t y2 = x[i];
352
353 return (y2-y1)/(x2-x1) * (accgam-x2) + y2;
354}
355
356// --------------------------------------------------------------------------
357//
358// Print the corresponding Gammas Acceptance for a hadron acceptance of
359// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
360// acceptance and the corresponding acceptances and hadroness value is
361// printed, together with the maximum Q-factor.
362//
363void MHHadronness::Print(Option_t *) const
364{
365 *fLog << all;
366 *fLog << "Hadronness histograms:" << endl;
367 *fLog << "---------------------" << endl;
368
369 if (fGraph->GetN()==0)
370 {
371 *fLog << " <No Entries>" << endl;
372 return;
373 }
374
375 *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
376 *fLog << "acc(hadron) acc(gamma) acc(g)/acc(h)" << endl <<endl;
377
378 *fLog << " 0.005 " << Form("%6.3f", GetGammaAcceptance(0.005)) << " " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << endl;
379 *fLog << " 0.02 " << Form("%6.3f", GetGammaAcceptance(0.02)) << " " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << endl;
380 *fLog << " 0.05 " << Form("%6.3f", GetGammaAcceptance(0.05)) << " " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << endl;
381 *fLog << " 0.1 " << Form("%6.3f", GetGammaAcceptance(0.1 )) << " " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << endl;
382 *fLog << " 0.2 " << Form("%6.3f", GetGammaAcceptance(0.2 )) << " " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << endl;
383 *fLog << " 0.3 " << Form("%6.3f", GetGammaAcceptance(0.3 )) << " " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << endl;
384 *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << " 0.1 " << endl;
385 *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << " 0.2 " << endl;
386 *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << " 0.3 " << endl;
387 *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << " 0.4 " << endl;
388 *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << " 0.5 " << endl;
389 *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << " 0.6 " << endl;
390 *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << " 0.7 " << endl;
391 *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << " 0.8 " << endl;
392 *fLog << endl;
393
394 *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
395 *fLog << " Acc Hadrons = " << Form("%5.1f", GetHadronAcceptance(0.5)*100) << "%" << endl;
396 *fLog << endl;
397}
398
399// --------------------------------------------------------------------------
400//
401// Draw clone of all histograms. (For the Meaning see class description)
402//
403TObject *MHHadronness::DrawClone(Option_t *opt) const
404{
405 if (fGraph->GetN()==0)
406 return NULL;
407
408 TCanvas &c = *MakeDefCanvas("Hadronness", fTitle);
409 c.Divide(2, 2);
410
411 gROOT->SetSelectedPad(NULL);
412
413 c.cd(1);
414 //gStyle->SetOptStat(10);
415 Getghness()->DrawCopy();
416 Getphness()->SetLineColor(kRed);
417 Getphness()->DrawCopy("same");
418
419 c.cd(2);
420 //gStyle->SetOptStat(0);
421 Getighness()->DrawCopy();
422 Getiphness()->SetLineColor(kRed);
423 Getiphness()->DrawCopy("same");
424
425 c.cd(3);
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
442 c.cd(4);
443 gPad->Modified();
444 gPad->Update();
445 TGraph &g = (TGraph&)*fGraph->DrawClone("AC");
446 g.SetBit(kCanDelete);
447 gPad->Modified();
448 gPad->Update();
449 if (g.GetHistogram())
450 {
451 g.GetXaxis()->SetRangeUser(0, 1);
452 g.GetXaxis()->SetTitle("Acceptance Hadrons");
453 g.GetYaxis()->SetTitle("Acceptance Gammas");
454 g.SetMarkerStyle(kFullDotSmall);
455 g.Draw("P");
456
457 gPad->Modified();
458 gPad->Update();
459 }
460 /*
461 const Int_t h = fMinDist->GetMaximumBin();
462 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
463 fIntGhness->GetBinContent(h), kStar);
464 m->SetMarkerColor(kBlue);
465 m->SetBit(kCanDelete);
466 m->Draw();
467 */
468
469 gStyle->SetOptStat(1111);
470
471 return &c;
472}
473
474// --------------------------------------------------------------------------
475//
476// Draw all histograms. (For the Meaning see class description)
477//
478void MHHadronness::Draw(Option_t *)
479{
480 if (fGraph->GetN()==0)
481 return;
482
483 if (!gPad)
484 MakeDefCanvas("Hadronness", fTitle);
485
486 gPad->Divide(2, 2);
487
488 gPad->cd(1);
489 //gStyle->SetOptStat(10);
490 Getghness()->Draw();
491 Getphness()->SetLineColor(kRed);
492 Getphness()->Draw("same");
493
494 gPad->cd(2);
495 //gStyle->SetOptStat(0);
496 Getighness()->Draw();
497 Getiphness()->SetLineColor(kRed);
498 Getiphness()->Draw("same");
499
500 gPad->cd(3);
501 fQfac->Draw("A*");
502 gPad->Modified();
503 gPad->Update();
504 if (fQfac->GetHistogram())
505 {
506 fQfac->GetXaxis()->SetRangeUser(0, 1);
507 fQfac->GetXaxis()->SetTitle("Acceptance Gammas");
508 fQfac->GetYaxis()->SetTitle("Quality");
509 fQfac->SetMarkerStyle(kFullDotSmall);
510 fQfac->Draw("P");
511
512 gPad->Modified();
513 gPad->Update();
514 }
515 gPad->cd(4);
516 gPad->Modified();
517 gPad->Update();
518 fGraph->Draw("AC");
519 gPad->Modified();
520 gPad->Update();
521 if (fGraph->GetHistogram())
522 {
523 fGraph->GetXaxis()->SetRangeUser(0, 1);
524 fGraph->GetXaxis()->SetTitle("Acceptance Hadrons");
525 fGraph->GetYaxis()->SetTitle("Acceptance Gammas");
526 fGraph->SetMarkerStyle(kFullDotSmall);
527 fGraph->Draw("P");
528 gPad->Modified();
529 gPad->Update();
530 }
531}
Note: See TracBrowser for help on using the repository browser.