source: trunk/MagicSoft/Mars/mtemp/mucm/classes/MHMyFindSignificanceONOFF.cc@ 6634

Last change on this file since 6634 was 6634, checked in by marcos, 20 years ago
*** empty log message ***
File size: 23.1 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): Marcols Lopez, August 2004 <mailto:marcos@gae.ucm.es>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MHMyFindSignificance
28//
29//
30// Non: # events in the signal region (0,fAlphasig) of the ON-Data hist.
31// Nbg: # events in the signal region (0,fAlphasig) of the OFF-Data hist.
32// Nex = Non - Nbg
33//
34/////////////////////////////////////////////////////////////////////////////
35#include "MHMyFindSignificanceONOFF.h"
36
37#include <TH1.h>
38#include <TF1.h>
39#include <TCanvas.h>
40#include <TPad.h>
41#include <TPaveText.h>
42#include <TGaxis.h>
43#include <TStyle.h>
44#include <TLine.h>
45#include <TString.h>
46#include <TGraph.h>
47#include <TPaveStats.h>
48
49#include "MLog.h"
50#include "MLogManip.h"
51
52
53ClassImp(MHMyFindSignificanceONOFF);
54
55using namespace std;
56
57
58// --------------------------------------------------------------------------
59//
60// Constructor
61//
62MHMyFindSignificanceONOFF::MHMyFindSignificanceONOFF(const char *name, const char *title)
63 : fNon(0.0), fdNon(0.0), fNbg(0.0), fdNbg(0.0), fNex(0.0), fdNex(0.0),
64 fSigLiMa(0.0),
65 fNbgFit(0.0), fdNbgFit(0.0), fNexFit(0.0), fdNexFit(0.0), fGamma(0.0),
66 fNoff(0.0), fSigLiMaFit(0.0),
67 fBestSigma(0.0), fBestAlphaCut(0.0), fBestExcess(0.0)
68{
69 fName = name ? name : "MHMyFindSignificanceONOFF";
70 fTitle = title ? title : "Find Significance in alpha plot";
71
72 fHistON = NULL;
73 fHistOrigON = NULL;
74 fHistOFF = NULL;
75 fHistOrigOFF = NULL;
76 fHistOFFNormalized = NULL;
77
78 // allow rebinning of the alpha plot
79 fRebin = kTRUE;
80
81 // allow reducing the degree of the polynomial
82 fReduceDegree = kTRUE;
83}
84
85
86
87
88// ---------------------------------------------------------------------------
89//
90// Calculate
91//
92// - fNon: number of events in the signal region from ON data
93//
94// and 3 different estimation of the background:
95//
96// - fNbgOFF: number of events in the signal region from OFF data
97// - fNbgFitOFF: estimation of # events in the signal region from a fit of
98// the OFF data in the Bkg region.
99// - fNbgFitON: estimation of # events in the signal region from a fit of
100// the ON data in the Bkg region.
101// - fNex
102//
103Bool_t MHMyFindSignificanceONOFF::FindSigmaONOFF(TH1 *histON, TH1 *histOFF,
104 Double_t alphasig, Double_t alphamin, Double_t alphamax, Int_t degree)
105{
106
107 fHistOrigON = histON;
108 fHistOrigOFF = histOFF;
109
110 fAlphasig = alphasig;
111 fAlphaminON = alphamin;
112 fAlphamaxON = alphamax;
113 fAlphaminOFF = 0; // <--------------------
114 fAlphamaxOFF = alphamax;
115 fDegree = degree;
116
117
118
119
120 // -----------------------------------------------------------------
121 //
122 // Calculate the normalization factor in range (alphanin,alphamax)
123 //
124 // FIXME: This can be done directly with:
125 // fFindsigON.GetNbtot()/fFindsigOFF.GetNbtot()
126 // but for that the polynomial for OFF-Data has to be fitted
127 // (alphamin, alphamax) instead of (0,alphamax)
128 //
129 CalcNormalizationFactor(fHistOrigON, fHistOrigOFF, fAlphaminON, fAlphamaxON);
130
131
132 // -----------------------------------------------------------------------
133 //
134 // Use MHFindSignificance to find the # evts in the signal and Bkg region
135 // and to fit the histograms.
136 // The ON hist is fitted to a polynomial in reg. (fAlphaminON,fAlphamaxON)
137 // and to a gaussian+polynomial.
138 // The OFF hist is only fitted to a polynomial in(fAlphaminOFF,fAlphamaxOFF
139 //
140 fFindsigON.SetReduceDegree(fReduceDegree);
141 fFindsigOFF.SetReduceDegree(fReduceDegree);
142 fFindsigON.SetRebin(fRebin);
143 fFindsigOFF.SetRebin(fRebin);
144
145 const Bool_t rc1 = fFindsigON.FindSigma(fHistOrigON, fAlphaminON, fAlphamaxON, fDegree, fAlphasig, 0,1,0);
146
147 const Bool_t rc2 = fFindsigOFF.FindSigma(fHistOrigOFF, fAlphaminOFF,fAlphamaxOFF, fDegree, fAlphasig, 0,0,0);
148
149 if (rc1==kFALSE || rc2==kFALSE)
150 {
151 *fLog << err << dbginf << "MHFindSignificance retuns kFALSE." << endl;
152 return kFALSE;
153 }
154
155
156 // Check consistency
157 if (fabs(fFindsigON.GetAlphasi() - fFindsigOFF.GetAlphasi()) > fEps)
158 {
159 *fLog << "fAlphasiOFF (" << fFindsigOFF.GetAlphasi()<< ") is not equal to fAlphasi ("
160 << fFindsigON.GetAlphasi() << "), this is something that should not happen"
161 << endl;
162 //return kFALSE; It might happen in pathological cases (very few OFF)
163 }
164
165
166 //
167 // Get results of the fits
168 //
169 fHistON = fFindsigON.GetHist(); // don't use the original hist to draw
170 fHistOFF = fFindsigOFF.GetHist();
171 // From ON-Data
172 fNon = fFindsigON.GetNon();
173 fdNon = fFindsigON.GetdNon();
174 // From OFF-Data
175 fNbg = fFindsigOFF.GetNon(); // counts in signal reg. OFF hist NO
176 fdNbg = fFindsigOFF.GetdNon(); //normalized !!!
177 // From OFF-Data polynomial fit
178 fNbgFit = fFindsigOFF.GetNbg();
179 fdNbgFit = fFindsigOFF.GetdNbg();
180 fGamma = fFindsigOFF.GetGamma() * fNormFactor; //<-----------
181 fNoff = fFindsigOFF.GetNoff();
182
183
184 //
185 // Compare the counter number of bg events, NbgMea, and the fitted, NbgFit
186 //
187 if (fabs(fNbg - fNbgFit) > 0.1 * fNbgFit)
188 {
189 *fLog << err << dbginf << "number of OFF events and Fitted number of OFF events in signal region do not agree (within 10 %)" << endl;
190
191 *fLog << "fNbg = " << fNbg << " ; fNbgFit = " << fNbgFit << endl;
192 }
193
194
195
196 // ---------------------------------------------------
197 //
198 // calculate the number of excess events in the signal region
199 //
200 fNex = fNon - fNbg * fNormFactor;
201 fNexFit = fNon - fNbgFit * fNormFactor;
202
203
204 //
205 // Calculate Significance
206 //
207 // Significance computed using counted number of OFF-data events in signal
208 // region (fNbg) and normalization factor (fNormFactor).
209 // This is the strictly speaking the significance in Li&Ma paper...
210 //
211 if ( !SigmaLiMa(fNon, fNbg, fNormFactor, &fSigLiMa) )
212 {
213 *fLog << err << dbginf << "SigmaLiMa failed" << endl;
214 return kFALSE;
215 }
216 //
217 // Significance computed using effective number of OFF events in signal
218 // region (fNoff) and gamma factor (fGama).
219 // This is Wolfgang approach to the calulation of significance
220 // using Li&Ma formula and estimated OFF events from polynomial fit.
221 //
222 if ( !SigmaLiMa(fNon, fNoff, fGamma , &fSigLiMaFit) )
223 {
224 *fLog << err << dbginf << "SigmaLiMaFit failed" << endl;
225 return kFALSE;
226 }
227
228
229 //
230 // calculate the error of the number of excess events
231 // using fitted quantities and counted quantities
232 //
233 fdNex = fNex / fSigLiMa;
234 fdNexFit = fNexFit / fSigLiMaFit;
235
236
237
238 //
239 // Create histogram of OFF data scaled by the normalization factor
240 //
241 NormalizedHistOFF();
242
243
244 return kTRUE;
245}
246
247
248
249// --------------------------------------------------------------------------
250//
251// Create a new histogram to store the hist of OFF data scale with by the
252// normalization factor
253//
254void MHMyFindSignificanceONOFF::NormalizedHistOFF()
255{
256
257 //
258 // We need the following factor in case that one of the histogram
259 // (ON or OFF) is rebin when calculatin the sifnigicance
260 //
261 const Double_t BinWidthRatioONOFF = fHistON->GetBinWidth(1)/fHistOFF->GetBinWidth(1) ;
262
263
264 //
265 // Clone OFF histogram and scale it
266 //
267 fHistOFFNormalized = (TH1*)fHistOFF->Clone();
268 fHistOFFNormalized->SetName(Form("%s (Normalized)",fHistOFF->GetName()));
269
270 fHistOFFNormalized->Scale(fNormFactor*BinWidthRatioONOFF);
271
272
273 //
274 // Scale the fit function of OFF data
275 //
276 TF1* fPolyOFFNormalized = (TF1*)fFindsigOFF.GetPoly()->Clone();
277
278 fPolyOFFNormalized->SetLineColor(2);
279 fPolyOFFNormalized->SetLineWidth(2);
280
281
282 const Int_t npar = fPolyOFFNormalized->GetNpar();
283 Double_t Parameter = 0.0;
284 Double_t ParameterError = 0.0;
285 for(Int_t i=0; i<npar; i++)
286 {
287 Parameter = fNormFactor*BinWidthRatioONOFF *fPolyOFFNormalized->GetParameter(i);
288 ParameterError = fNormFactor*BinWidthRatioONOFF *fPolyOFFNormalized->GetParError(i);
289
290 fPolyOFFNormalized -> SetParameter(i, Parameter);
291 fPolyOFFNormalized -> SetParError(i,ParameterError);
292 }
293
294 TList *funclistOFF = fHistOFFNormalized->GetListOfFunctions();
295 funclistOFF->Clear(); // delete the functions cloned from fHistOFF
296 funclistOFF->Add(fPolyOFFNormalized);
297
298
299}
300
301
302
303// --------------------------------------------------------------------------
304//
305// Draw hte Alpha plot for the ON data, and on top of it, the normalized Alpha
306// plot for the OFF data.
307//
308void MHMyFindSignificanceONOFF::Draw(Option_t* option)
309{
310
311 if(!fHistON || !fHistOFF)
312 return;
313
314
315
316 //
317 // If a canvas exists use it, otherwise Creates a new canvas
318 //
319 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
320 pad->SetBorderMode(0);
321
322 // pad->Divide(2,1);
323
324 AppendPad("");
325
326 pad->cd(1);
327
328 //
329 // Draw histograms and fit functions
330 //
331 fHistON->SetTitle("");
332 fHistON->SetMarkerStyle(21);
333 fHistON->Draw("p");
334
335 fHistOFFNormalized->SetLineColor(2);
336 fHistOFFNormalized->SetMarkerStyle(25);
337 fHistOFFNormalized->SetMarkerColor(2);
338 fHistOFFNormalized->Draw("Psame");
339
340 TList *funclistON = fHistON->GetListOfFunctions();
341 TF1* polyON = (TF1*)funclistON->FindObject("Poly");
342 TF1* gaussON = (TF1*)funclistON->FindObject("PolyGauss");
343 TF1* backgON = (TF1*)funclistON->FindObject("Backg");
344
345 if(gaussON)
346 {
347 gaussON->SetLineColor(4);
348 gaussON->SetLineWidth(4);
349 }
350 if(backgON)
351 funclistON->Remove(backgON);
352 if(polyON)
353 {
354 polyON->SetLineStyle(2);
355 polyON->SetLineColor(2);
356 polyON->SetLineWidth(2);
357 }
358
359
360 //
361 // Draw results onto the figure
362 //
363 DrawResults(option);
364
365
366 gPad->Modified();
367 gPad->Update();
368
369
370
371
372 // //
373// // Draw the difference of ON-OFF
374// //
375// pad->cd(2);
376// TH1D* h = (TH1D*)fHistON->Clone();
377// h->GetListOfFunctions()->Clear();
378// h->SetTitle("Alpha Plot ON-OFF");
379// h->Add(fHistOFFNormalized,-1.);
380// h->SetFillColor(17);
381// h->Draw("hist");
382
383// TLine* line = new TLine(0, 0, 90, 0);
384// line->Draw();
385
386// DrawResults("short");
387
388}
389
390
391
392// --------------------------------------------------------------------------
393//
394// Draw the results on top the histograms pad
395//
396void MHMyFindSignificanceONOFF::DrawResults(Option_t* option)
397{
398
399 TString opt(option);
400
401 TPaveText *pt;
402 char tx[100];
403
404
405 if(opt.Contains("short",TString::kIgnoreCase))
406 {
407 pt= new TPaveText(0.24, 0.59, 0.86, 0.87, "NDC");
408
409
410 // sprintf(tx, " Nbg measured (Normalized) = %8.1f #pm %8.1f",
411// fNbg * fNormFactor, fdNbg * fNormFactor);
412// pt->AddText(tx);
413
414// sprintf(tx, " Nex (ON - OFF measured) = %8.1f #pm %8.1f",fNex,fdNex);
415// pt->AddText(tx);
416
417// Double_t ratio = fNbg>0.0 ? fNex/(fNbg*fNormFactor) : 0.0;
418// sprintf(tx, " Significance = %6.2f, Nex/(Nbg*NormFactor) = %6.2f", fSigLiMa, ratio);
419// pt->AddText(tx);
420 // sprintf(tx, " Measured OFF background:");
421// TText *t7 = pt->AddText(tx);
422// t7->SetTextSize(0.02);
423// t7->SetTextColor(8);
424
425// sprintf(tx, "Results for |alpha|< %6.2f [\\circ] using fit to OFF:",
426// fFindsigON.GetAlphasi());
427// TText *t6 = pt->AddText(tx);
428// t6->SetTextSize(0.03);
429// t6->SetTextColor(8);
430
431 // sprintf(tx, "Noff = %6.1f", fNoff);
432// pt->AddText(tx);
433
434// sprintf(tx, "Nbg = %8.1f",
435// fNbgFit*fNormFactor);
436// pt->AddText(tx);
437
438// sprintf(tx, "Nex = %8.1f", fNexFit);
439// pt->AddText(tx);
440
441
442 sprintf(tx, "T_{obs} = 64 min.");
443 pt->AddText(tx);
444
445 sprintf(tx, "Zenith Angle = 7-17 #circ");
446 pt->AddText(tx);
447
448 sprintf(tx, "Significance = %6.2f #sigma", fSigLiMaFit);
449 TText *t8 =pt->AddText(tx);
450 //t8->SetTextSize(0.1);
451 t8->SetTextColor(8);
452 }
453
454 else if(opt.Contains("all",TString::kIgnoreCase))
455 {
456 pt = new TPaveText(0.20, 0.33, 0.88, 0.90, "NDC");
457
458 //
459 // Fit to OFF data
460 //
461 sprintf(tx, "Results of polynomial fit to OFF (order %2d) :", fDegree);
462 TText *t1 = pt->AddText(tx);
463 t1->SetTextSize(0.03);
464 t1->SetTextColor(2);
465
466 sprintf(tx, " (%6.2f< |alpha| <%6.2f [\\circ])", fAlphaminOFF,
467 fAlphamaxOFF);
468 pt->AddText(tx);
469
470 sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f",
471 fFindsigOFF.GetChisq(), fFindsigOFF.GetNdf(), fFindsigOFF.GetProb());
472 pt->AddText(tx);
473
474 sprintf(tx, " Nbgtot(fit)= %8.1f #pm %8.1f",
475 fFindsigOFF.GetNbgtotFitted(), fFindsigOFF.GetdNbgtotFitted());
476 pt->AddText(tx);
477
478 sprintf(tx, " Nbgtot(meas) = %8.1f", fFindsigOFF.GetNbgtot());
479 pt->AddText(tx);
480
481
482 //
483 // Results
484 //
485 sprintf(tx, "Results for |alpha|< %6.2f [\\circ] :",
486 fFindsigON.GetAlphasi());
487 TText *t6 = pt->AddText(tx);
488 t6->SetTextSize(0.03);
489 t6->SetTextColor(8);
490
491 sprintf(tx, " Non = %8.1f #pm %8.1f", fNon, fdNon);
492 pt->AddText(tx);
493
494
495 // PRINT INFORMATION ABOUT MEASURED QUANTITIES
496 sprintf(tx, " Measured OFF background:");
497 TText *t7 = pt->AddText(tx);
498 t7->SetTextSize(0.02);
499 t7->SetTextColor(8);
500
501 sprintf(tx, " Nbg(meas)*NormFactor = %8.1f #pm %8.1f",
502 fNbg*fNormFactor, fdNbg*fNormFactor);
503 pt->AddText(tx);
504
505 sprintf(tx, " Nex (meas) = %8.1f #pm %8.1f", fNex, fdNex);
506 pt->AddText(tx);
507
508 sprintf(tx," Normalization Factor (= Non/Noff) = %4.4f",fNormFactor);
509 pt->AddText(tx);
510
511 Double_t ratio = fNbg>0.0 ? fNex/(fNbg*fNormFactor) : 0.0;
512 sprintf(tx, " Significance(mea) = %6.2f, Nex/(Nbg*NormFactor) = %6.2f", fSigLiMa, ratio);
513 pt->AddText(tx);
514
515
516 // PRINT INFORMATION ABOUT FITTED QUANTITIES
517 sprintf(tx, " Fitted OFF background:");
518 TText *t8 = pt->AddText(tx);
519 t8->SetTextSize(0.02);
520 t8->SetTextColor(8);
521
522 sprintf(tx, " Nbg(fit)*NormFactor = %8.1f #pm %8.1f",
523 fNbgFit*fNormFactor, fdNbgFit*fNormFactor);
524 pt->AddText(tx);
525
526 sprintf(tx, " Nex(fit) = %8.1f #pm %8.1f", fNexFit, fdNexFit);
527 pt->AddText(tx);
528
529 sprintf(tx, " Effective Noff (i.e. fNoff) = %6.1f, Gamma = %4.4f",
530 fNoff, fGamma);
531 pt->AddText(tx);
532
533 ratio = fNbgFit>0.0 ? fNexFit/(fNbgFit*fNormFactor):0.0;
534 sprintf(tx, " Significance(fit) = %6.2f, Nex/(Nbg*NormFactor) = %6.2f", fSigLiMaFit, ratio);
535 pt->AddText(tx);
536
537
538 //
539 // Fit to ON data (polynomial+Gaus)
540 //
541 sprintf(tx, "Results of (polynomial+Gauss) fit to ON:");
542 TText *t9 = pt->AddText(tx);
543 t9->SetTextSize(0.03);
544 t9->SetTextColor(4);
545
546 //sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f",
547 //fGChisq, fGNdf, fGProb);
548 //pt->AddText(tx);
549
550 sprintf(tx, " Sigma of Gauss = %8.1f #pm %8.1f [\\circ]",
551 fFindsigON.GetSigmaGauss(), fFindsigON.GetdSigmaGauss());
552 pt->AddText(tx);
553
554 sprintf(tx, " total no.of excess events = %8.1f #pm %8.1f",
555 fFindsigON.GetNexGauss(), fFindsigON.GetdNexGauss());
556 pt->AddText(tx);
557 }
558 else
559 return;
560
561
562 pt->SetFillStyle(0);
563 pt->SetBorderSize(0);
564 pt->SetTextAlign(12);
565
566 pt->Draw();
567}
568
569
570// --------------------------------------------------------------------------
571//
572// Set flag fRebin
573//
574void MHMyFindSignificanceONOFF::SetRebin(Bool_t b)
575{
576 fRebin = b;
577
578 *fLog << "MHMyFindSignificanceONOFF::SetRebin; flag fRebin set to "
579 << (b? "kTRUE" : "kFALSE") << endl;
580}
581
582// --------------------------------------------------------------------------
583//
584// Set flag fReduceDegree
585//
586void MHMyFindSignificanceONOFF::SetReduceDegree(Bool_t b)
587{
588 fReduceDegree = b;
589
590 *fLog << "MHMyFindSignificanceONOFF::SetReduceDegree; flag fReduceDegree set to "
591 << (b? "kTRUE" : "kFALSE") << endl;
592}
593
594// ---------------------------------------------------------------------------
595//
596// Given an histogram, count the number of events within a given range.
597// Is equivaltent to TH1::Integral(binlo,binhi) but making sure that the bins
598// are fully contain in the interval (within a given tolerance).
599//
600Int_t MHMyFindSignificanceONOFF::CountEventsInRange(TH1* hist, Double_t alphalo, Double_t alphaup, Double_t* numevents, Double_t* numeventserror)
601{
602
603 const Int_t nbins = hist->GetNbinsX();
604 const Double_t binwidth = hist->GetBinWidth(1);
605
606 Double_t nevt = 0;
607 Double_t nevtErr = 0;
608
609 //
610 // Loop over all the bins
611 //
612 for(Int_t i=1; i<=nbins; i++)
613 {
614 Double_t xlo = hist->GetBinLowEdge(i);
615 Double_t xup = hist->GetBinLowEdge(i+1);
616
617 // bin must be completely contained in the signal region
618 if( xlo >= (alphalo-fEps) && xup <= (alphaup+fEps))
619 {
620 // check all bins all of same width
621 const Double_t width = fabs(xup-xlo);
622 if (fabs(width-binwidth) > fEps)
623 {
624 *fLog << err << GetName() << " alpha plot has variable binning, which is not allowed...exit." << endl;
625 return kFALSE;
626 }
627
628
629 nevt += hist->GetBinContent(i);
630 nevtErr += hist->GetBinError(i) * hist->GetBinError(i);
631 }
632 }
633
634 nevtErr = TMath::Sqrt(nevtErr);
635
636
637 // Return results
638 *numevents = nevt;
639 *numeventserror = nevtErr;
640}
641
642
643// ---------------------------------------------------------------------------
644//
645// Counts the number of events in the Bkg region in the ON_Data and OFF-Data
646// histograms. Then compute the normalization factor as:
647// NormFactor = NoffON / NoffNOFF
648//
649Int_t MHMyFindSignificanceONOFF::CalcNormalizationFactor(TH1* hon, TH1* hoff, Double_t alphalo, Double_t alphaup)
650{
651 Double_t nevton = 0;
652 Double_t nevtoff = 0;
653 Double_t dnevton = 0;
654 Double_t dnevtoff = 0;
655
656 // counts events in the region (alphalo,alphaup) for the ON-Data hist
657 CountEventsInRange(hon, alphalo, alphaup,&nevton,&dnevton);
658
659 // counts events in the region (alphalo,alphaup) for the OFF-Data hist
660 CountEventsInRange(hoff, alphalo, alphaup,&nevtoff,&dnevtoff);
661
662
663 // Normalization factor and its error
664 fNormFactor = nevton/nevtoff;
665
666 fNormFactorError = 1./nevton + 1./nevtoff;
667 fNormFactorError = TMath::Sqrt(fNormFactorError);
668 fNormFactorError *= fNormFactor; // <--------------
669
670
671 *fLog << "non = " << nevton << " noff = " << nevtoff << endl;
672 *fLog << "Normalization factor = " << fNormFactor << " +- "
673 << fNormFactorError << endl;
674
675 return kTRUE;
676}
677
678
679
680// --------------------------------------------------------------------------
681//
682// SigmaLiMa
683//
684// calculates the significance according to Li & Ma
685// ApJ 272 (1983) 317
686//
687Bool_t MHMyFindSignificanceONOFF::SigmaLiMa(Double_t non, Double_t noff,
688 Double_t gamma, Double_t *siglima)
689{
690 if (gamma <= 0.0 || non <= 0.0 || noff <= 0.0)
691 {
692 *siglima = 0.0;
693 return kFALSE;
694 }
695
696 Double_t help1 = non * log( (1.0+gamma)*non / (gamma*(non+noff)) );
697 Double_t help2 = noff * log( (1.0+gamma)*noff / ( non+noff ) );
698 *siglima = sqrt( 2.0 * (help1+help2) );
699
700 Double_t nex = non - gamma*noff;
701 if (nex < 0.0)
702 *siglima = - *siglima;
703
704 //*fLog << "MHMyFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = "
705 // << non << ", " << noff << ", " << gamma << ", " << *siglima << endl;
706
707 return kTRUE;
708}
709
710
711// --------------------------------------------------------------------------
712//
713// Sigma Vs Alpha
714//
715//
716Bool_t MHMyFindSignificanceONOFF::SigmaVsAlpha(TH1 *histON, TH1 *histOFF, Double_t alphamin, Double_t alphamax, Int_t degree, Bool_t draw)
717{
718
719
720 //--------------------------------------------
721 // Create histogram
722
723 //Int_t nsteps = 15;
724 Int_t nsteps = histON->FindBin(40);
725
726 TH1D* fSigVsAlpha = new TH1D("SigVsAlpha","Sigma vs Alpha", nsteps, 0.0, 30);
727 MH::SetBinning(fSigVsAlpha, histON->GetXaxis());
728 fSigVsAlpha->SetXTitle("upper edge of signal region in |alpha| [\\circ]");
729 fSigVsAlpha->SetYTitle("Significance of gamma signal");
730
731
732 TGraph* grSigVsAlpha = new TGraph;
733 TGraph* grNexVsAlpha = new TGraph;
734
735
736
737
738 //--------------------------------------------
739 // loop over different signal regions
740 for (Int_t i=1; i<=nsteps; i++)
741 {
742 //Double_t alphasig = fSigVsAlpha->GetBinCenter(i);
743 //Double_t alphasig = fSigVsAlpha->GetBinLowEdge(i+1);
744 Double_t alphasig = histON->GetBinLowEdge(i+1);
745
746
747 gLog.SetNullOutput(1);
748 //FindSigmaONOFF(histON, histOFF, alphasig, alphamin, alphamax, degree);
749 MHMyFindSignificanceONOFF findsig;
750 findsig.SetRebin(fRebin);
751 findsig.SetReduceDegree(fReduceDegree);
752
753 findsig.FindSigmaONOFF( histON, histOFF, alphasig, alphamin, alphamax, degree);
754 const Double_t sig = findsig.GetSigLiMa();
755 const Double_t Nex = findsig.GetNex();
756
757 fSigVsAlpha->SetBinContent(i, sig);
758
759 if( sig>fBestSigma )
760 {
761 fBestSigma = sig;
762 fBestAlphaCut = alphasig;
763 fBestExcess = Nex;
764 }
765
766 grSigVsAlpha->SetPoint(i-1,alphasig,sig);
767 grNexVsAlpha->SetPoint(i-1,alphasig,Nex);
768
769 cout << "****************** " << i<< " " << alphasig << " " << sig<< " " << Nex << "*************"<< endl;
770 gLog.SetNullOutput(0);
771 }
772
773
774// fSigVsAlpha->DrawCopy();
775
776
777 if (!draw)
778 return kTRUE;
779
780 // -----------------------------------------------------------------
781 //
782 // Draw
783 //
784 TCanvas *c = new TCanvas("SigVsAlpha", "Sigma vs Alpha", 600, 600);
785 c->SetGridx();
786 c->SetGridy();
787
788 Double_t margin = 0.12; // give more space for drawing the axis titles
789 c->SetLeftMargin(margin);
790 c->SetRightMargin(margin);
791
792
793 //
794 // Draw Sig Vs Alpha
795 //
796 grSigVsAlpha->SetMarkerStyle(8);
797 //grSigVsAlpha->SetMarkerColor(2);
798 //grSigVsAlpha->SetLineColor(2);
799 grSigVsAlpha->Draw("apl");
800 grSigVsAlpha->GetHistogram()->SetXTitle("Alpha Cut [\\circ]");
801 grSigVsAlpha->GetHistogram()->SetYTitle("Significance");
802 grSigVsAlpha->GetHistogram()->GetYaxis()->SetTitleOffset(1.3);
803
804
805 //
806 // Draw second graph on a transparent pad
807 //
808 TPad *overlay = new TPad("overlay","",0,0,1,1);
809 overlay->SetFillStyle(4000);
810 overlay->SetFillColor(0);
811 overlay->SetFrameFillStyle(4000);
812 overlay->SetLeftMargin(margin);
813 overlay->SetRightMargin(margin);
814 overlay->Draw();
815 overlay->cd();
816
817 grNexVsAlpha->SetMarkerStyle(8);
818 grNexVsAlpha->SetMarkerColor(kRed);
819 grNexVsAlpha->SetLineColor(kRed);
820 grNexVsAlpha->Draw("apl");
821
822 // Trick to not draw the graph Y-Axis
823 grNexVsAlpha->GetHistogram()->GetYaxis()->SetNdivisions(0);
824
825
826 //
827 // Draw an axis on the right side
828 //
829 c->Update(); //<-- Imprescindible for updating the gPad coordinates
830
831 TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
832 gPad->GetUxmax(), gPad->GetUymax(),
833 gPad->GetUymin(),gPad->GetUymax() ,510,"+L");
834 axis->SetTitle("Excess events");
835 axis->SetTitleOffset(1.3);
836 //axis->SetLineColor(kRed);
837 axis->SetLabelColor(kRed);
838 axis->SetTextColor(kRed);
839 axis->Draw();
840
841
842 c->Modified();
843 c->Update();
844
845 return kTRUE;
846}
Note: See TracBrowser for help on using the repository browser.