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

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