source: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MHFindSignificanceONOFF.cc@ 6724

Last change on this file since 6724 was 6355, checked in by mazin, 20 years ago
*** empty log message ***
File size: 95.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): Wolfgang Wittek, July 2003 <mailto:wittek@mppmu.mpg.de>
19! David Paneque, Nov 2003 <mailto:dpaneque@mppmu.mpg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MHFindSignificanceONOFF
29//
30// determines the significance of a gamma signal in an |alpha| plot
31// it uses real OFF events in the computation of excess events
32//
33// Input : 2 TH1 histogram (ON and OFF) of |alpha| : with 0 < |alpha| < 90 degrees
34//
35//************TEMP ************************************************
36
37
38// alphamin, alphamax : defining the background region
39// alphasig : defining the signal region for which
40// the significance is calculated
41// degree : the degree of the polynomial to be fitted to the background
42// ( a0 + a1*x + a2*x**2 + a3*x**3 + ...)
43//
44// Output :
45//
46//
47// - the number of events in the signal region (Non)
48// the number of background events in the signal region (Nbg)
49// (counted number of events 'NoffSig' and fitted number of events 'NoffSigFitted
50// - the number of excess events in the signal region (Nex = Non - Nbg)
51// (again, counted 'NexONOFF' and fitted 'NexONOFFFitted'
52// - thew effective number of background events (Noff), and gamma :
53// Nbg = gamma * Noff
54// -
55// - the significance of the gamma signal according to Li & Ma
56// 3 significances are computed
57// a) LiMa formula (17) using fitted quantities; fSigLiMa
58// b) LiMa formula (17) using counted quantities; fSigLiMa2
59// c) LiMa formula (5) using counted quantities.
60//
61// call member function 'FindSigmaONOFF'
62// to fit the background and to determine the significance
63//
64//
65//
66/////////////////////////////////////////////////////////////////////////////
67
68////////// *********** ENDTEMPORAL *************************
69
70
71#include "MHFindSignificanceONOFF.h"
72
73#include <fstream>
74#include <math.h>
75
76#include <TArrayD.h>
77#include <TArrayI.h>
78#include <TH1.h>
79#include <TF1.h>
80#include <TCanvas.h>
81#include <TFitter.h>
82#include <TMinuit.h>
83#include <TPaveText.h>
84#include <TStyle.h>
85#include <TPostScript.h>
86
87#include "MLog.h"
88#include "MLogManip.h"
89#include "MMinuitInterface.h"
90
91
92ClassImp(MHFindSignificanceONOFF);
93
94using namespace std;
95
96const TString MHFindSignificanceONOFF::gsDefName = "MHFindSignificanceONOFF";
97const TString MHFindSignificanceONOFF::gsDefTitle = "Find Significance in alpha plot";
98
99
100
101// --------------------------------------------------------------------------
102//
103// fcnpoly
104//
105// calculates the chi2 for the fit of the polynomial function 'poly'
106// to the histogram 'fhist'
107//
108// it is called by CallMinuit() (which is called in FitPolynomial())
109//
110// bins of fhist with huge errors are ignored in the calculation of the chi2
111// (the huge errors were set in 'FitPolynomial()')
112//
113
114static void fcnpoly(Int_t &npar, Double_t *gin, Double_t &f,
115 Double_t *par, Int_t iflag)
116{
117 TH1 *fhist = (TH1*)gMinuit->GetObjectFit();
118 TF1 *fpoly = fhist->GetFunction("Poly");
119
120
121 //-------------------------------------------
122
123 Double_t chi2 = 0.0;
124
125 Int_t nbins = fhist->GetNbinsX();
126 Int_t mbins = 0;
127 for (Int_t i=1; i<=nbins; i++)
128 {
129 Double_t content = fhist->GetBinContent(i);
130 Double_t error = fhist->GetBinError(i);
131 Double_t center = fhist->GetBinCenter(i);
132
133 //-----------------------------
134 // ignore unwanted points if (error > 1.e19)
135 continue;
136
137 if (content <= 0.0)
138 {
139 gLog << "fcnpoly : bin with zero content; i, content, error = "
140 << i << ", " << content << ", " << error << endl;
141 continue;
142 }
143
144 if (error <= 0.0)
145 {
146 gLog << "fcnpoly : bin with zero error; i, content, error = "
147 << i << ", " << content << ", " << error << endl;
148 continue;
149 }
150
151 //-----------------------------
152 mbins++;
153
154 Double_t fu;
155 fu = fpoly->EvalPar(&center, par);
156
157 // the fitted function must not be negative
158 if (fu <= 0.0)
159 {
160 chi2 = 1.e10;
161 break;
162 }
163
164 Double_t temp = (content - fu) / error;
165 chi2 += temp*temp;
166 }
167
168 //-------------------------------------------
169
170 f = chi2;
171
172 //-------------------------------------------
173 // final calculations
174 //if (iflag == 3)
175 //{
176 //}
177
178 //-------------------------------------------------------------
179}
180
181
182// --------------------------------------------------------------------------
183//
184// fcnpolyOFF
185//
186// calculates the chi2 for the fit of the polynomial function 'poly'
187// to the histogram 'fhist'
188//
189// it is called by CallMinuit() (which is called in FitPolynomial())
190//
191// bins of fhist with huge errors are ignored in the calculation of the chi2
192// (the huge errors were set in 'FitPolynomial()')
193//
194
195static void fcnpolyOFF(Int_t &npar, Double_t *gin, Double_t &f,
196 Double_t *par, Int_t iflag)
197{
198 TH1 *fhist = (TH1*)gMinuit->GetObjectFit();
199 TF1 *fpoly = fhist->GetFunction("PolyOFF");
200
201
202 //-------------------------------------------
203
204 Double_t chi2 = 0.0;
205
206 Int_t nbins = fhist->GetNbinsX();
207 Int_t mbins = 0;
208 for (Int_t i=1; i<=nbins; i++)
209 {
210 Double_t content = fhist->GetBinContent(i);
211 Double_t error = fhist->GetBinError(i);
212 Double_t center = fhist->GetBinCenter(i);
213
214 //-----------------------------
215 // ignore unwanted points
216 if (error > 1.e19)
217 continue;
218
219 if (content <= 0.0)
220 {
221 gLog << "fcnpoly : bin with zero content; i, content, error = "
222 << i << ", " << content << ", " << error << endl;
223 continue;
224 }
225
226 if (error <= 0.0)
227 {
228 gLog << "fcnpoly : bin with zero error; i, content, error = "
229 << i << ", " << content << ", " << error << endl;
230 continue;
231 }
232
233 //-----------------------------
234 mbins++;
235
236 Double_t fu;
237 fu = fpoly->EvalPar(&center, par);
238
239 // the fitted function must not be negative
240 if (fu <= 0.0)
241 {
242 chi2 = 1.e10;
243 break;
244 }
245
246 Double_t temp = (content - fu) / error;
247 chi2 += temp*temp;
248 }
249
250 //-------------------------------------------
251
252 f = chi2;
253
254 //-------------------------------------------
255 // final calculations
256 //if (iflag == 3)
257 //{
258 //}
259
260 //-------------------------------------------------------------
261}
262
263
264
265
266// --------------------------------------------------------------------------
267//
268// fcnpolygauss
269//
270// calculates the chi2 for the fit of the (polynomial+Gauss) function
271// 'PolyGauss' to the histogram 'fhist'
272//
273// it is called by CallMinuit() (which is called in FitGaussPoly())
274//
275// bins of fhist with huge errors are ignored in the calculation of the chi2
276// (the huge errors were set in 'FitGaussPoly()')
277//
278
279static void fcnpolygauss(Int_t &npar, Double_t *gin, Double_t &f,
280 Double_t *par, Int_t iflag)
281{
282 TH1 *fhist = (TH1*)gMinuit->GetObjectFit();
283 TF1 *fpolygauss = fhist->GetFunction("PolyGauss");
284
285
286 //-------------------------------------------
287
288 Double_t chi2 = 0.0;
289
290 Int_t nbins = fhist->GetNbinsX();
291 Int_t mbins = 0;
292 for (Int_t i=1; i<=nbins; i++)
293 {
294 Double_t content = fhist->GetBinContent(i);
295 Double_t error = fhist->GetBinError(i);
296 Double_t center = fhist->GetBinCenter(i);
297
298 //-----------------------------
299 // ignore unwanted points
300 if (error > 1.e19)
301 continue;
302
303 if (content <= 0.0)
304 {
305 gLog << "fcnpolygauss : bin with zero content; i, content, error = "
306 << i << ", " << content << ", " << error << endl;
307 continue;
308 }
309
310 if (error <= 0.0)
311 {
312 gLog << "fcnpolygauss : bin with zero error; i, content, error = "
313 << i << ", " << content << ", " << error << endl;
314 continue;
315 }
316
317 //-----------------------------
318 mbins++;
319
320 Double_t fu;
321 fu = fpolygauss->EvalPar(&center, par);
322
323 // the fitted function must not be negative
324 if (fu <= 0.0)
325 {
326 chi2 = 1.e10;
327 break;
328 }
329
330 Double_t temp = (content - fu) / error;
331 chi2 += temp*temp;
332 }
333
334 //-------------------------------------------
335
336 f = chi2;
337
338 //-------------------------------------------
339 // final calculations
340 //if (iflag == 3)
341 //{
342 //}
343
344 //-------------------------------------------------------------
345}
346
347
348
349// --------------------------------------------------------------------------
350//
351// Constructor
352//
353MHFindSignificanceONOFF::MHFindSignificanceONOFF(const char *name, const char *title)
354{
355 fName = name ? name : gsDefName.Data();
356 fTitle = title ? title : gsDefTitle.Data();
357
358// fSigVsAlpha = NULL;
359
360 fPoly = NULL;
361 fPolyOFF = NULL;
362 fPolyOFFNormalized = NULL;
363 fGPoly = NULL;
364 fGBackg = NULL;
365
366 fHist = NULL;
367 fHistOrig = NULL;
368
369 fHistOFF = NULL;
370 fHistOrigOFF = NULL;
371 fHistOFFNormalized = NULL;
372
373 // allow rebinning of the alpha plot
374// fRebin = kTRUE;
375
376 fRebin = kFALSE;
377
378 // allow reducing the degree of the polynomial
379 fReduceDegree = kTRUE;
380
381 // Low and Upper limits for the OFF alpha distribution fit
382 // are set to 0 and 90 degrees respectively
383
384 fAlphaminOFF = 0.0;
385 fAlphamaxOFF = 90.0;
386
387 // use quantities computed from the fits
388 // The variable allows the user to NOT use these quantities when there is not
389 // enough statistics and fit not always is possible.
390 // Default value is kTRUE.
391 fUseFittedQuantities = kTRUE;
392
393 // Bool variable used to decide wether to print or not the results
394 // of the fit, significance, Nex... onto the final alpha plot.
395 // for the time being, this variable is set in the constructor.
396 // At some point, I might make it such it can be set externally...
397
398 fPrintResultsOntoAlphaPlot = kTRUE;
399 //fPrintResultsOntoAlphaPlot = kFALSE;
400
401
402
403 fCanvas = NULL;
404
405 fSavePlots = kFALSE; // By default plots are not saved in Postscriptfiles
406
407 fPsFilename = NULL;
408
409
410
411}
412
413// --------------------------------------------------------------------------
414//
415// Destructor.
416//
417// =====> it is not clear why one obtains sometimes a segmentation violation
418// when the destructor is active <=======================
419//
420// therefore the 'return'statement
421//
422
423MHFindSignificanceONOFF::~MHFindSignificanceONOFF()
424{
425
426
427
428 return;
429
430 *fLog << "destructor of MHFindSignificanceONOFF is called" << endl;
431
432
433 fPsFilename = NULL;
434
435 cout << "PS set to null... " << endl;
436
437 delete fHist;
438 delete fHistOFF;
439
440 delete fHistOFFNormalized;
441
442 cout << "Histos removed..." << endl;
443
444 delete fPoly;
445 delete fPolyOFF;
446 delete fPolyOFFNormalized;
447 delete fGPoly;
448 delete fGBackg;
449
450 cout << "Functions are also removed..." << endl;
451
452 // delete fCanvas; if I removed fCanvas pointed memory address the
453 // program crashes ???
454
455 *fLog << "destructor of MHFindSignificanceONOFF finished successfully" << endl;
456
457
458}
459
460// --------------------------------------------------------------------------
461//
462// Set flag fRebin
463//
464// if flag is kTRUE rebinning of the alpha plot is allowed
465//
466//
467void MHFindSignificanceONOFF::SetRebin(Bool_t b)
468{
469 fRebin = b;
470
471 *fLog << "MHFindSignificanceONOFF::SetRebin; flag fRebin set to "
472 << (b? "kTRUE" : "kFALSE") << endl;
473}
474
475// --------------------------------------------------------------------------
476//
477// Set flag fReduceDegree
478//
479// if flag is kTRUE reducing of the degree of the polynomial is allowed
480//
481//
482void MHFindSignificanceONOFF::SetReduceDegree(Bool_t b)
483{
484 fReduceDegree = b;
485
486 *fLog << "MHFindSignificanceONOFF::SetReduceDegree; flag fReduceDegree set to "
487 << (b? "kTRUE" : "kFALSE") << endl;
488}
489
490// Function that returns one of the 3 LiMa sigmas.
491// The returned value is the one used in the optimization
492// and final alpha plots.
493
494// For the time being, if fUseFittedQuantities = kTRUE (default)
495// fSigLiMa is used, otherwise, fSigLiMa2 is used.
496
497Double_t MHFindSignificanceONOFF::GetSignificance()
498{
499
500 if(fUseFittedQuantities)
501 {
502 return fSigLiMa;
503 }
504
505 return fSigLiMa2;
506
507
508
509}
510
511
512// --------------------------------------------------------------------------
513//
514// FindSigmaONOFF
515//
516// calls FitPolynomialOFF it gets bkg events in "signal" region
517// from histogram histOFF which is the alpha
518// distribution of OFF data NON normalized.
519// Normalization factor is also one of the
520// arguments.
521// calls DetExcess to determine the number of excess events
522// using the previously computed bkg events
523
524// calls SigmaLiMa to determine the significance of the gamma signal
525// in the range |alpha| < alphasig
526// calls FitGaussPoly to fit a (polynomial+Gauss) function in the
527// whole |alpha| region of ON - OFF diagram
528//
529//
530
531Bool_t MHFindSignificanceONOFF::FindSigmaONOFF(TH1 *fhistON,
532 TH1 *fhistOFF,
533 Double_t NormFactor,
534 Double_t alphamin,
535 Double_t alphamax,
536 Int_t degreeON,
537 Int_t degreeOFF,
538 Double_t alphasig,
539 Bool_t drawpoly,
540 Bool_t fitgauss,
541 Bool_t print,
542 Bool_t saveplots,
543 //TPostScript* PsFile
544 const TString psfilename)
545{
546
547 //*fLog << "MHFindSignificanceONOFF::FindSigma;" << endl;
548
549 // Pointer to object TPostScript where plots will be stored
550 // is copied into member varialbe
551 // NOT WORKING !!!
552 // fPsFilename = PsFile;
553
554 // Temporally (while TPostSctipt option is not working) psfiles
555 // will be produced by the standard way (Canvas.SaveAs())
556
557 fPsFilenameString = psfilename;
558
559
560
561 // "3 Li and Ma significances" are initialized to 0.0
562
563 fSigLiMa = 0.0;
564 fSigLiMa2 = 0.0;
565 fSigLiMa3 = 0.0;
566
567
568 // fNormFactor is set
569
570 fNormFactor = NormFactor;
571
572 // Report when this histograms given in the
573 // arguments are empty
574
575 Double_t tmpdouble= -1.0;
576
577 tmpdouble = double(fhistON->GetEntries());
578 if (tmpdouble < 0.5)
579 {
580 cout << "MHFindSignificanceONOFF::FindSigmaONOFF; ERROR " << endl
581 << "fhistON has ZERO entries" << endl;
582 }
583
584 tmpdouble = double(fhistOFF->GetEntries());
585 if (tmpdouble < 0.5)
586 {
587 cout << "MHFindSignificanceONOFF::FindSigmaONOFF; ERROR " << endl
588 << "fhistOFF has ZERO entries" << endl;
589 }
590
591
592
593
594 // Variables set for alpha from ON events
595 fHistOrig = fhistON;
596 fHist = (TH1*)fHistOrig->Clone();
597 fHist->SetName(fhistON->GetName());
598 fDegree = degreeON;
599
600
601 // Variables set for alpha from OFF events
602
603 fHistOrigOFF = fhistOFF;
604 fHistOFF = (TH1*)fHistOrigOFF->Clone();
605 fHistOFF->SetName(fhistOFF->GetName());
606 fDegreeOFF = degreeOFF;
607
608 if ( !fHist )
609 {
610 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; Clone of histogram could not be generated"
611 << endl;
612 return kFALSE;
613 }
614
615
616
617 if ( !fHistOFF )
618 {
619 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; Clone of OFF histogram could not be generated"
620 << endl;
621 return kFALSE;
622 }
623
624
625
626
627 fHist->Sumw2();
628 fHist->SetXTitle("|alpha| [\\circ]");
629 fHist->SetYTitle("Counts");
630 fHist->UseCurrentStyle();
631
632
633 fHistOFF->Sumw2(); // if error has been set (via function SetBinError(j))
634 // the errors set remain, i.e. are not overwritten with the sum of the square of weights.
635 // Which means that this function will not have any effect.
636
637 fHistOFF->SetXTitle("|alpha| [\\circ]");
638 fHistOFF->SetYTitle("Counts");
639 fHistOFF->UseCurrentStyle();
640
641
642
643
644
645/////////////////////////////////////
646
647 fAlphamin = alphamin;
648 fAlphamax = alphamax;
649 fAlphammm = (alphamin+alphamax)/2.0;
650 fAlphasig = alphasig;
651
652
653
654
655 // UYpper limits for fit to OFF data set are taken also from alphamax
656 // fAlphaminOFF is set in constructor. Inn principle, it WILL ALWAYS BE ZERO.
657
658 fAlphamaxOFF = alphamax;
659
660
661 fDraw = drawpoly;
662 fSavePlots = saveplots;
663 fFitGauss = fitgauss;
664
665
666 //--------------------------------------------
667 // fit a polynomial in the backgr region
668
669 //*fLog << "MHFindSignificanceONOFF::FindSigma; calling FitPolynomial()" << endl;
670 if ( !FitPolynomialOFF())
671 {
672 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; PolynomialOFF failed"
673 << endl;
674 return kFALSE;
675 }
676
677
678 //--------------------------------------------
679 // calculate the number of excess events in the signal region
680
681 //*fLog << "MHFindSignificanceONOFF::FindSigma; calling DetExcess()" << endl;
682 if ( !DetExcessONOFF())
683 {
684 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; DetExcessONOFF failed"
685 << endl;
686 return kFALSE;
687 }
688
689
690
691
692
693 //--------------------------------------------
694 // calculate the significance of the excess
695
696 //*fLog << "MHFindSignificanceONOFF::FindSigma; calling SigmaLiMa()" << endl;
697
698 // For testing purposes "3 Li&Ma significances" will be computed
699 // At some point, only one will remain
700
701
702 Double_t siglima = 0.0;
703
704
705 // Significance computed using effective number of OFF events in signal
706 // region (fNoff) and gamma factor (fGama).
707 // This is Wolfgang approach to the calulation of significance
708 // using Li&Ma formula and estimated OFF events from polynomial fit.
709
710 if(fUseFittedQuantities)
711 {
712 if ( !SigmaLiMa(fNon, fNoff, fGamma, &siglima) )
713 {
714 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; SigmaLiMa failed"
715 << endl;
716 return kFALSE;
717 }
718
719 fSigLiMa = siglima;
720 }
721 else
722 {
723 fSigLiMa = 0.0;
724 }
725
726
727
728 // Significance computed using counted number of OFF events in signal
729 // region (fNoffSig) and normalization factor (fNormFactor).
730 // This is the strictly speaking the significance in Li&Ma paper...
731
732
733 if ( !SigmaLiMa(fNon, fNoffSig, fNormFactor, &siglima) )
734 {
735 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; SigmaLiMa2 failed"
736 << endl;
737 return kFALSE;
738 }
739 fSigLiMa2 = siglima;
740
741
742 // Significance computed using counted number of OFF events in signal
743 // region (fNoffSig) and normalization factor (fNormFactor).
744 // significance of gamma signal according to Li & Ma using
745 // formula (5)
746
747
748
749 if ( !SigmaLiMaForm5(fNon, fNoffSig, fNormFactor, &siglima) )
750 {
751 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; SigmaLiMa failed"
752 << endl;
753 return kFALSE;
754 }
755
756 fSigLiMa3 = siglima;
757
758
759
760
761
762//--------------------------------------------
763// calculate the error of the number of excess events
764// using fitted quantities and counted quantities
765
766
767// from fit to OFF histogram
768
769 if (fSigLiMa > 0.0)
770 {fdNexONOFFFitted = fNexONOFFFitted / fSigLiMa;}
771 else
772 {fdNexONOFFFitted = 0.0;}
773
774
775 // from counted OFF events
776 if (fSigLiMa2 > 0.0)
777 {
778 fdNexONOFF = fNexONOFF / fSigLiMa2;
779 }
780 else
781 {
782 fdNexONOFF = 0.0;
783 }
784
785 if (fDraw || fFitGauss)
786 {
787
788 // ------------------------------------------------
789 // Polynomial fit to bkg region from ON data is performed,
790 // Because some of the quantities will be used as
791 // initial values in the PolyGAuss fit.
792
793 // Besides, this function will modify the binning of fHist
794 // so that the number of events in each bin is big enough
795
796 // This might change in future...
797
798 if ( !FitPolynomial())
799 {
800 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; Polynomial fit failed"
801 << endl;
802 return kFALSE;
803 }
804 }
805
806
807 if (fDraw)
808 {
809
810 // Compute fHistOFFNormalized
811
812 if (!ComputeHistOFFNormalized())
813 {
814 *fLog << "MHFindSignificanceONOFF::ComputeHistOFFNormalized; Normalization of fHistOFF was not possible"
815 << endl;
816 return kFALSE;
817
818 }
819 }
820
821
822
823
824
825 //--------------------------------------------
826
827 //*fLog << "MHFindSignificanceONOFF::FindSigma; calling PrintPoly()" << endl;
828 if (print)
829 PrintPolyOFF();
830
831
832
833
834
835
836
837
838
839 //--------------------------------------------
840 // fit a (polynomial + Gauss) function to the ON-OFF alpha distribution
841
842 if (fFitGauss)
843 {
844
845 //--------------------------------------------------
846 // delete objects from this fit
847 // in order to have independent starting conditions for the next fit
848
849 delete gMinuit;
850 gMinuit = NULL;
851 //--------------------------------------------------
852
853 //*fLog << "MHFindSignificanceONOFF::FindSigma; calling FitGaussPoly()"
854 // << endl;
855 if ( !FitGaussPoly() )
856 {
857 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; FitGaussPoly failed"
858 << endl;
859 return kFALSE;
860 }
861
862 if (print)
863 {
864 //*fLog << "MHFindSignificanceONOFF::FindSigma; calling PrintPolyGauss()"
865 // << endl;
866 PrintPolyGauss();
867 }
868 }
869
870 //--------------------------------------------------
871 // draw the histogram if requested
872
873 if (fDraw)
874 {
875 // TEMPORALLY I will plot fHistOFF and fHistOFFNormalized (with the fits)
876
877
878 if (!DrawHistOFF())
879 {
880 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; Drawing of fHistOFF was not possible"
881 << endl;
882 // return kFALSE;
883
884 }
885
886
887 if (!DrawHistOFFNormalized())
888 {
889 *fLog << "MHFindSignificanceONOFF::DrawHistOFFNormalized; Drawing of fHistOFFNormalized was not possible"
890 << endl;
891 // return kFALSE;
892
893 }
894
895
896
897
898 //*fLog << "MHFindSignificanceONOFF::FindSigma; calling DrawFit()" << endl;
899 if ( !DrawFit() )
900 {
901 *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; DrawFit failed"
902 << endl;
903 return kFALSE;
904 }
905 }
906
907
908 //--------------------------------------------------
909 // delete objects from this fit
910 // in order to have independent starting conditions for the next fit
911
912 delete gMinuit;
913 gMinuit = NULL;
914 //--------------------------------------------------
915
916 return kTRUE;
917
918}
919
920
921//// ********************* end sigmaonoff ************************
922
923
924
925
926
927
928
929// UNDER CONSTRUCTION
930
931Bool_t MHFindSignificanceONOFF::SigmaVsAlphaONOFF(TH1 *fhistON, TH1 *fhistOFF,
932 Double_t alphamin, Double_t alphamax,
933 Int_t degree, Bool_t print)
934{
935 *fLog << " MHFindSignificanceONOFF::SigmaVsAlphaONOFF still under construction !!!" << endl;
936
937 return kFALSE;
938
939}
940
941
942
943
944// --------------------------------------------------------------------------
945//
946// FitPolynomialOFF
947// - create a clone 'fHistOFF' of the |alpha| distribution 'fHistOrigOFF'
948// - fit a polynomial of degree 'fDegreeOFF' to the alpha distribution
949// 'fHistOFF' in the region alphaminOFF < |alpha| < alphamaxOFF
950
951
952
953Bool_t MHFindSignificanceONOFF::FitPolynomialOFF()
954
955{
956 //--------------------------------------------------
957 // check the histogram :
958 // - calculate initial values of the parameters
959 // - check for bins with zero entries
960 // - set minimum errors
961 // - save the original errors
962 // - set errors huge outside the fit range
963 // (in 'fcnpolyOFF' points with huge errors will be ignored)
964
965
966 Double_t dummy = 1.e20;
967
968 Double_t mean;
969 Double_t rms;
970 Double_t nclose;
971 Double_t nfar;
972 Double_t a2init = 0.0;
973 TArrayD saveError;
974
975 Int_t nbins;
976 Int_t nrebin = 1;
977
978
979 //---------------- start while loop for rebinning -----------------
980 while(1)
981 {
982
983 fNzeroOFF = 0;
984 fMbinsOFF = 0;
985 fMlowOFF = 0;
986 fNoffTot = 0.0;
987
988 // same variables (as in fitpoly to ON data) are used
989 // here for naming the actual values for low/up limit for fit
990 fAlphami = 10000.0;
991 fAlphamm = 10000.0;
992 fAlphama = -10000.0;
993
994 mean = 0.0;
995 rms = 0.0;
996 nclose = 0.0;
997 nfar = 0.0;
998
999 nbins = fHistOFF->GetNbinsX();
1000 saveError.Set(nbins);
1001
1002 for (Int_t i=1; i<=nbins; i++)
1003 {
1004 saveError[i-1] = fHistOFF->GetBinError(i);
1005
1006 // bin should be completely contained in the fit range
1007 // (fAlphaminOFF, fAlphamaxOFF)
1008 Double_t xlo = fHistOFF->GetBinLowEdge(i);
1009 Double_t xup = fHistOFF->GetBinLowEdge(i+1);
1010
1011 if ( xlo >= fAlphaminOFF-fEps && xlo <= fAlphamaxOFF+fEps &&
1012 xup >= fAlphaminOFF-fEps && xup <= fAlphamaxOFF+fEps )
1013 {
1014 fMbinsOFF++;
1015
1016 if ( xlo < fAlphami )
1017 fAlphami = xlo;
1018
1019 if ( xup > fAlphama )
1020 fAlphama = xup;
1021
1022 Double_t content = fHistOFF->GetBinContent(i);
1023 // fNoffTot += content;
1024
1025 mean += content;
1026 rms += content*content;
1027
1028// count events in low-alpha and high-alpha region
1029 if ( xlo >= fAlphammm-fEps && xup >= fAlphammm-fEps)
1030 {
1031 nfar += content;
1032 if ( xlo < fAlphamm )
1033 fAlphamm = xlo;
1034 if ( xup < fAlphamm )
1035 fAlphamm = xup;
1036 }
1037 else
1038 {
1039 nclose += content;
1040 if ( xlo > fAlphamm )
1041 fAlphamm = xlo;
1042 if ( xup > fAlphamm )
1043 fAlphamm = xup;
1044 }
1045
1046 // count bins with zero entry
1047 if (content < 0.0001)
1048 fNzeroOFF++;
1049
1050 // set minimum error
1051 // modified by K. Mase
1052 if (content < 0.0)
1053 {
1054 fMlowOFF += 1;
1055 fHistOFF->SetBinError(i, 3.0);
1056 }
1057
1058 //*fLog << "Take : i, content, error = " << i << ", "
1059 // << fHist->GetBinContent(i) << ", "
1060 // << fHist->GetBinError(i) << endl;
1061
1062 continue;
1063 }
1064 // bin is not completely contained in the fit range : set error huge
1065
1066 fHistOFF->SetBinError(i, dummy);
1067
1068 //*fLog << "Omit : i, content, error = " << i << ", "
1069 // << fHist->GetBinContent(i) << ", " << fHist->GetBinError(i)
1070 // << endl;
1071
1072 }
1073
1074 // mean of entries/bin in the fit range
1075 if (fMbinsOFF > 0)
1076 {
1077 mean /= ((Double_t) fMbinsOFF);
1078 rms /= ((Double_t) fMbinsOFF);
1079 rms = sqrt( rms - mean*mean );
1080 }
1081
1082
1083
1084 // if there are no events in the background region
1085 // there is no reason for rebinning
1086 // and this is the condition for assuming a constant background (= 0)
1087 if (mean <= 0.0)
1088 break;
1089
1090 Double_t helpmi = fAlphami*fAlphami*fAlphami;
1091 Double_t helpmm = fAlphamm*fAlphamm*fAlphamm;
1092 Double_t helpma = fAlphama*fAlphama*fAlphama;
1093 Double_t help = (helpma-helpmm) * (fAlphamm-fAlphami)
1094 - (helpmm-helpmi) * (fAlphama-fAlphamm);
1095 if (help != 0.0)
1096 a2init = ( (fAlphamm-fAlphami)*nfar - (fAlphama-fAlphamm)*nclose )
1097 * 1.5 * fHistOFF->GetBinWidth(1) / help;
1098 else
1099 a2init = 0.0;
1100
1101
1102 //--------------------------------------------
1103 // rebin the histogram
1104 // - if a bin has no entries
1105 // - or if there are too many bins with too few entries
1106 // - or if the new bin width would exceed half the size of the
1107 // signal region
1108
1109 if ( !fRebin ||
1110 ( fNzeroOFF <= 0 && (Double_t)fMlowOFF<0.05*(Double_t)fMbinsOFF ) ||
1111 (Double_t)(nrebin+1)/(Double_t)nrebin * fHistOFF->GetBinWidth(1)
1112 > fAlphasig/2.0 )
1113 {
1114 //*fLog << "before break" << endl;
1115 break;
1116 }
1117
1118 nrebin += 1;
1119 TString histname = fHistOFF->GetName();
1120 delete fHistOFF;
1121 fHistOFF = NULL;
1122
1123 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; rebin the |alpha|OFF plot, grouping "
1124 << nrebin << " bins together" << endl;
1125
1126 // TH1::Rebin doesn't work properly
1127 //fHist = fHistOrig->Rebin(nrebin, "Rebinned");
1128 // use private routine RebinHistogram()
1129 fHistOFF = new TH1F;
1130 fHistOFF->Sumw2();
1131 fHistOFF->SetNameTitle(histname, histname);
1132 fHistOFF->UseCurrentStyle();
1133
1134 // do rebinning such that x0 remains a lower bin edge
1135 Double_t x0 = 0.0;
1136 if ( !RebinHistogramOFF(x0, nrebin) )
1137 {
1138 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; RebinHistgramOFF() failed"
1139 << endl;
1140 return kFALSE;
1141 }
1142
1143 fHistOFF->SetXTitle("|alpha| [\\circ]");
1144 fHistOFF->SetYTitle("Counts");
1145
1146 }
1147 //---------------- end of while loop for rebinning -----------------
1148
1149
1150 // if there are still too many bins with too few entries don't fit
1151 // and assume a constant background
1152
1153 fConstantBackg = kFALSE;
1154 if ( fNzeroOFF > 0 || (Double_t)fMlowOFF>0.05*(Double_t)fMbinsOFF )
1155 {
1156 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; polynomial fit not possible, fNzeroOFF, fMlowOFF, fMbinsOFF = "
1157 << fNzeroOFF << ", " << fMlowOFF << ", " << fMbinsOFF << endl;
1158 *fLog << " assume a constant background" << endl;
1159
1160 fConstantBackg = kTRUE;
1161 fDegreeOFF = 0;
1162
1163 TString funcname = "PolyOFF";
1164 Double_t xmin = 0.0;
1165 Double_t xmax = 90.0;
1166
1167 TString formula = "[0]";
1168
1169 fPolyOFF = new TF1(funcname, formula, xmin, xmax);
1170 TList *funclist = fHistOFF->GetListOfFunctions();
1171 funclist->Add(fPolyOFF);
1172
1173 //--------------------
1174 Int_t nparfree = 1;
1175 fChisqOFF = 0.0;
1176 fNdfOFF = fMbinsOFF - nparfree;
1177 fProbOFF = 0.0;
1178 fIstatOFF = 0;
1179
1180 fValuesOFF.Set(1);
1181 fErrorsOFF.Set(1);
1182
1183 Double_t val, err;
1184 val = mean;
1185 err = rms;
1186
1187
1188 fPolyOFF->SetParameter(0, val);
1189 fPolyOFF->SetParError (0, err);
1190
1191 fValuesOFF[0] = val;
1192 fErrorsOFF[0] = err;
1193
1194 fEmaOFF[0][0] = err*err;
1195 fCorrOFF[0][0] = 1.0;
1196 //--------------------
1197 //--------------------------------------------------
1198 // reset the errors of the points in the histogram
1199 for (Int_t i=1; i<=nbins; i++)
1200 {
1201 fHistOFF->SetBinError(i, saveError[i-1]);
1202 }
1203
1204
1205 return kTRUE;
1206 }
1207
1208
1209 //=========== start loop for reducing the degree ==================
1210 // of the polynomial
1211 while (1)
1212 {
1213 //--------------------------------------------------
1214 // prepare fit of a polynomial : (a0 + a1*x + a2*x**2 + a3*x**3 + ...)
1215
1216
1217 TString funcname = "PolyOFF";
1218 Double_t xmin = 0.0;
1219 Double_t xmax = 90.0;
1220
1221 TString formula = "[0]";
1222 TString bra1 = "+[";
1223 TString bra2 = "]";
1224 TString xpower = "*x";
1225 TString newpower = "*x";
1226 for (Int_t i=1; i<=fDegreeOFF; i++)
1227 {
1228 formula += bra1;
1229 formula += i;
1230 formula += bra2;
1231 formula += xpower;
1232
1233 xpower += newpower;
1234 }
1235
1236 //*fLog << "FitPolynomial : formula = " << formula << endl;
1237
1238 fPolyOFF = new TF1(funcname, formula, xmin, xmax);
1239 TList *funclist = fHistOFF->GetListOfFunctions();
1240 funclist->Add(fPolyOFF);
1241
1242 //------------------------
1243 // attention : the dimensions must agree with those in CallMinuit()
1244 const UInt_t npar = fDegreeOFF+1;
1245
1246 TString parname[npar];
1247 TArrayD vinit(npar);
1248 TArrayD step(npar);
1249 TArrayD limlo(npar);
1250 TArrayD limup(npar);
1251 TArrayI fix(npar);
1252
1253 vinit[0] = mean;
1254 vinit[2] = a2init;
1255
1256 for (UInt_t j=0; j<npar; j++)
1257 {
1258 parname[j] = "p";
1259 parname[j] += j+1;
1260
1261 step[j] = vinit[j] != 0.0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
1262 }
1263
1264 // limit the first coefficient of the polynomial to positive values
1265 // because the background must not be negative
1266 limup[0] = fHistOFF->GetEntries();
1267
1268 // use the subsequernt loop if you want to apply the
1269 // constraint : uneven derivatives (at alpha=0) = zero
1270 for (UInt_t j=1; j<npar; j+=2)
1271 {
1272 vinit[j] = 0;
1273 step[j] = 0;
1274 fix[j] = 1;
1275 }
1276
1277 //*fLog << "FitPolynomial : before CallMinuit()" << endl;
1278
1279 MMinuitInterface inter;
1280 const Bool_t rc = inter.CallMinuit(fcnpolyOFF, parname, vinit, step,
1281 limlo, limup, fix, fHistOFF, "Migrad",
1282 kFALSE);
1283
1284 //*fLog << "FitPolynomial : after CallMinuit()" << endl;
1285
1286 if (rc != 0)
1287 {
1288 // *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit failed"
1289 // << endl;
1290 // return kFALSE;
1291 }
1292
1293 //-------------------
1294 // get status of minimization
1295 Double_t fmin = 0;
1296 Double_t fedm = 0;
1297 Double_t errdef = 0;
1298 Int_t npari = 0;
1299 Int_t nparx = 0;
1300
1301 if (gMinuit)
1302 gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fIstatOFF);
1303
1304 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; fmin, fedm, errdef, npari, nparx, fIstat = "
1305 << fmin << ", " << fedm << ", " << errdef << ", " << npari
1306 << ", " << nparx << ", " << fIstatOFF << endl;
1307
1308
1309 //-------------------
1310 // store the results
1311
1312 Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
1313 fChisqOFF = fmin;
1314 fNdfOFF = fMbinsOFF - nparfree;
1315 fProbOFF = TMath::Prob(fChisqOFF, fNdfOFF);
1316
1317
1318 // get fitted parameter values and errors
1319 fValuesOFF.Set(npar);
1320 fErrorsOFF.Set(npar);
1321
1322 for (Int_t j=0; j<=fDegreeOFF; j++)
1323 {
1324 Double_t val, err;
1325 if (gMinuit)
1326 gMinuit->GetParameter(j, val, err);
1327
1328 fPolyOFF->SetParameter(j, val);
1329 fPolyOFF->SetParError(j, err);
1330
1331 fValuesOFF[j] = val;
1332 fErrorsOFF[j] = err;
1333 }
1334
1335 // if the highest coefficient (j0) of the polynomial
1336 // is consistent with zero reduce the degree of the polynomial
1337
1338 Int_t j0 = 0;
1339 for (Int_t j=fDegreeOFF; j>1; j--)
1340 {
1341 // ignore fixed parameters
1342 if (fErrorsOFF[j] == 0)
1343 continue;
1344
1345 // this is the highest coefficient
1346 j0 = j;
1347 break;
1348 }
1349
1350 if (!fReduceDegree || j0==0 || TMath::Abs(fValuesOFF[j0]) > fErrorsOFF[j0])
1351 break;
1352
1353 // reduce the degree of the polynomial
1354 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; reduce the degree of the polynomial from "
1355 << fDegreeOFF << " to " << (j0-2) << endl;
1356 fDegreeOFF = j0 - 2;
1357
1358 funclist->Remove(fPolyOFF);
1359 //if (fPoly)
1360 delete fPolyOFF;
1361 fPolyOFF = NULL;
1362
1363 // delete the Minuit object in order to have independent starting
1364 // conditions for the next minimization
1365 //if (gMinuit)
1366 delete gMinuit;
1367 gMinuit = NULL;
1368 }
1369 //=========== end of loop for reducing the degree ==================
1370 // of the polynomial
1371
1372
1373 //--------------------------------------------------
1374 // get the error matrix of the fitted parameters
1375
1376 if (fIstatOFF >= 1)
1377 {
1378 // error matrix was calculated
1379 if (gMinuit)
1380 gMinuit->mnemat(&fEmatOFF[0][0], fNdimOFF);
1381
1382 // copy covariance matrix into a matrix which includes also the fixed
1383 // parameters
1384 TString name;
1385 Double_t bnd1, bnd2, val, err;
1386 Int_t jvarbl;
1387 Int_t kvarbl;
1388 for (Int_t j=0; j<=fDegreeOFF; j++)
1389 {
1390 if (gMinuit)
1391 gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
1392
1393 for (Int_t k=0; k<=fDegreeOFF; k++)
1394 {
1395 if (gMinuit)
1396 gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
1397
1398 fEmaOFF[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fEmatOFF[jvarbl-1][kvarbl-1];
1399 }
1400 }
1401 }
1402 else
1403 {
1404 // error matrix was not calculated, construct it
1405 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; error matrix not defined"
1406 << endl;
1407 for (Int_t j=0; j<=fDegreeOFF; j++)
1408 {
1409 for (Int_t k=0; k<=fDegreeOFF; k++)
1410 fEmaOFF[j][k] = 0;
1411
1412 fEmaOFF[j][j] = fErrorsOFF[j]*fErrorsOFF[j];
1413 }
1414 }
1415
1416
1417
1418
1419 //--------------------------------------------------
1420 // calculate correlation matrix
1421 for (Int_t j=0; j<=fDegreeOFF; j++)
1422 {
1423 for (Int_t k=0; k<=fDegreeOFF; k++)
1424 {
1425 const Double_t sq = fEmaOFF[j][j]*fEmaOFF[k][k];
1426 fCorrOFF[j][k] =
1427 sq == 0 ? 0 : (fEmaOFF[j][k] / TMath::Sqrt(fEmaOFF[j][j]*fEmaOFF[k][k]));
1428 }
1429 }
1430
1431 //--------------------------------------------------
1432 // reset the errors of the points in the histogram
1433 for (Int_t i=1; i<=nbins; i++)
1434 {
1435 fHistOFF->SetBinError(i, saveError[i-1]);
1436 }
1437
1438 return kTRUE;
1439
1440
1441}
1442
1443
1444
1445
1446
1447
1448
1449
1450// --------------------------------------------------------------------------
1451//
1452// FitPolynomial
1453//
1454// - create a clone 'fHist' of the |alpha| distribution 'fHistOrig'
1455// - fit a polynomial of degree 'fDegree' to the alpha distribution
1456// 'fHist' in the region alphamin < |alpha| < alphamax
1457//
1458// in pathological cases the histogram is rebinned before fitting
1459// (this is done only if fRebin is kTRUE)
1460//
1461// if the highest coefficient of the polynomial is compatible with zero
1462// the fit is repeated with a polynomial of lower degree
1463// (this is done only if fReduceDegree is kTRUE)
1464//
1465//
1466
1467Bool_t MHFindSignificanceONOFF::FitPolynomial()
1468{
1469 //--------------------------------------------------
1470 // check the histogram :
1471 // - calculate initial values of the parameters
1472 // - check for bins with zero entries
1473 // - set minimum errors
1474 // - save the original errors
1475 // - set errors huge outside the fit range
1476 // (in 'fcnpoly' points with huge errors will be ignored)
1477
1478
1479 Double_t dummy = 1.e20;
1480
1481 Double_t mean;
1482 Double_t rms;
1483 Double_t nclose;
1484 Double_t nfar;
1485 Double_t a2init = 0.0;
1486 TArrayD saveError;
1487
1488 Int_t nbins;
1489 Int_t nrebin = 1;
1490
1491 //---------------- start while loop for rebinning -----------------
1492 while(1)
1493 {
1494
1495 fNzero = 0;
1496 fMbins = 0;
1497 fMlow = 0;
1498 fNbgtot = 0.0;
1499
1500 fAlphami = 10000.0;
1501 fAlphamm = 10000.0;
1502 fAlphama = -10000.0;
1503
1504 mean = 0.0;
1505 rms = 0.0;
1506 nclose = 0.0;
1507 nfar = 0.0;
1508
1509 nbins = fHist->GetNbinsX();
1510 saveError.Set(nbins);
1511
1512 for (Int_t i=1; i<=nbins; i++)
1513 {
1514 saveError[i-1] = fHist->GetBinError(i);
1515
1516 // bin should be completely contained in the fit range
1517 // (fAlphamin, fAlphamax)
1518 Double_t xlo = fHist->GetBinLowEdge(i);
1519 Double_t xup = fHist->GetBinLowEdge(i+1);
1520
1521 if ( xlo >= fAlphamin-fEps && xlo <= fAlphamax+fEps &&
1522 xup >= fAlphamin-fEps && xup <= fAlphamax+fEps )
1523 {
1524 fMbins++;
1525
1526 if ( xlo < fAlphami )
1527 fAlphami = xlo;
1528
1529 if ( xup > fAlphama )
1530 fAlphama = xup;
1531
1532 Double_t content = fHist->GetBinContent(i);
1533 fNbgtot += content;
1534
1535 mean += content;
1536 rms += content*content;
1537
1538 // count events in low-alpha and high-alpha region
1539 if ( xlo >= fAlphammm-fEps && xup >= fAlphammm-fEps)
1540 {
1541 nfar += content;
1542 if ( xlo < fAlphamm )
1543 fAlphamm = xlo;
1544 if ( xup < fAlphamm )
1545 fAlphamm = xup;
1546 }
1547 else
1548 {
1549 nclose += content;
1550 if ( xlo > fAlphamm )
1551 fAlphamm = xlo;
1552 if ( xup > fAlphamm )
1553 fAlphamm = xup;
1554 }
1555
1556 // count bins with zero entry
1557 /*if (content <= 0.0)
1558 fNzero++;*/
1559
1560 // set minimum error
1561 /*if (content < 9.0)
1562 {
1563 fMlow += 1;
1564 fHist->SetBinError(i, 3.0);
1565 }*/
1566
1567 //*fLog << "Take : i, content, error = " << i << ", "
1568 // << fHist->GetBinContent(i) << ", "
1569 // << fHist->GetBinError(i) << endl;
1570
1571 continue;
1572 }
1573 // bin is not completely contained in the fit range : set error huge
1574
1575 fHist->SetBinError(i, dummy);
1576
1577 //*fLog << "Omit : i, content, error = " << i << ", "
1578 // << fHist->GetBinContent(i) << ", " << fHist->GetBinError(i)
1579 // << endl;
1580
1581 }
1582
1583 // mean of entries/bin in the fit range
1584 if (fMbins > 0)
1585 {
1586 mean /= ((Double_t) fMbins);
1587 rms /= ((Double_t) fMbins);
1588 }
1589
1590 rms = sqrt( rms - mean*mean );
1591
1592 // if there are no events in the background region
1593 // there is no reason for rebinning
1594 // and this is the condition for assuming a constant background (= 0)
1595 if (mean <= 0.0)
1596 break;
1597
1598 Double_t helpmi = fAlphami*fAlphami*fAlphami;
1599 Double_t helpmm = fAlphamm*fAlphamm*fAlphamm;
1600 Double_t helpma = fAlphama*fAlphama*fAlphama;
1601 Double_t help = (helpma-helpmm) * (fAlphamm-fAlphami)
1602 - (helpmm-helpmi) * (fAlphama-fAlphamm);
1603 if (help != 0.0)
1604 a2init = ( (fAlphamm-fAlphami)*nfar - (fAlphama-fAlphamm)*nclose )
1605 * 1.5 * fHist->GetBinWidth(1) / help;
1606 else
1607 a2init = 0.0;
1608
1609
1610 //--------------------------------------------
1611 // rebin the histogram
1612 // - if a bin has no entries
1613 // - or if there are too many bins with too few entries
1614 // - or if the new bin width would exceed half the size of the
1615 // signal region
1616
1617 if ( !fRebin ||
1618 ( fNzero <= 0 && (Double_t)fMlow<0.05*(Double_t)fMbins ) ||
1619 (Double_t)(nrebin+1)/(Double_t)nrebin * fHist->GetBinWidth(1)
1620 > fAlphasig/2.0 )
1621 {
1622 //*fLog << "before break" << endl;
1623 break;
1624 }
1625
1626 nrebin += 1;
1627 TString histname = fHist->GetName();
1628 delete fHist;
1629 fHist = NULL;
1630
1631 *fLog << "MHFindSignificanceONOFF::FitPolynomial; rebin the |alpha| plot, grouping "
1632 << nrebin << " bins together" << endl;
1633
1634 // TH1::Rebin doesn't work properly
1635 //fHist = fHistOrig->Rebin(nrebin, "Rebinned");
1636 // use private routine RebinHistogram()
1637 fHist = new TH1F;
1638 fHist->Sumw2();
1639 fHist->SetNameTitle(histname, histname);
1640 fHist->UseCurrentStyle();
1641
1642 // do rebinning such that x0 remains a lower bin edge
1643 Double_t x0 = 0.0;
1644 if ( !RebinHistogram(x0, nrebin) )
1645 {
1646 *fLog << "MHFindSignificanceONOFF::FitPolynomial; RebinHistgram() failed"
1647 << endl;
1648 return kFALSE;
1649 }
1650
1651 fHist->SetXTitle("|alpha| [\\circ]");
1652 fHist->SetYTitle("Counts");
1653
1654 }
1655 //---------------- end of while loop for rebinning -----------------
1656
1657
1658 // if there are still too many bins with too few entries don't fit
1659 // and assume a constant background
1660
1661 fConstantBackg = kFALSE;
1662 if ( fNzero > 0 || (Double_t)fMlow>0.05*(Double_t)fMbins )
1663 {
1664 *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit not possible, fNzero, fMlow, fMbins = "
1665 << fNzero << ", " << fMlow << ", " << fMbins << endl;
1666 *fLog << " assume a constant background" << endl;
1667
1668 fConstantBackg = kTRUE;
1669 fDegree = 0;
1670
1671 TString funcname = "Poly";
1672 Double_t xmin = 0.0;
1673 Double_t xmax = 90.0;
1674
1675 TString formula = "[0]";
1676
1677 fPoly = new TF1(funcname, formula, xmin, xmax);
1678 TList *funclist = fHist->GetListOfFunctions();
1679 funclist->Add(fPoly);
1680
1681 //--------------------
1682 Int_t nparfree = 1;
1683 fChisq = 0.0;
1684 fNdf = fMbins - nparfree;
1685 fProb = 0.0;
1686 fIstat = 0;
1687
1688 fValues.Set(1);
1689 fErrors.Set(1);
1690
1691 Double_t val, err;
1692 val = mean;
1693 err = sqrt( mean / (Double_t)fMbins );
1694
1695 fPoly->SetParameter(0, val);
1696 fPoly->SetParError (0, err);
1697
1698 fValues[0] = val;
1699 fErrors[0] = err;
1700
1701 fEma[0][0] = err*err;
1702 fCorr[0][0] = 1.0;
1703 //--------------------
1704
1705 //--------------------------------------------------
1706 // reset the errors of the points in the histogram
1707 for (Int_t i=1; i<=nbins; i++)
1708 {
1709 fHist->SetBinError(i, saveError[i-1]);
1710 }
1711
1712
1713 return kTRUE;
1714 }
1715
1716
1717 //=========== start loop for reducing the degree ==================
1718 // of the polynomial
1719 while (1)
1720 {
1721 //--------------------------------------------------
1722 // prepare fit of a polynomial : (a0 + a1*x + a2*x**2 + a3*x**3 + ...)
1723
1724 TString funcname = "Poly";
1725 Double_t xmin = 0.0;
1726 Double_t xmax = 90.0;
1727
1728 TString formula = "[0]";
1729 TString bra1 = "+[";
1730 TString bra2 = "]";
1731 TString xpower = "*x";
1732 TString newpower = "*x";
1733 for (Int_t i=1; i<=fDegree; i++)
1734 {
1735 formula += bra1;
1736 formula += i;
1737 formula += bra2;
1738 formula += xpower;
1739
1740 xpower += newpower;
1741 }
1742
1743 //*fLog << "FitPolynomial : formula = " << formula << endl;
1744
1745 fPoly = new TF1(funcname, formula, xmin, xmax);
1746 TList *funclist = fHist->GetListOfFunctions();
1747 funclist->Add(fPoly);
1748
1749 //------------------------
1750 // attention : the dimensions must agree with those in CallMinuit()
1751 const UInt_t npar = fDegree+1;
1752
1753 TString parname[npar];
1754 TArrayD vinit(npar);
1755 TArrayD step(npar);
1756 TArrayD limlo(npar);
1757 TArrayD limup(npar);
1758 TArrayI fix(npar);
1759
1760 vinit[0] = mean;
1761 vinit[2] = a2init;
1762
1763 for (UInt_t j=0; j<npar; j++)
1764 {
1765 parname[j] = "p";
1766 parname[j] += j+1;
1767
1768 step[j] = vinit[j] != 0.0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
1769 }
1770
1771 // limit the first coefficient of the polynomial to positive values
1772 // because the background must not be negative
1773 limup[0] = fHist->GetEntries();
1774
1775 // use the subsequernt loop if you want to apply the
1776 // constraint : uneven derivatives (at alpha=0) = zero
1777 for (UInt_t j=1; j<npar; j+=2)
1778 {
1779 vinit[j] = 0;
1780 step[j] = 0;
1781 fix[j] = 1;
1782 }
1783
1784 //*fLog << "FitPolynomial : before CallMinuit()" << endl;
1785
1786 MMinuitInterface inter;
1787 const Bool_t rc = inter.CallMinuit(fcnpoly, parname, vinit, step,
1788 limlo, limup, fix, fHist, "Migrad",
1789 kFALSE);
1790
1791 //*fLog << "FitPolynomial : after CallMinuit()" << endl;
1792
1793 if (rc != 0)
1794 {
1795 // *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit failed"
1796 // << endl;
1797 // return kFALSE;
1798 }
1799
1800
1801 //-------------------
1802 // get status of minimization
1803 Double_t fmin = 0;
1804 Double_t fedm = 0;
1805 Double_t errdef = 0;
1806 Int_t npari = 0;
1807 Int_t nparx = 0;
1808
1809 if (gMinuit)
1810 gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fIstat);
1811
1812 *fLog << "MHFindSignificanceONOFF::FitPolynomial; fmin, fedm, errdef, npari, nparx, fIstat = "
1813 << fmin << ", " << fedm << ", " << errdef << ", " << npari
1814 << ", " << nparx << ", " << fIstat << endl;
1815
1816
1817 //-------------------
1818 // store the results
1819
1820 Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
1821 fChisq = fmin;
1822 fNdf = fMbins - nparfree;
1823 fProb = TMath::Prob(fChisq, fNdf);
1824
1825
1826 // get fitted parameter values and errors
1827 fValues.Set(npar);
1828 fErrors.Set(npar);
1829
1830 for (Int_t j=0; j<=fDegree; j++)
1831 {
1832 Double_t val, err;
1833 if (gMinuit)
1834 gMinuit->GetParameter(j, val, err);
1835
1836 fPoly->SetParameter(j, val);
1837 fPoly->SetParError(j, err);
1838
1839 fValues[j] = val;
1840 fErrors[j] = err;
1841 }
1842
1843
1844 //--------------------------------------------------
1845 // if the highest coefficient (j0) of the polynomial
1846 // is consistent with zero reduce the degree of the polynomial
1847
1848 Int_t j0 = 0;
1849 for (Int_t j=fDegree; j>1; j--)
1850 {
1851 // ignore fixed parameters
1852 if (fErrors[j] == 0)
1853 continue;
1854
1855 // this is the highest coefficient
1856 j0 = j;
1857 break;
1858 }
1859
1860 if (!fReduceDegree || j0==0 || TMath::Abs(fValues[j0]) > fErrors[j0])
1861 break;
1862
1863 // reduce the degree of the polynomial
1864 *fLog << "MHFindSignificanceONOFF::FitPolynomial; reduce the degree of the polynomial from "
1865 << fDegree << " to " << (j0-2) << endl;
1866 fDegree = j0 - 2;
1867
1868 funclist->Remove(fPoly);
1869 //if (fPoly)
1870 delete fPoly;
1871 fPoly = NULL;
1872
1873 // delete the Minuit object in order to have independent starting
1874 // conditions for the next minimization
1875 //if (gMinuit)
1876 delete gMinuit;
1877 gMinuit = NULL;
1878 }
1879 //=========== end of loop for reducing the degree ==================
1880 // of the polynomial
1881
1882
1883 //--------------------------------------------------
1884 // get the error matrix of the fitted parameters
1885
1886
1887 if (fIstat >= 1)
1888 {
1889 // error matrix was calculated
1890 if (gMinuit)
1891 gMinuit->mnemat(&fEmat[0][0], fNdim);
1892
1893 // copy covariance matrix into a matrix which includes also the fixed
1894 // parameters
1895 TString name;
1896 Double_t bnd1, bnd2, val, err;
1897 Int_t jvarbl;
1898 Int_t kvarbl;
1899 for (Int_t j=0; j<=fDegree; j++)
1900 {
1901 if (gMinuit)
1902 gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
1903
1904 for (Int_t k=0; k<=fDegree; k++)
1905 {
1906 if (gMinuit)
1907 gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
1908
1909 fEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fEmat[jvarbl-1][kvarbl-1];
1910 }
1911 }
1912 }
1913 else
1914 {
1915 // error matrix was not calculated, construct it
1916 *fLog << "MHFindSignificanceONOFF::FitPolynomial; error matrix not defined"
1917 << endl;
1918 for (Int_t j=0; j<=fDegree; j++)
1919 {
1920 for (Int_t k=0; k<=fDegree; k++)
1921 fEma[j][k] = 0;
1922
1923 fEma[j][j] = fErrors[j]*fErrors[j];
1924 }
1925 }
1926
1927
1928 //--------------------------------------------------
1929 // calculate correlation matrix
1930 for (Int_t j=0; j<=fDegree; j++)
1931 for (Int_t k=0; k<=fDegree; k++)
1932 {
1933 const Double_t sq = fEma[j][j]*fEma[k][k];
1934 fCorr[j][k] = sq==0 ? 0 : fEma[j][k] / TMath::Sqrt(fEma[j][j]*fEma[k][k]);
1935 }
1936
1937
1938 //--------------------------------------------------
1939 // reset the errors of the points in the histogram
1940 for (Int_t i=1; i<=nbins; i++)
1941 fHist->SetBinError(i, saveError[i-1]);
1942
1943
1944 return kTRUE;
1945}
1946
1947// --------------------------------------------------------------------------
1948//
1949// ReBinHistogramOFF
1950//
1951// rebin the histogram 'fHistOrig' by grouping 'nrebin' bins together
1952// put the result into the histogram 'fHist'
1953// the rebinning is made such that 'x0' remains a lower bound of a bin
1954//
1955
1956Bool_t MHFindSignificanceONOFF::RebinHistogramOFF(Double_t x0, Int_t nrebin)
1957{
1958 //-----------------------------------------
1959 // search bin i0 which has x0 as lower edge
1960
1961 Int_t i0 = -1;
1962 Int_t nbold = fHistOrigOFF->GetNbinsX();
1963 for (Int_t i=1; i<=nbold; i++)
1964 {
1965 if (TMath::Abs(fHistOrigOFF->GetBinLowEdge(i) - x0) < 1.e-4 )
1966 {
1967 i0 = i;
1968 break;
1969 }
1970 }
1971
1972 if (i0 == -1)
1973 {
1974 i0 = 1;
1975 *fLog << "MHFindsignificanceOFF::RebinOFF; no bin found with " << x0
1976 << " as lower edge, start rebinning with bin 1" << endl;
1977 }
1978
1979 Int_t istart = i0 - nrebin * ( (i0-1)/nrebin );
1980
1981 //-----------------------------------------
1982 // get new bin edges
1983
1984 const Int_t nbnew = (nbold-istart+1) / nrebin;
1985 const Double_t xmin = fHistOrigOFF->GetBinLowEdge(istart);
1986 const Double_t xmax = xmin + (Double_t)nbnew * nrebin * fHistOrigOFF->GetBinWidth(1);
1987 fHistOFF->SetBins(nbnew, xmin, xmax);
1988
1989 *fLog << "MHFindSignificanceONOFF::ReBinOFF; x0, i0, nbold, nbnew, xmin, xmax = "
1990 << x0 << ", " << i0 << ", " << nbold << ", " << nbnew << ", "
1991 << xmin << ", " << xmax << endl;
1992
1993 //-----------------------------------------
1994 // get new bin entries
1995
1996 for (Int_t i=1; i<=nbnew; i++)
1997 {
1998 Int_t j = nrebin*(i-1) + istart;
1999
2000 Double_t content = 0;
2001 Double_t error2 = 0;
2002 for (Int_t k=0; k<nrebin; k++)
2003 {
2004 content += fHistOrigOFF->GetBinContent(j+k);
2005 error2 += fHistOrigOFF->GetBinError(j+k) * fHistOrigOFF->GetBinError(j+k);
2006 }
2007 fHistOFF->SetBinContent(i, content);
2008 fHistOFF->SetBinError (i, sqrt(error2));
2009 }
2010 fHistOFF->SetEntries( fHistOrigOFF->GetEntries() );
2011
2012 return kTRUE;
2013}
2014
2015
2016
2017
2018// --------------------------------------------------------------------------
2019
2020
2021
2022//
2023// ReBinHistogram
2024//
2025// rebin the histogram 'fHistOrig' by grouping 'nrebin' bins together
2026// put the result into the histogram 'fHist'
2027// the rebinning is made such that 'x0' remains a lower bound of a bin
2028//
2029
2030Bool_t MHFindSignificanceONOFF::RebinHistogram(Double_t x0, Int_t nrebin)
2031{
2032 //-----------------------------------------
2033 // search bin i0 which has x0 as lower edge
2034
2035 Int_t i0 = -1;
2036 Int_t nbold = fHistOrig->GetNbinsX();
2037 for (Int_t i=1; i<=nbold; i++)
2038 {
2039 if (TMath::Abs(fHistOrig->GetBinLowEdge(i) - x0) < 1.e-4 )
2040 {
2041 i0 = i;
2042 break;
2043 }
2044 }
2045
2046 if (i0 == -1)
2047 {
2048 i0 = 1;
2049 *fLog << "MHFindsignificance::Rebin; no bin found with " << x0
2050 << " as lower edge, start rebinning with bin 1" << endl;
2051 }
2052
2053 Int_t istart = i0 - nrebin * ( (i0-1)/nrebin );
2054
2055 //-----------------------------------------
2056 // get new bin edges
2057
2058 const Int_t nbnew = (nbold-istart+1) / nrebin;
2059 const Double_t xmin = fHistOrig->GetBinLowEdge(istart);
2060 const Double_t xmax = xmin + (Double_t)nbnew * nrebin * fHistOrig->GetBinWidth(1);
2061 fHist->SetBins(nbnew, xmin, xmax);
2062
2063 *fLog << "MHFindSignificanceONOFF::ReBin; x0, i0, nbold, nbnew, xmin, xmax = "
2064 << x0 << ", " << i0 << ", " << nbold << ", " << nbnew << ", "
2065 << xmin << ", " << xmax << endl;
2066
2067 //-----------------------------------------
2068 // get new bin entries
2069
2070 for (Int_t i=1; i<=nbnew; i++)
2071 {
2072 Int_t j = nrebin*(i-1) + istart;
2073
2074 Double_t content = 0;
2075 Double_t error2 = 0;
2076 for (Int_t k=0; k<nrebin; k++)
2077 {
2078 content += fHistOrig->GetBinContent(j+k);
2079 error2 += fHistOrig->GetBinError(j+k) * fHistOrig->GetBinError(j+k);
2080 }
2081 fHist->SetBinContent(i, content);
2082 fHist->SetBinError (i, sqrt(error2));
2083 }
2084 fHist->SetEntries( fHistOrig->GetEntries() );
2085
2086 return kTRUE;
2087}
2088
2089// --------------------------------------------------------------------------
2090//
2091// FitGaussPoly
2092//
2093// fits a (Gauss + polynomial function) to the alpha distribution 'fhist'
2094//
2095//
2096Bool_t MHFindSignificanceONOFF::FitGaussPoly()
2097{
2098 *fLog << "Entry FitGaussPoly" << endl;
2099
2100 //--------------------------------------------------
2101 // check the histogram :
2102 // - calculate initial values of the parameters
2103 // - check for bins with zero entries
2104 // - set minimum errors
2105 // - save the original errors
2106 // - set errors huge outside the fit range
2107 // (in 'fcnpoly' points with huge errors will be ignored)
2108
2109
2110 Double_t dummy = 1.e20;
2111
2112 fGNzero = 0;
2113 fGMbins = 0;
2114
2115 //------------------------------------------
2116 // if a constant background has been assumed (due to low statistics)
2117 // fit only in the signal region
2118 if ( !fConstantBackg )
2119 {
2120 fAlphalow = 0.0;
2121 fAlphahig = fAlphamax;
2122 }
2123 else
2124 {
2125 fAlphalow = 0.0;
2126 fAlphahig = 2.0*fAlphasig>25.0 ? 25.0 : 2.0*fAlphasig;
2127 }
2128 //------------------------------------------
2129
2130
2131 fAlphalo = 10000.0;
2132 fAlphahi = -10000.0;
2133
2134
2135 Int_t nbins = fHist->GetNbinsX();
2136 TArrayD saveError(nbins);
2137
2138 for (Int_t i=1; i<=nbins; i++)
2139 {
2140 saveError[i-1] = fHist->GetBinError(i);
2141
2142 // bin should be completely contained in the fit range
2143 // (fAlphalow, fAlphahig)
2144 Double_t xlo = fHist->GetBinLowEdge(i);
2145 Double_t xup = fHist->GetBinLowEdge(i+1);
2146
2147 if ( xlo >= fAlphalow-fEps && xlo <= fAlphahig+fEps &&
2148 xup >= fAlphalow-fEps && xup <= fAlphahig+fEps )
2149 {
2150 fGMbins++;
2151
2152 if ( xlo < fAlphalo )
2153 fAlphalo = xlo;
2154
2155 if ( xup > fAlphahi )
2156 fAlphahi = xup;
2157
2158 Double_t content = fHist->GetBinContent(i);
2159
2160
2161 // count bins with zero entry
2162 if (content <= 0.0)
2163 fGNzero++;
2164
2165 // set minimum error
2166 if (content < 9.0)
2167 fHist->SetBinError(i, 3.0);
2168
2169 //*fLog << "Take : i, content, error = " << i << ", "
2170 // << fHist->GetBinContent(i) << ", "
2171 // << fHist->GetBinError(i) << endl;
2172
2173 continue;
2174 }
2175 // bin is not completely contained in the fit range : set error huge
2176
2177 fHist->SetBinError(i, dummy);
2178
2179 //*fLog << "Omit : i, content, error = " << i << ", "
2180 // << fHist->GetBinContent(i) << ", " << fHist->GetBinError(i)
2181 // << endl;
2182
2183 }
2184
2185
2186 // if a bin has no entries don't fit
2187 if (fGNzero > 0)
2188 {
2189 *fLog << "MHFindSignificanceONOFF::FitGaussPoly; out of " << fGMbins
2190 << " bins there are " << fGNzero
2191 << " bins with zero entry" << endl;
2192
2193 fGPoly = NULL;
2194 return kFALSE;
2195 }
2196
2197
2198 //--------------------------------------------------
2199 // prepare fit of a (polynomial+Gauss) :
2200 // (a0 + a1*x + a2*x**2 + a3*x**3 + ...) + A*exp( -0.5*((x-x0)/sigma)**2 )
2201
2202 TString funcname = "PolyGauss";
2203 Double_t xmin = 0.0;
2204 Double_t xmax = 90.0;
2205
2206 TString xpower = "*x";
2207 TString newpower = "*x";
2208
2209 TString formulaBackg = "[0]";
2210 for (Int_t i=1; i<=fDegree; i++)
2211 formulaBackg += Form("+[%d]*x^%d", i, i);
2212
2213 const TString formulaGauss =
2214 Form("[%d]/[%d]*exp(-0.5*((x-[%d])/[%d])^2)",
2215 fDegree+1, fDegree+3, fDegree+2, fDegree+3);
2216
2217 TString formula = formulaBackg;
2218 formula += "+";
2219 formula += formulaGauss;
2220
2221 *fLog << "FitGaussPoly : formulaBackg = " << formulaBackg << endl;
2222 *fLog << "FitGaussPoly : formulaGauss = " << formulaGauss << endl;
2223 *fLog << "FitGaussPoly : formula = " << formula << endl;
2224
2225 fGPoly = new TF1(funcname, formula, xmin, xmax);
2226 TList *funclist = fHist->GetListOfFunctions();
2227 funclist->Add(fGPoly);
2228
2229 fGBackg = new TF1("Backg", formulaBackg, xmin, xmax);
2230 //funclist->Add(fGBackg);
2231
2232 //------------------------
2233 // attention : the dimensions must agree with those in CallMinuit()
2234 Int_t npar = fDegree+1 + 3;
2235
2236 TString parname[npar];
2237 TArrayD vinit(npar);
2238 TArrayD step(npar);
2239 TArrayD limlo(npar);
2240 TArrayD limup(npar);
2241 TArrayI fix(npar);
2242
2243
2244 // take as initial values for the polynomial
2245 // the result from the polynomial fit
2246 for (Int_t j=0; j<=fDegree; j++)
2247 vinit[j] = fPoly->GetParameter(j);
2248
2249 Double_t sigma = 8;
2250 vinit[fDegree+1] = 2.0 * fNexONOFF * fHist->GetBinWidth(1) / TMath::Sqrt(TMath::Pi()*2);
2251 vinit[fDegree+2] = 0;
2252 vinit[fDegree+3] = sigma;
2253
2254 *fLog << "FitGaussPoly : starting value for Gauss-amplitude = "
2255 << vinit[fDegree+1] << endl;
2256
2257 for (Int_t j=0; j<npar; j++)
2258 {
2259 parname[j] = "p";
2260 parname[j] += j+1;
2261
2262 step[j] = vinit[j]!=0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
2263 }
2264
2265 // limit the first coefficient of the polynomial to positive values
2266 // because the background must not be negative
2267 limup[0] = fHist->GetEntries()*10;
2268
2269 // limit the sigma of the Gauss function
2270 limup[fDegree+3] = 20;
2271
2272
2273 // use the subsequernt loop if you want to apply the
2274 // constraint : uneven derivatives (at alpha=0) = zero
2275 for (Int_t j=1; j<=fDegree; j+=2)
2276 {
2277 vinit[j] = 0;
2278 step[j] = 0;
2279 fix[j] = 1;
2280 }
2281
2282 // fix position of Gauss function
2283 vinit[fDegree+2] = 0;
2284 step[fDegree+2] = 0;
2285 fix[fDegree+2] = 1;
2286
2287 // if a constant background has been assumed (due to low statistics)
2288 // fix the background
2289 if (fConstantBackg)
2290 {
2291 step[0] = 0;
2292 fix[0] = 1;
2293 }
2294
2295 MMinuitInterface inter;
2296 const Bool_t rc = inter.CallMinuit(fcnpolygauss, parname, vinit, step,
2297 limlo, limup, fix, fHist, "Migrad",
2298 kFALSE);
2299
2300 if (rc != 0)
2301 {
2302 // *fLog << "MHFindSignificanceONOFF::FitGaussPoly; (polynomial+Gauss) fit failed"
2303 // << endl;
2304 // return kFALSE;
2305 }
2306
2307
2308 //-------------------
2309 // get status of the minimization
2310 Double_t fmin;
2311 Double_t fedm;
2312 Double_t errdef;
2313 Int_t npari;
2314 Int_t nparx;
2315
2316 if (gMinuit)
2317 gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fGIstat);
2318
2319 *fLog << "MHFindSignificanceONOFF::FitGaussPoly; fmin, fedm, errdef, npari, nparx, fGIstat = "
2320 << fmin << ", " << fedm << ", " << errdef << ", " << npari
2321 << ", " << nparx << ", " << fGIstat << endl;
2322
2323
2324 //-------------------
2325 // store the results
2326
2327 Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
2328 fGChisq = fmin;
2329 fGNdf = fGMbins - nparfree;
2330 fGProb = TMath::Prob(fGChisq, fGNdf);
2331
2332
2333 // get fitted parameter values and errors
2334 fGValues.Set(npar);
2335 fGErrors.Set(npar);
2336
2337 for (Int_t j=0; j<npar; j++)
2338 {
2339 Double_t val, err;
2340 if (gMinuit)
2341 gMinuit->GetParameter(j, val, err);
2342
2343 fGPoly->SetParameter(j, val);
2344 fGPoly->SetParError(j, err);
2345
2346 fGValues[j] = val;
2347 fGErrors[j] = err;
2348
2349 if (j <=fDegree)
2350 {
2351 fGBackg->SetParameter(j, val);
2352 fGBackg->SetParError(j, err);
2353 }
2354 }
2355
2356 fSigmaGauss = fGValues[fDegree+3];
2357 fdSigmaGauss = fGErrors[fDegree+3];
2358 // fitted total number of excess events
2359 fNexGauss = fGValues[fDegree+1] * TMath::Sqrt(TMath::Pi()*2) /
2360 (fHist->GetBinWidth(1)*2 );
2361 fdNexGauss = fNexGauss * fGErrors[fDegree+1]/fGValues[fDegree+1];
2362
2363 //--------------------------------------------------
2364 // get the error matrix of the fitted parameters
2365
2366
2367 if (fGIstat >= 1)
2368 {
2369 // error matrix was calculated
2370 if (gMinuit)
2371 gMinuit->mnemat(&fGEmat[0][0], fGNdim);
2372
2373 // copy covariance matrix into a matrix which includes also the fixed
2374 // parameters
2375 TString name;
2376 Double_t bnd1, bnd2, val, err;
2377 Int_t jvarbl;
2378 Int_t kvarbl;
2379 for (Int_t j=0; j<npar; j++)
2380 {
2381 if (gMinuit)
2382 gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
2383
2384 for (Int_t k=0; k<npar; k++)
2385 {
2386 if (gMinuit)
2387 gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
2388
2389 fGEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fGEmat[jvarbl-1][kvarbl-1];
2390 }
2391 }
2392 }
2393 else
2394 {
2395 // error matrix was not calculated, construct it
2396 *fLog << "MHFindSignificanceONOFF::FitPolynomial; error matrix not defined"
2397 << endl;
2398 for (Int_t j=0; j<npar; j++)
2399 {
2400 for (Int_t k=0; k<npar; k++)
2401 fGEma[j][k] = 0;
2402
2403 fGEma[j][j] = fGErrors[j]*fGErrors[j];
2404 }
2405 }
2406
2407
2408 //--------------------------------------------------
2409 // calculate correlation matrix
2410 for (Int_t j=0; j<npar; j++)
2411 {
2412 for (Int_t k=0; k<npar; k++)
2413 {
2414 const Double_t sq = fGEma[j][j]*fGEma[k][k];
2415 fGCorr[j][k] = sq==0 ? 0 : fGEma[j][k] / sqrt( fGEma[j][j]*fGEma[k][k] );
2416 }
2417 }
2418
2419
2420 //--------------------------------------------------
2421 // reset the errors of the points in the histogram
2422 for (Int_t i=1; i<=nbins; i++)
2423 fHist->SetBinError(i, saveError[i-1]);
2424
2425 return kTRUE;
2426
2427}
2428
2429
2430
2431// --------------------------------------------------------------------------
2432//
2433// DetExcessONOFF
2434//
2435// using the result of the polynomial fit (fValuesOFF), DetExcessONOFF determines
2436//
2437// - the total number of events in the signal region (fNon)
2438// - the number of backgound events (fitted) in the signal region (fNoffSigFitted)
2439// - the number of OFF (normalized evetns) in the alpha OFF distribution (fNoffSig)
2440// - the number of excess events (fNexONOFF)
2441// - the number of excess events using the fitted OFF (fNexONOFFFitted)
2442// - the effective number of background events (fNoff), and fGamma :
2443// fNoffSigFitted = fGamma * fNoff; fdNoffSigFitted = fGamma * sqrt(fNoff);
2444//
2445// It assumed that the polynomial is defined as
2446// a0 + a1*x + a2*x**2 + a3*x**3 + ..
2447//
2448// and that the alpha distribution has the range 0 < alpha < 90 degrees
2449//
2450
2451Bool_t MHFindSignificanceONOFF::DetExcessONOFF()
2452{
2453 //*fLog << "MHFindSignificanceONOFF::DetExcessONOFF;" << endl;
2454
2455//--------------------------------------------
2456 // calculate the total number of events (fNon) in the signal region
2457
2458 fNon = 0.0;
2459 fdNon = 0.0;
2460
2461 Double_t alphaup = -1000.0;
2462 Double_t binwidth = fHist->GetBinWidth(1);
2463
2464 Int_t nbins = fHist->GetNbinsX();
2465 for (Int_t i=1; i<=nbins; i++)
2466 {
2467 Double_t xlo = fHist->GetBinLowEdge(i);
2468 Double_t xup = fHist->GetBinLowEdge(i+1);
2469
2470 // bin must be completely contained in the signal region
2471 if ( xlo <= (fAlphasig+fEps) && xup <= (fAlphasig+fEps) )
2472 {
2473 Double_t width = fabs(xup-xlo);
2474 if (fabs(width-binwidth) > fEps)
2475 {
2476 *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alpha plot has variable binning, which is not allowed"
2477 << endl;
2478 return kFALSE;
2479 }
2480
2481 if (xup > alphaup)
2482 alphaup = xup;
2483
2484 fNon += fHist->GetBinContent(i);
2485 fdNon += fHist->GetBinError(i) * fHist->GetBinError(i);
2486 }
2487 }
2488 fdNon = sqrt(fdNon);
2489
2490 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF:" << endl
2491 << "fNon = " << fNon << " ; fdNon = " << fdNon << endl;
2492
2493
2494 // tmp
2495 //if (fNon == 0)
2496 //{
2497 // cout << "ERROR... WITH COUNTED OF ON EVETNS: fAlphasig, fEps = "
2498// << fAlphasig << ", " << fEps << endl;
2499 // fHist-> DrawCopy();
2500 // gPad -> SaveAs("HistON.ps");
2501 //}
2502 // endtmp
2503
2504
2505 // the actual signal range is :
2506 if (alphaup == -1000.0)
2507 return kFALSE;
2508
2509 fAlphasi = alphaup;
2510
2511
2512 //*fLog << "fAlphasi, fNon, fdNon, binwidth, fDegree = " << fAlphasi << ", "
2513 // << fNon << ", " << fdNon << ", " << binwidth << ", "
2514 // << fDegree << endl;
2515
2516
2517 // Calculate the number of OFF events in the signal region
2518 // ___________________________________________________
2519
2520
2521 fNoffSig = 0.0;
2522 fdNoffSig = 0.0;
2523
2524 Double_t alphaupOFF = -1000.0;
2525 Double_t binwidthOFF = fHistOFF->GetBinWidth(1);
2526
2527 Int_t nbinsOFF = fHistOFF->GetNbinsX();
2528
2529 for (Int_t i=1; i<=nbinsOFF; i++)
2530 {
2531 Double_t xlo = fHistOFF->GetBinLowEdge(i);
2532 Double_t xup = fHistOFF->GetBinLowEdge(i+1);
2533
2534 // bin must be completely contained in the signal region
2535 if ( xlo <= (fAlphasig+fEps) && xup <= (fAlphasig+fEps) )
2536 {
2537 Double_t width = fabs(xup-xlo);
2538 if (fabs(width-binwidthOFF) > fEps)
2539 {
2540 *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alphaOFF plot has variable binning, which is not allowed"
2541 << endl;
2542 return kFALSE;
2543 }
2544
2545 if (xup > alphaupOFF)
2546 alphaup = xup;
2547
2548 fNoffSig += fHistOFF->GetBinContent(i);
2549 fdNoffSig += fHistOFF->GetBinError(i) * fHistOFF->GetBinError(i);
2550 }
2551 }
2552 fdNoffSig = sqrt(fdNoffSig);
2553
2554 // tmp
2555 //if (fNoffSig == 0)
2556 //{
2557 // cout << "ERROR... WITH COUNTED OF OFF EVETNS: fAlphasig, fEps = "
2558// << fAlphasig << ", " << fEps << endl;
2559 // fHistOFF-> DrawCopy();
2560 // gPad -> SaveAs("HistOFF.ps");
2561 //}
2562 //endtmp
2563
2564 // the actual signal range is :
2565 if (alphaup == -1000.0)
2566 return kFALSE;
2567
2568 fAlphasiOFF = alphaup;
2569
2570 if (fabs(fAlphasiOFF - fAlphasi) > fEps)
2571 {
2572 *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; fAlphasiOFF ("
2573 << fAlphasiOFF << ") is not equal to fAlphasi ("
2574 << fAlphasi << "), this is something that should not happen"
2575 << endl;
2576
2577 //return kFALSE; It might happen in pathological cases (very few OFF)
2578 // and I want to see the results, anyhow
2579 }
2580
2581
2582
2583 // Calculate the number of OFF events in the total OFF region
2584 // defined by fAlphaminOFF and fAlphamaxOFF
2585 // ___________________________________________________
2586
2587 fNoffTot = 0.0;
2588 fdNoffTot = 0.0;
2589
2590
2591
2592 for (Int_t i=1; i<=nbinsOFF; i++)
2593 {
2594 Double_t xlo = fHistOFF->GetBinLowEdge(i);
2595 Double_t xup = fHistOFF->GetBinLowEdge(i+1);
2596
2597 // bin must be completely contained in the signal region
2598 if ( xlo >= (fAlphaminOFF-fEps) && xup <= (fAlphamaxOFF+fEps) )
2599 {
2600 Double_t width = fabs(xup-xlo);
2601 if (fabs(width-binwidthOFF) > fEps)
2602 {
2603 *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alphaOFF plot has variable binning, which is not allowed"
2604 << endl;
2605 return kFALSE;
2606 }
2607
2608 fNoffTot += fHistOFF->GetBinContent(i);
2609 fdNoffTot += fHistOFF->GetBinError(i) * fHistOFF->GetBinError(i);
2610 }
2611 }
2612 fdNoffTot = sqrt(fdNoffTot);
2613
2614
2615
2616
2617
2618
2619
2620
2621
2622 //--------------------------------------------
2623 // calculate the number of OFF fitted events (fNoffSigFitted) in the signal region
2624 // and its error (fdNoffSigFitted)
2625
2626
2627
2628 if (fUseFittedQuantities)
2629 {
2630 //--------------------------------------------
2631 // calculate the number of OFF fitted events (fNoffSigFitted) in the signal region
2632 // and its error (fdNoffSigFitted)
2633
2634 Double_t fac = 1.0/binwidthOFF;
2635
2636 fNoffSigFitted = 0.0;
2637 Double_t altothejplus1 = fAlphasi; // Limit for signal found for ON data is used.
2638 for (Int_t j=0; j<=fDegreeOFF; j++)
2639 {
2640 fNoffSigFitted += fValuesOFF[j] * altothejplus1 / ((Double_t)(j+1));
2641 altothejplus1 *= fAlphasi;
2642 }
2643 fNoffSigFitted *= fac;
2644
2645 // derivative of fNoffSigFitted
2646 Double_t facj;
2647 Double_t fack;
2648
2649 Double_t sum = 0.0;
2650 altothejplus1 = fAlphasi;
2651 for (Int_t j=0; j<=fDegreeOFF; j++)
2652 {
2653 facj = altothejplus1 / ((Double_t)(j+1));
2654
2655 Double_t altothekplus1 = fAlphasi;
2656 for (Int_t k=0; k<=fDegreeOFF; k++)
2657 {
2658 fack = altothekplus1 / ((Double_t)(k+1));
2659 sum += facj * fack * fEmaOFF[j][k];
2660 altothekplus1 *= fAlphasi;
2661 }
2662 altothejplus1 *= fAlphasi;
2663 }
2664 sum *= fac*fac;
2665
2666 if (sum < 0.0)
2667 {
2668 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; error squared is negative"
2669 << endl;
2670 return kFALSE;
2671 }
2672
2673 fdNoffSigFitted = sqrt(sum);
2674
2675
2676 // We can now compare fNoffSig with fNoffSigFitted (and their errors)
2677 // NUmbers should agree within 10 % (errors within 20%)
2678
2679 if (fabs(fNoffSig - fNoffSigFitted) > 0.1 * fNoffSigFitted)
2680 {
2681 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; number of OFF events and Fitted number of OFF events in signal region do not agree (within 10 %)" << endl;
2682
2683 *fLog << "fNoffSig = " << fNoffSig << " ; fNoffSigFitted = " << fNoffSigFitted << endl;
2684
2685
2686 // return kFALSE; NOt yet...
2687 }
2688
2689 /*
2690 if (fabs(fdNoffSig - fdNoffSigFitted) > 0.2 * fdNoffSigFitted)
2691 {
2692 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; error in number of OFF events and error in Fitted number of OFF events in signal region do not agree (within 20 %)"
2693 << endl;
2694
2695
2696 *fLog << "fdNoffSig = " << fdNoffSig << " ; fdNoffSigFitted = " << fdNoffSigFitted << endl;
2697
2698
2699 //return kFALSE; NOt yet...
2700 }
2701
2702 */
2703
2704
2705
2706 // Calculate the number of OFF events in the whole fit region (fAlphaminOFF-fAlphamaxOFF)
2707 // ___________________________________________________
2708
2709
2710
2711 fNoffTotFitted = 0.0;
2712
2713 altothejplus1 = fAlphamaxOFF; // Limit for OFF data fit
2714 for (Int_t j=0; j<=fDegreeOFF; j++)
2715 {
2716 fNoffTotFitted += fValuesOFF[j] * altothejplus1 / ((Double_t)(j+1));
2717 altothejplus1 *= fAlphamaxOFF;
2718 }
2719 fNoffTotFitted *= fac;
2720
2721 // derivative of fdNoffTotFitted
2722
2723
2724 sum = 0.0;
2725 altothejplus1 = fAlphamaxOFF;
2726 for (Int_t j=0; j<=fDegreeOFF; j++)
2727 {
2728 facj = altothejplus1 / ((Double_t)(j+1));
2729
2730 Double_t altothekplus1 = fAlphamaxOFF;
2731 for (Int_t k=0; k<=fDegreeOFF; k++)
2732 {
2733 fack = altothekplus1 / ((Double_t)(k+1));
2734
2735 sum += facj * fack * fEmaOFF[j][k];
2736 altothekplus1 *= fAlphamaxOFF;
2737 }
2738 altothejplus1 *= fAlphamaxOFF;
2739 }
2740 sum *= fac*fac;
2741
2742 if (sum < 0.0)
2743 {
2744 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; error squared is negative"
2745 << endl;
2746 return kFALSE;
2747 }
2748
2749 fdNoffTotFitted = sqrt(sum);
2750
2751
2752
2753 }
2754
2755else
2756 {
2757 fNoffSigFitted = 0.0;
2758 fdNoffSigFitted = 0.0;
2759 fNoffTotFitted = 0.0;
2760 fdNoffTotFitted = 0.0;
2761 }
2762
2763
2764
2765 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF; INFO ABOOUT COMPUTED OFF EVENTS..." << endl
2766 << "fNoffSig = " << fNoffSig << "; fdNoffSig = " << fdNoffSig << endl
2767 << "fNoffSigFitted = " << fNoffSigFitted
2768 << "; fdNoffSigFitted = " << fdNoffSigFitted << endl;
2769
2770
2771
2772
2773
2774 //--------------------------------------------
2775 // calculate the number of excess events in the signal region
2776
2777 fNexONOFF = fNon - fNoffSig*fNormFactor;
2778 fNexONOFFFitted = fNon - fNoffSigFitted*fNormFactor;
2779
2780
2781 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF;" << endl
2782 << "fNexONOFF (= fNon - fNoffSig*fNormFactor) = " << fNexONOFF << endl
2783 << "fNexONOFFFitted (= fNon - fNoffSigFitted*fNormFactor) = " << fNexONOFFFitted << endl;
2784
2785
2786 //--------------------------------------------
2787 // calculate the effective number of background events (fNoff) , and fGamma :
2788 // fNbg = fGamma * fNoff = fNormFactor* fNoffSigFitted;
2789 // dfNbg = fGamma * sqrt(fNoff) = fNormFactor * fdNoffSigFitted;
2790
2791 if (fNoffSigFitted < 0.0)
2792 {
2793 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF; number of fitted OFF events in signal region is negative, fNoffSigFitted, fdNoffSigFitted = "
2794 << fNoffSigFitted << ", " << fdNoffSigFitted << endl;
2795
2796 fGamma = 1.0;
2797 fNoff = 0.0;
2798 return kFALSE;
2799 }
2800
2801 if (fNoffSigFitted > 0.0)
2802 {
2803 fGamma = fNormFactor * fdNoffSigFitted*fdNoffSigFitted/fNoffSigFitted;
2804 fNoff = fNormFactor * fNoffSigFitted/fGamma;
2805 }
2806 else
2807 {
2808 fGamma = 1.0;
2809 fNoff = 0.0;
2810 }
2811
2812 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF: " << endl
2813 << "fGamma = " << fGamma << " ; fNoff = " << fNoff << endl;
2814
2815
2816 return kTRUE;
2817}
2818
2819
2820
2821
2822
2823// --------------------------------------------------------------------------
2824//
2825// SigmaLiMa
2826//
2827// calculates the significance according to Li & Ma
2828// ApJ 272 (1983) 317; formula 17
2829//
2830Bool_t MHFindSignificanceONOFF::SigmaLiMa(Double_t non, Double_t noff,
2831 Double_t gamma, Double_t *siglima)
2832{
2833 if (gamma <= 0.0 || non <= 0.0 || noff <= 0.0)
2834 {
2835 *siglima = 0.0;
2836 return kFALSE;
2837 }
2838
2839 Double_t help1 = non * log( (1.0+gamma)*non / (gamma*(non+noff)) );
2840 Double_t help2 = noff * log( (1.0+gamma)*noff / ( non+noff ) );
2841 *siglima = sqrt( 2.0 * (help1+help2) );
2842
2843 Double_t nex = non - gamma*noff;
2844 if (nex < 0.0)
2845 *siglima = - *siglima;
2846
2847 //*fLog << "MHFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = "
2848 // << non << ", " << noff << ", " << gamma << ", " << *siglima << endl;
2849
2850 return kTRUE;
2851}
2852
2853
2854
2855// calculates the significance according to Li & Ma
2856// ApJ 272 (1983) 317; formula 5
2857//
2858Bool_t MHFindSignificanceONOFF::SigmaLiMaForm5(Double_t non, Double_t noff,
2859 Double_t gamma, Double_t *siglima)
2860{
2861 if (gamma <= 0.0 || non <= 0.0 || noff <= 0.0)
2862 {
2863 *siglima = 0.0;
2864 return kFALSE;
2865 }
2866
2867 Double_t nex = non - (gamma*noff);
2868 Double_t tmp = non + (gamma*gamma)*noff;
2869 tmp = TMath::Sqrt(tmp);
2870
2871 *siglima = nex/tmp;
2872
2873 if (nex < 0.0)
2874 *siglima = - *siglima;
2875
2876 //*fLog << "MHFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = "
2877 // << non << ", " << noff << ", " << gamma << ", " << *siglima << endl;
2878
2879 return kTRUE;
2880}
2881
2882
2883
2884
2885// --------------------------------------------------------------------------
2886//
2887
2888// Following function computes a clone of fHistOFF and normalizes
2889// contents, errors and fPolyOFF (if exists) with the fNormFactor.
2890// This normalized OFF hist will be used when plotting OFF data
2891// together with ON data.
2892
2893Bool_t MHFindSignificanceONOFF::ComputeHistOFFNormalized()
2894{
2895
2896
2897 if (!fHist)
2898 {
2899 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fHist does not exist, normalization of HistOFF can not be performed properly..."
2900 << endl;
2901 return kFALSE;
2902
2903 }
2904
2905
2906
2907 if (!fHistOFF)
2908 {
2909 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fHistOFF does not exist, hence can not be normalized"
2910 << endl;
2911 return kFALSE;
2912
2913 }
2914
2915
2916 if (fNormFactor <= 0)
2917 {
2918 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fNormFactor is ZERO or NEGATIVE, it might have not been defined yet..."
2919 << endl;
2920 return kFALSE;
2921
2922 }
2923
2924
2925 Double_t BinWidthAlphaON = fHist -> GetBinWidth(1);
2926 Double_t BinWidthAlphaOFF = fHistOFF -> GetBinWidth(1);
2927 Double_t BinWidthRatioONOFF = BinWidthAlphaON/BinWidthAlphaOFF;
2928
2929
2930 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; INFO about alpha ON, OFF histo bins"
2931 << endl
2932 // << "fHist bin width = " << BinWidthAlphaON
2933 << " fHistOFF bin width = " << BinWidthAlphaOFF << endl;
2934
2935
2936
2937 TString fHistOFFNormalizedName = fHistOFF -> GetName();
2938 fHistOFFNormalizedName += (" (Normalized)");
2939 // fHistOFFNormalized = (TH1*) fHistOFF -> Clone();
2940 // fHistOFFNormalized -> SetNameTitle(fHistOFFNormalizedName, fHistOFFNormalizedName);
2941
2942
2943 Int_t nbinsOFFNormalized = 0;
2944 Int_t nbinsOFF = 0;
2945 Double_t xlow = 0.0;
2946 Double_t xup = 0.0;
2947 Double_t content = 0.0;
2948 Double_t error = 0.0;
2949 Double_t BinCenter = 0.0;
2950
2951
2952 // Bins for normalized OFF histo will be the ones of ON histo
2953
2954
2955 nbinsOFF = fHistOFF -> GetNbinsX();
2956 nbinsOFFNormalized = nbinsOFF;
2957 xlow = fHistOFF -> GetBinLowEdge(1);
2958 xup = fHistOFF -> GetBinLowEdge(nbinsOFFNormalized);
2959 xup = xup + BinWidthAlphaON;
2960
2961 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; Limits for fHistOFFNormalized: "
2962 << "nbins, xlow, xup: " << nbinsOFFNormalized << ", "
2963 << xlow << ", " << xup << endl;
2964
2965 fHistOFFNormalized = new TH1F (fHistOFFNormalizedName, fHistOFFNormalizedName,
2966 nbinsOFFNormalized, xlow, xup);
2967
2968
2969 // fHistOFFNormalized is filled with data from fHistOFF,
2970 // taken into account the possible bin difference
2971
2972
2973 for (Int_t i = 1; i <= nbinsOFF; i++)
2974 {
2975 BinCenter = fHistOFF -> GetBinCenter(i);
2976 fHistOFFNormalized -> Fill (BinCenter, fHistOFF -> GetBinContent(i));
2977 fHistOFFNormalized -> SetBinError(i, fHistOFF -> GetBinError(i));
2978 }
2979
2980
2981
2982
2983 for (Int_t i = 1; i <= nbinsOFFNormalized; i++)
2984 {
2985 content = fNormFactor * fHistOFFNormalized -> GetBinContent(i);
2986 error = fNormFactor * fHistOFFNormalized -> GetBinError(i);
2987
2988 fHistOFFNormalized -> SetBinContent (i, content);
2989 fHistOFFNormalized -> SetBinError (i, error);
2990 }
2991
2992
2993 // Number of entries is obtained from histOFF.
2994 // and set to histOFFNoramlized; otherwise, the number
2995 // of entries in histOFFNoramlized would be "nbins"
2996
2997 Double_t entries = fNormFactor * (fHistOFF -> GetEntries());
2998 fHistOFFNormalized -> SetEntries(entries);
2999
3000
3001
3002
3003 // If polynomial fit has been performed for fHistOFF,
3004 // it is defined a new polyfunction for fHistOFFNormalized,
3005 // which will be the polyfunction of fHistOFF normalized.
3006 // Function will be added to the function list of fHistOFFNormalized
3007
3008
3009
3010 if (fPolyOFF == NULL)
3011 {
3012 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fPolyOFF does not exist..."
3013 << endl;
3014 }
3015
3016 if (fPolyOFF)
3017 {
3018 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fPolyOFF exists... will be also normalized and included in the list of functions of fHistOFFNormalized"
3019 << endl;
3020
3021
3022
3023
3024 // Normalization of the function using fNormFactor and
3025 // BinWidthON/BinWidthOFF relation of alpha ON/OFF histograms
3026 // This makes possible to plot it together with ON alpha histo
3027
3028 TString FunctionName("PolyOFFNormalized");
3029
3030 Double_t xmin = fAlphaminOFF;
3031 Double_t xmax = fAlphamaxOFF;
3032
3033 TString formula = "[0]";
3034 TString bra1 = "+[";
3035 TString bra2 = "]";
3036 TString xpower = "*x";
3037 TString newpower = "*x";
3038 for (Int_t i=1; i<=fDegreeOFF; i++)
3039 {
3040 formula += bra1;
3041 formula += i;
3042 formula += bra2;
3043 formula += xpower;
3044
3045 xpower += newpower;
3046 }
3047
3048 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; formula = " << formula << endl;
3049
3050
3051
3052 fPolyOFFNormalized = new TF1 (FunctionName, formula, xmin, xmax);
3053
3054
3055 Double_t Parameter = 0.0;
3056 Double_t ParameterError = 0.0;
3057
3058 *fLog << " MHFindsignificanceONOFF::ComputeHistOFFNormalized; Fit parameters info: " << endl;
3059 for (Int_t i = 0; i <= fDegreeOFF; i++)
3060 {
3061 Parameter = fNormFactor * BinWidthRatioONOFF * fValuesOFF[i];
3062 ParameterError = fNormFactor * BinWidthRatioONOFF * fErrorsOFF[i];
3063
3064 fPolyOFFNormalized -> SetParameter(i, Parameter);
3065 fPolyOFFNormalized -> SetParError(i,ParameterError);
3066
3067
3068
3069 // Parameters are shown :
3070
3071 *fLog << " fValuesOFF[" << i<< "] = " << fValuesOFF[i]
3072 << " ; Parameter for fPolyOFFNormalized = " << Parameter << endl;
3073
3074 }
3075
3076
3077 TList *funclist = fHistOFFNormalized->GetListOfFunctions();
3078
3079 // temporal...
3080 //*fLog << "INFO concerning list of functions of fHistOFFNormalized :" << endl
3081 // << "List before adding OFF Normal., after adding it and after removing fPolyOFF..."
3082 // << endl;
3083
3084 //funclist-> Print();
3085 funclist-> Add(fPolyOFFNormalized);
3086
3087 //funclist-> Print();
3088
3089
3090
3091 }
3092
3093 return kTRUE;
3094
3095}
3096
3097// --------------------------------------------------------------------------
3098//
3099
3100Bool_t MHFindSignificanceONOFF::DrawHistOFF()
3101{
3102 if (fHistOFF == NULL )
3103 {
3104 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; fHistOFF = NULL" << endl;
3105 return kFALSE;
3106 }
3107
3108// fPsFilename -> NewPage();
3109
3110 // PLOT DISABLE
3111 /*
3112 TCanvas* CanvasHistOFF = new TCanvas(fHistOFF->GetName(), fHistOFF->GetName(), 600, 600);
3113
3114 //gStyle->SetOptFit(1011);
3115
3116 gROOT->SetSelectedPad(NULL);
3117 gStyle->SetPadLeftMargin(0.1);
3118 gStyle -> SetOptStat(1);
3119
3120 CanvasHistOFF->cd();
3121
3122
3123 if (fHistOFF)
3124 {
3125 fHistOFF->DrawCopy();
3126 }
3127
3128 // TF1 *fpolyOFF = fHistOFF->GetFunction("PolyOFF");
3129 if (fPolyOFF == NULL)
3130 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; fpolyOFF = NULL" << endl;
3131
3132 if (fPolyOFF)
3133 {
3134 // 2, 1 is red and solid
3135 fPolyOFF->SetLineColor(2);
3136 fPolyOFF->SetLineStyle(1);
3137 fPolyOFF->SetLineWidth(2);
3138 fPolyOFF->DrawCopy("same");
3139 }
3140
3141 CanvasHistOFF -> Update();
3142
3143
3144 */
3145
3146 return kTRUE;
3147}
3148
3149
3150// --------------------------------------------------------------------------
3151//
3152
3153Bool_t MHFindSignificanceONOFF::DrawHistOFFNormalized()
3154{
3155 if (fHistOFFNormalized == NULL )
3156 {
3157 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; fHistOFFNormalized = NULL" << endl;
3158 return kFALSE;
3159 }
3160
3161 // fPsFilename -> NewPage();
3162
3163 // PLOT DISABLE TO PERFORM GRID ANALYSIS
3164 /*
3165 TCanvas* CanvasHistOFFNormalized = new TCanvas(fHistOFFNormalized->GetName(),
3166 fHistOFFNormalized->GetName(), 600, 600);
3167
3168 //gStyle->SetOptFit(1011);
3169
3170 // gROOT->SetSelectedPad(NULL);
3171 gStyle->SetPadLeftMargin(0.1);
3172 gStyle -> SetOptStat(1);
3173
3174 CanvasHistOFFNormalized->cd();
3175
3176
3177 if (fHistOFFNormalized)
3178 {
3179 fHistOFFNormalized->DrawCopy();
3180 }
3181
3182 // TF1 *fpolyOFFNormalized = fHistOFFNormalized->GetFunction("PolyOFFNormalized");
3183 if (fPolyOFFNormalized == NULL)
3184 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; fPolyOFFNormalized = NULL" << endl;
3185
3186 if (fPolyOFFNormalized)
3187 {
3188 // 2, 1 is red and solid
3189 fPolyOFFNormalized->SetLineColor(2);
3190 fPolyOFFNormalized->SetLineStyle(1);
3191 fPolyOFFNormalized->SetLineWidth(2);
3192 fPolyOFFNormalized->DrawCopy("same");
3193 }
3194
3195 CanvasHistOFFNormalized -> Update();
3196*/
3197
3198 return kTRUE;
3199
3200}
3201
3202
3203// --------------------------------------------------------------------------
3204//
3205
3206Bool_t MHFindSignificanceONOFF::DrawFit(const Option_t *opt)
3207{
3208 if (fHistOFFNormalized == NULL || fHist == NULL)
3209 *fLog << "MHFindSignificanceONOFF::DrawFit; fHistOFFNormalized = NULL or fHist == NULL" << endl;
3210
3211 //fPsFilename -> NewPage();
3212
3213
3214 // PLOT DISABLE TO PERFORM GRID ANALYSIS
3215 // I DO SAVE PS FILE, BUT CANVAS IS DELETED AFTERWARDS
3216
3217 fCanvas = new TCanvas("Alpha", "Alpha plot", 600, 600);
3218 fCanvas -> SetFillColor(10);
3219
3220
3221
3222
3223 //gStyle->SetOptFit(1011);
3224
3225 gROOT->SetSelectedPad(NULL);
3226 gStyle -> SetFrameFillColor(10);
3227 gStyle->SetPadLeftMargin(0.15);
3228 gStyle -> SetOptStat(1);
3229
3230 fCanvas->cd();
3231
3232 if (fHist)
3233 {
3234 fHist -> SetTitle("Alpha Plot");
3235 fHist-> SetTitleOffset(1.5, "Y");
3236 fHist-> SetLineColor(kRed);
3237 fHist-> DrawCopy();
3238
3239 }
3240
3241
3242 if (fHistOFFNormalized)
3243 {
3244 TF1 *fpoly = fHistOFFNormalized->GetFunction("PolyOFFNormalized");
3245 if (fpoly == NULL)
3246 *fLog << "MHFindSignificanceONOFF::DrawFit; fPolyOFFNormalized = NULL" << endl;
3247
3248 if (fpoly)
3249 {
3250 // 2, 1 is red and solid
3251 fpoly->SetLineColor(2);
3252 fpoly->SetLineStyle(1);
3253 fpoly->SetLineWidth(2);
3254 fpoly->DrawCopy("same");
3255 }
3256
3257
3258// fHistOFFNormalized->SetLineColor(kRed);
3259// fHistOFFNormalized->DrawCopy("same");
3260
3261 }
3262
3263 if (fFitGauss)
3264 {
3265 TF1 *fpolygauss = fHist->GetFunction("PolyGauss");
3266 if (fpolygauss == NULL)
3267 *fLog << "MHFindSignificanceONOFF::DrawFit; fpolygauss = NULL" << endl;
3268
3269 if (fpolygauss)
3270 {
3271 // 4, 1 is blue and solid
3272 fpolygauss->SetLineColor(4);
3273 fpolygauss->SetLineStyle(1);
3274 fpolygauss->SetLineWidth(4);
3275 fpolygauss->DrawCopy("same");
3276 }
3277
3278 TF1 *fbackg = fHist->GetFunction("Backg");
3279 if (fbackg == NULL)
3280 *fLog << "MHFindSignificanceONOFF::DrawFit; fbackg = NULL" << endl;
3281
3282 if (fHistOFFNormalized)
3283 {
3284 fHistOFFNormalized->SetFillStyle(3002);
3285 fHistOFFNormalized->SetFillColor(9);
3286 fHistOFFNormalized->DrawCopy("Psamebar");}
3287
3288
3289 if (fbackg)
3290 {
3291 // 10, 1 is white and solid
3292 //fbackg->SetLineColor(kRed);
3293 fbackg->SetLineColor(196);
3294 fbackg->SetLineStyle(1);
3295 fbackg->SetLineWidth(4);
3296 fbackg->DrawCopy("same"); // I do not want to draw it... already too many things.
3297 }
3298 }
3299
3300// fCanvas->SaveAs("plot.root");
3301// fCanvas->SaveAs("plot.ps");
3302
3303 //-------------------------------
3304 // print results onto the figure
3305
3306
3307
3308 TPaveText *pt = new TPaveText(0.30, 0.35, 0.70, 0.90, "NDC");
3309 char tx[100];
3310
3311 sprintf(tx, "Results of polynomial fit to OFF (order %2d) :", fDegreeOFF);
3312 TText *t1 = pt->AddText(tx);
3313 t1->SetTextSize(0.03);
3314 t1->SetTextColor(2);
3315
3316 sprintf(tx, " (%6.2f< |alpha| <%6.2f [\\circ])", fAlphaminOFF, fAlphamaxOFF);
3317 pt->AddText(tx);
3318
3319 sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f",
3320 fChisqOFF, fNdfOFF, fProbOFF);
3321 pt->AddText(tx);
3322
3323 sprintf(tx, " OFF events (fit)= %8.1f #pm %8.1f",
3324 fNoffTotFitted, fdNoffTotFitted);
3325 pt->AddText(tx);
3326
3327 sprintf(tx, " OFF events (meas) = %8.1f #pm %8.1f", fNoffTot, fdNoffTot);
3328 pt->AddText(tx);
3329
3330 sprintf(tx, " OFF Normalization Factor (= Non/Noff) = %4.4f", fNormFactor);
3331 pt->AddText(tx);
3332
3333
3334
3335
3336 //sprintf(tx, " ");
3337 //pt->AddText(tx);
3338
3339 //--------------
3340 sprintf(tx, "Results for |alpha|< %6.2f [\\circ] :", fAlphasi);
3341 TText *t6 = pt->AddText(tx);
3342 t6->SetTextSize(0.03);
3343 t6->SetTextColor(8);
3344
3345 sprintf(tx, " Non = %8.1f #pm %8.1f", fNon, fdNon);
3346 pt->AddText(tx);
3347
3348
3349 if(fUseFittedQuantities)
3350 {
3351 //// **************************************************
3352 ///// PRINT INFORMATION ABOUT FITTED QUANTITIES /////////
3353
3354
3355
3356 Double_t NoffFitNormalized = fNoffSigFitted * fNormFactor;
3357 Double_t ErrorNoffFitNormalized = fdNoffSigFitted * fNormFactor;
3358 Double_t SignificanceUsed = GetSignificance();
3359
3360 sprintf(tx, " Noff Fitted (Normalized) = %8.1f #pm %8.1f",
3361 NoffFitNormalized, ErrorNoffFitNormalized);
3362 pt->AddText(tx);
3363
3364
3365 sprintf(tx, " Nex (ON - OFF Fitted) = %8.1f #pm %8.1f",
3366 fNexONOFFFitted, fdNexONOFFFitted);
3367 pt->AddText(tx);
3368
3369
3370 sprintf(tx, " Gamma = %4.4f, Effective Noff (i.e. fNoff) = %6.1f",
3371 fGamma, fNoff);
3372 pt->AddText(tx);
3373
3374
3375 Double_t ratio = fNoffSigFitted>0.0 ? fNexONOFFFitted/(fNoffSigFitted*fNormFactor) : 0.0;
3376 sprintf(tx, " Significance = %6.2f, Nex/(Nbg*NormFactor) = %6.2f",
3377 SignificanceUsed, ratio);
3378 pt->AddText(tx);
3379
3380
3381 }
3382
3383 else
3384 {
3385 //// **************************************************
3386 ///// PRINT INFORMATION ABOUT MEASURED QUANTITIES /////////
3387
3388
3389 Double_t NoffNormalized = fNoffSig * fNormFactor;
3390 Double_t ErrorNoffNormalized = fdNoffSig * fNormFactor;
3391 Double_t SignificanceUsed = GetSignificance();
3392
3393 sprintf(tx, " Noff measured (Normalized) = %8.1f #pm %8.1f",
3394 NoffNormalized, ErrorNoffNormalized);
3395 pt->AddText(tx);
3396
3397 sprintf(tx, " Nex (ON - OFF measured) = %8.1f #pm %8.1f",
3398 fNexONOFF, fdNexONOFF);
3399 pt->AddText(tx);
3400
3401 Double_t ratio = fNoffSig>0.0 ? fNexONOFF/(fNoffSig*fNormFactor) : 0.0;
3402 sprintf(tx, " Significance = %6.2f, Nex/(Nbg*NormFactor) = %6.2f",
3403 SignificanceUsed, ratio);
3404 pt->AddText(tx);
3405
3406 }
3407
3408 /*
3409 // Temporally I will also show ALL SIGMALIMA COMPUTED.
3410
3411 sprintf(tx,
3412 " fSigLiMa1 = %6.2f, fSigLiMa2 = %6.2f, fSigLiMa3 = %6.2f",
3413 fSigLiMa,fSigLiMa2, fSigLiMa3);
3414 pt->AddText(tx);
3415 */
3416
3417
3418 //--------------
3419 if (fFitGauss)
3420 {
3421 sprintf(tx, "Results of (polynomial+Gauss) fit :");
3422 TText *t7 = pt->AddText(tx);
3423 t7->SetTextSize(0.03);
3424 t7->SetTextColor(4);
3425
3426 sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f",
3427 fGChisq, fGNdf, fGProb);
3428 pt->AddText(tx);
3429
3430 sprintf(tx, " Sigma of Gauss = %8.1f #pm %8.1f [\\circ]",
3431 fSigmaGauss, fdSigmaGauss);
3432 pt->AddText(tx);
3433
3434 sprintf(tx, " total no.of excess events = %8.1f #pm %8.1f",
3435 fNexGauss, fdNexGauss);
3436 pt->AddText(tx);
3437 }
3438 //--------------
3439
3440 pt->SetFillStyle(0);
3441 pt->SetBorderSize(0);
3442 pt->SetTextAlign(12);
3443
3444
3445 if(fPrintResultsOntoAlphaPlot)
3446 {
3447 pt->Draw();
3448 }
3449 fCanvas->Modified();
3450 fCanvas->Update();
3451
3452 // fPsFilename -> NewPage();
3453
3454
3455
3456 if (fSavePlots)
3457 {
3458 // ********************************************
3459 // TMP solution while the TPostScript thing is not working.
3460 // PsFileName for storing these histograms is derived
3461 // from fPsFilenameString.
3462
3463
3464 cout << "Alpha plot with ON-OFF data will be saved in PostScript file " ;
3465
3466
3467 if (!fPsFilenameString.IsNull())
3468 {
3469 TString filename = (fPsFilenameString);
3470 TString filename2 = (fPsFilenameString);
3471 // Train or Test Sample is specified outside
3472 // class MHFindSignificanceONOFF, and included in
3473 // fPsFilenameString
3474
3475 filename += ("AlphaPlot.ps");
3476 filename2 += ("AlphaPlot.C");
3477 cout << filename << endl;
3478 fCanvas -> SaveAs(filename);
3479 fCanvas -> SaveAs(filename2);
3480 }
3481
3482 // END OF TEMPORAL SOLUTION
3483 // ********************************************
3484
3485 }
3486
3487 // Canvvas deleted to allow for GRID analysis
3488
3489 delete fCanvas;
3490
3491
3492 return kTRUE;
3493}
3494
3495
3496// --------------------------------------------------------------------------
3497//
3498// Print the results of the polynomial fit to the alpha OFF distribution
3499//
3500//
3501void MHFindSignificanceONOFF::PrintPolyOFF(Option_t *o)
3502{
3503 *fLog << "---------------------------" << endl;
3504 *fLog << "MHFindSignificanceONOFF::PrintPolyOFF :" << endl;
3505
3506 *fLog << "fAlphaminOFF, fAlphamaxOFF, fDegreeOFF "
3507 << fAlphaminOFF << ", " << fAlphamaxOFF << ", " << fDegreeOFF << endl;
3508
3509 *fLog << "fMbinsOFF, fNzeroOFF, fIstatOFF = " << fMbinsOFF << ", "
3510 << fNzeroOFF << ", " << fIstatOFF << endl;
3511
3512 *fLog << "fChisqOFF, fNdfOFF, fProbOFF = " << fChisqOFF << ", "
3513 << fNdfOFF << ", " << fProbOFF << endl;
3514
3515 *fLog << "fNon; fNoffSigFitted, fdNoffSigFitted; fNoffSig, fdNoffSig = "
3516 << fNon << "; " << fNoffSigFitted << ", " << fdNoffSigFitted
3517 << "; " << fNoffSig << ", " << fdNoffSig << endl;
3518
3519 Double_t sigtoback = fNoffSigFitted >0.0 ? fNexONOFFFitted/(fNoffSigFitted*fNormFactor) : 0.0;
3520
3521
3522 *fLog << "fNexONOFFFitted, fdNexONOFFFitted, fGamma, fNoff, fSigLiMa, sigtoback = "
3523 << fNexONOFFFitted << ", " << fdNexONOFFFitted << ", "
3524 << fGamma << ", " << fNoff
3525 << ", " << fSigLiMa << ", "
3526 << sigtoback << endl;
3527
3528
3529 Double_t sigtoback2 = fNoffSig >0.0 ? fNexONOFF/(fNoffSig*fNormFactor) : 0.0;
3530
3531 *fLog << "fNexONOFF, fdNexONOFF, fNormFactor, fNoffSig, fSigLiMa2, sigtoback2 = "
3532 << fNexONOFF << ", " << fdNexONOFF << ", "
3533 << fNormFactor << ", " << fNoffSig
3534 << ", " << fSigLiMa2 << ", "
3535 << sigtoback2 << endl;
3536
3537 *fLog << "---------------------------" << endl;
3538}
3539
3540// --------------------------------------------------------------------------
3541//
3542// Print the results of the polynomial fit to the alpha distribution
3543//
3544//
3545void MHFindSignificanceONOFF::PrintPoly(Option_t *o)
3546{
3547 *fLog << "---------------------------" << endl;
3548 *fLog << "MHFindSignificanceONOFF::PrintPoly :" << endl;
3549
3550 *fLog << "fAlphami, fAlphama, fDegree, fAlphasi = "
3551 << fAlphami << ", " << fAlphama << ", " << fDegree << ", "
3552 << fAlphasi << endl;
3553
3554 *fLog << "fMbins, fNzero, fIstat = " << fMbins << ", "
3555 << fNzero << ", " << fIstat << endl;
3556
3557 *fLog << "fChisq, fNdf, fProb = " << fChisq << ", "
3558 << fNdf << ", " << fProb << endl;
3559
3560 *fLog << "fNon, fNbg, fdNbg, fNbgtot, fNbgtotFitted, fdNbgtotFitted = "
3561 << fNon << ", " << fNbg << ", " << fdNbg << ", " << fNbgtot
3562 << ", " << fNbgtotFitted << ", " << fdNbgtotFitted << endl;
3563
3564 Double_t sigtoback = fNbg>0.0 ? fNex/fNbg : 0.0;
3565 *fLog << "fNex, fdNex, fGamma, fNoff, fSigLiMa, sigtoback = "
3566 << fNex << ", " << fdNex << ", " << fGamma << ", " << fNoff
3567 << ", " << fSigLiMa << ", " << sigtoback << endl;
3568
3569 //------------------------------------
3570 // get errors
3571
3572 /*
3573 Double_t eplus;
3574 Double_t eminus;
3575 Double_t eparab;
3576 Double_t gcc;
3577 Double_t errdiag;
3578
3579
3580 if ( !fConstantBackg )
3581 {
3582 *fLog << "parameter value error eplus eminus eparab errdiag gcc"
3583 << endl;
3584
3585 for (Int_t j=0; j<=fDegree; j++)
3586 {
3587 if (gMinuit)
3588 gMinuit->mnerrs(j, eplus, eminus, eparab, gcc);
3589 errdiag = sqrt(fEma[j][j]);
3590 *fLog << j << " " << fValues[j] << " " << fErrors[j] << " "
3591 << eplus << " " << eminus << " " << eparab << " "
3592 << errdiag << " " << gcc << endl;
3593 }
3594 }
3595 else
3596 {
3597 *fLog << "parameter value error errdiag "
3598 << endl;
3599
3600 for (Int_t j=0; j<=fDegree; j++)
3601 {
3602 errdiag = sqrt(fEma[j][j]);
3603 *fLog << j << " " << fValues[j] << " " << fErrors[j] << " "
3604 << errdiag << endl;
3605 }
3606 }
3607 */
3608
3609 //----------------------------------------
3610 /*
3611 *fLog << "Covariance matrix :" << endl;
3612 for (Int_t j=0; j<=fDegree; j++)
3613 {
3614 *fLog << "j = " << j << " : ";
3615 for (Int_t k=0; k<=fDegree; k++)
3616 {
3617 *fLog << fEma[j][k] << " ";
3618 }
3619 *fLog << endl;
3620 }
3621
3622 *fLog << "Correlation matrix :" << endl;
3623 for (Int_t j=0; j<=fDegree; j++)
3624 {
3625 *fLog << "j = " << j << " : ";
3626 for (Int_t k=0; k<=fDegree; k++)
3627 {
3628 *fLog << fCorr[j][k] << " ";
3629 }
3630 *fLog << endl;
3631 }
3632 */
3633
3634 *fLog << "---------------------------" << endl;
3635}
3636
3637// --------------------------------------------------------------------------
3638//
3639// Print the results of the (polynomial+Gauss) fit to the alpha distribution
3640//
3641//
3642void MHFindSignificanceONOFF::PrintPolyGauss(Option_t *o)
3643{
3644 *fLog << "---------------------------" << endl;
3645 *fLog << "MHFindSignificanceONOFF::PrintPolyGauss :" << endl;
3646
3647 *fLog << "fAlphalo, fAlphahi = "
3648 << fAlphalo << ", " << fAlphahi << endl;
3649
3650 *fLog << "fGMbins, fGNzero, fGIstat = " << fGMbins << ", "
3651 << fGNzero << ", " << fGIstat << endl;
3652
3653 *fLog << "fGChisq, fGNdf, fGProb = " << fGChisq << ", "
3654 << fGNdf << ", " << fGProb << endl;
3655
3656
3657 //------------------------------------
3658 // get errors
3659
3660 Double_t eplus;
3661 Double_t eminus;
3662 Double_t eparab;
3663 Double_t gcc;
3664 Double_t errdiag;
3665
3666 *fLog << "parameter value error eplus eminus eparab errdiag gcc"
3667 << endl;
3668 for (Int_t j=0; j<=(fDegree+3); j++)
3669 {
3670 if (gMinuit)
3671 gMinuit->mnerrs(j, eplus, eminus, eparab, gcc);
3672 errdiag = sqrt(fGEma[j][j]);
3673 *fLog << j << " " << fGValues[j] << " " << fGErrors[j] << " "
3674 << eplus << " " << eminus << " " << eparab << " "
3675 << errdiag << " " << gcc << endl;
3676 }
3677
3678
3679 *fLog << "Covariance matrix :" << endl;
3680 for (Int_t j=0; j<=(fDegree+3); j++)
3681 {
3682 *fLog << "j = " << j << " : ";
3683 for (Int_t k=0; k<=(fDegree+3); k++)
3684 {
3685 *fLog << fGEma[j][k] << " ";
3686 }
3687 *fLog << endl;
3688 }
3689
3690 *fLog << "Correlation matrix :" << endl;
3691 for (Int_t j=0; j<=(fDegree+3); j++)
3692 {
3693 *fLog << "j = " << j << " : ";
3694 for (Int_t k=0; k<=(fDegree+3); k++)
3695 {
3696 *fLog << fGCorr[j][k] << " ";
3697 }
3698 *fLog << endl;
3699 }
3700
3701 *fLog << "---------------------------" << endl;
3702}
3703
3704//============================================================================
3705
3706
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718
3719
3720
3721
3722
3723
3724
3725
3726
3727
3728
3729
3730
3731
3732
3733
3734
3735
3736
3737
3738
3739
3740
Note: See TracBrowser for help on using the repository browser.