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

Last change on this file since 5140 was 4411, checked in by paneque, 20 years ago
*** empty log message ***
File size: 94.7 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 fNzeroOFF++;
1050
1051 // set minimum error
1052 if (content < 9.0)
1053 {
1054 fMlowOFF += 1;
1055 fHistOFF->SetBinError(i, 3.0);
1056 }
1057
1058 //*fLog << "Take : i, content, error = " << i << ", "
1059 // << fHist->GetBinContent(i) << ", "
1060 // << fHist->GetBinError(i) << endl;
1061
1062 continue;
1063 }
1064 // bin is not completely contained in the fit range : set error huge
1065
1066 fHistOFF->SetBinError(i, dummy);
1067
1068 //*fLog << "Omit : i, content, error = " << i << ", "
1069 // << fHist->GetBinContent(i) << ", " << fHist->GetBinError(i)
1070 // << endl;
1071
1072 }
1073
1074 // mean of entries/bin in the fit range
1075 if (fMbinsOFF > 0)
1076 {
1077 mean /= ((Double_t) fMbinsOFF);
1078 rms /= ((Double_t) fMbinsOFF);
1079 }
1080
1081 rms = sqrt( rms - mean*mean );
1082
1083 // if there are no events in the background region
1084 // there is no reason for rebinning
1085 // and this is the condition for assuming a constant background (= 0)
1086 if (mean <= 0.0)
1087 break;
1088
1089 Double_t helpmi = fAlphami*fAlphami*fAlphami;
1090 Double_t helpmm = fAlphamm*fAlphamm*fAlphamm;
1091 Double_t helpma = fAlphama*fAlphama*fAlphama;
1092 Double_t help = (helpma-helpmm) * (fAlphamm-fAlphami)
1093 - (helpmm-helpmi) * (fAlphama-fAlphamm);
1094 if (help != 0.0)
1095 a2init = ( (fAlphamm-fAlphami)*nfar - (fAlphama-fAlphamm)*nclose )
1096 * 1.5 * fHistOFF->GetBinWidth(1) / help;
1097 else
1098 a2init = 0.0;
1099
1100
1101 //--------------------------------------------
1102 // rebin the histogram
1103 // - if a bin has no entries
1104 // - or if there are too many bins with too few entries
1105 // - or if the new bin width would exceed half the size of the
1106 // signal region
1107
1108 if ( !fRebin ||
1109 ( fNzeroOFF <= 0 && (Double_t)fMlowOFF<0.05*(Double_t)fMbinsOFF ) ||
1110 (Double_t)(nrebin+1)/(Double_t)nrebin * fHistOFF->GetBinWidth(1)
1111 > fAlphasig/2.0 )
1112 {
1113 //*fLog << "before break" << endl;
1114 break;
1115 }
1116
1117 nrebin += 1;
1118 TString histname = fHistOFF->GetName();
1119 delete fHistOFF;
1120 fHistOFF = NULL;
1121
1122 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; rebin the |alpha|OFF plot, grouping "
1123 << nrebin << " bins together" << endl;
1124
1125 // TH1::Rebin doesn't work properly
1126 //fHist = fHistOrig->Rebin(nrebin, "Rebinned");
1127 // use private routine RebinHistogram()
1128 fHistOFF = new TH1F;
1129 fHistOFF->Sumw2();
1130 fHistOFF->SetNameTitle(histname, histname);
1131 fHistOFF->UseCurrentStyle();
1132
1133 // do rebinning such that x0 remains a lower bin edge
1134 Double_t x0 = 0.0;
1135 if ( !RebinHistogramOFF(x0, nrebin) )
1136 {
1137 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; RebinHistgramOFF() failed"
1138 << endl;
1139 return kFALSE;
1140 }
1141
1142 fHistOFF->SetXTitle("|alpha| [\\circ]");
1143 fHistOFF->SetYTitle("Counts");
1144
1145 }
1146 //---------------- end of while loop for rebinning -----------------
1147
1148
1149 // if there are still too many bins with too few entries don't fit
1150 // and assume a constant background
1151
1152 fConstantBackg = kFALSE;
1153 if ( fNzeroOFF > 0 || (Double_t)fMlowOFF>0.05*(Double_t)fMbinsOFF )
1154 {
1155 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; polynomial fit not possible, fNzeroOFF, fMlowOFF, fMbinsOFF = "
1156 << fNzeroOFF << ", " << fMlowOFF << ", " << fMbinsOFF << endl;
1157 *fLog << " assume a constant background" << endl;
1158
1159 fConstantBackg = kTRUE;
1160 fDegreeOFF = 0;
1161
1162 TString funcname = "PolyOFF";
1163 Double_t xmin = 0.0;
1164 Double_t xmax = 90.0;
1165
1166 TString formula = "[0]";
1167
1168 fPolyOFF = new TF1(funcname, formula, xmin, xmax);
1169 TList *funclist = fHistOFF->GetListOfFunctions();
1170 funclist->Add(fPolyOFF);
1171
1172 //--------------------
1173 Int_t nparfree = 1;
1174 fChisqOFF = 0.0;
1175 fNdfOFF = fMbinsOFF - nparfree;
1176 fProbOFF = 0.0;
1177 fIstatOFF = 0;
1178
1179 fValuesOFF.Set(1);
1180 fErrorsOFF.Set(1);
1181
1182 Double_t val, err;
1183 val = mean;
1184 err = sqrt( mean / (Double_t)fMbinsOFF );
1185
1186 fPolyOFF->SetParameter(0, val);
1187 fPolyOFF->SetParError (0, err);
1188
1189 fValuesOFF[0] = val;
1190 fErrorsOFF[0] = err;
1191
1192 fEmaOFF[0][0] = err*err;
1193 fCorrOFF[0][0] = 1.0;
1194 //--------------------
1195 //--------------------------------------------------
1196 // reset the errors of the points in the histogram
1197 for (Int_t i=1; i<=nbins; i++)
1198 {
1199 fHistOFF->SetBinError(i, saveError[i-1]);
1200 }
1201
1202
1203 return kTRUE;
1204 }
1205
1206
1207 //=========== start loop for reducing the degree ==================
1208 // of the polynomial
1209 while (1)
1210 {
1211 //--------------------------------------------------
1212 // prepare fit of a polynomial : (a0 + a1*x + a2*x**2 + a3*x**3 + ...)
1213
1214 TString funcname = "PolyOFF";
1215 Double_t xmin = 0.0;
1216 Double_t xmax = 90.0;
1217
1218 TString formula = "[0]";
1219 TString bra1 = "+[";
1220 TString bra2 = "]";
1221 TString xpower = "*x";
1222 TString newpower = "*x";
1223 for (Int_t i=1; i<=fDegreeOFF; i++)
1224 {
1225 formula += bra1;
1226 formula += i;
1227 formula += bra2;
1228 formula += xpower;
1229
1230 xpower += newpower;
1231 }
1232
1233 //*fLog << "FitPolynomial : formula = " << formula << endl;
1234
1235 fPolyOFF = new TF1(funcname, formula, xmin, xmax);
1236 TList *funclist = fHistOFF->GetListOfFunctions();
1237 funclist->Add(fPolyOFF);
1238
1239 //------------------------
1240 // attention : the dimensions must agree with those in CallMinuit()
1241 const UInt_t npar = fDegreeOFF+1;
1242
1243 TString parname[npar];
1244 TArrayD vinit(npar);
1245 TArrayD step(npar);
1246 TArrayD limlo(npar);
1247 TArrayD limup(npar);
1248 TArrayI fix(npar);
1249
1250 vinit[0] = mean;
1251 vinit[2] = a2init;
1252
1253 for (UInt_t j=0; j<npar; j++)
1254 {
1255 parname[j] = "p";
1256 parname[j] += j+1;
1257
1258 step[j] = vinit[j] != 0.0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
1259 }
1260
1261 // limit the first coefficient of the polynomial to positive values
1262 // because the background must not be negative
1263 limup[0] = fHistOFF->GetEntries();
1264
1265 // use the subsequernt loop if you want to apply the
1266 // constraint : uneven derivatives (at alpha=0) = zero
1267 for (UInt_t j=1; j<npar; j+=2)
1268 {
1269 vinit[j] = 0;
1270 step[j] = 0;
1271 fix[j] = 1;
1272 }
1273
1274 //*fLog << "FitPolynomial : before CallMinuit()" << endl;
1275
1276 MMinuitInterface inter;
1277 const Bool_t rc = inter.CallMinuit(fcnpolyOFF, parname, vinit, step,
1278 limlo, limup, fix, fHistOFF, "Migrad",
1279 kFALSE);
1280
1281 //*fLog << "FitPolynomial : after CallMinuit()" << endl;
1282
1283 if (rc != 0)
1284 {
1285 // *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit failed"
1286 // << endl;
1287 // return kFALSE;
1288 }
1289
1290 //-------------------
1291 // get status of minimization
1292 Double_t fmin = 0;
1293 Double_t fedm = 0;
1294 Double_t errdef = 0;
1295 Int_t npari = 0;
1296 Int_t nparx = 0;
1297
1298 if (gMinuit)
1299 gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fIstatOFF);
1300
1301 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; fmin, fedm, errdef, npari, nparx, fIstat = "
1302 << fmin << ", " << fedm << ", " << errdef << ", " << npari
1303 << ", " << nparx << ", " << fIstatOFF << endl;
1304
1305
1306 //-------------------
1307 // store the results
1308
1309 Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
1310 fChisqOFF = fmin;
1311 fNdfOFF = fMbinsOFF - nparfree;
1312 fProbOFF = TMath::Prob(fChisqOFF, fNdfOFF);
1313
1314
1315 // get fitted parameter values and errors
1316 fValuesOFF.Set(npar);
1317 fErrorsOFF.Set(npar);
1318
1319 for (Int_t j=0; j<=fDegreeOFF; j++)
1320 {
1321 Double_t val, err;
1322 if (gMinuit)
1323 gMinuit->GetParameter(j, val, err);
1324
1325 fPolyOFF->SetParameter(j, val);
1326 fPolyOFF->SetParError(j, err);
1327
1328 fValuesOFF[j] = val;
1329 fErrorsOFF[j] = err;
1330 }
1331
1332 // if the highest coefficient (j0) of the polynomial
1333 // is consistent with zero reduce the degree of the polynomial
1334
1335 Int_t j0 = 0;
1336 for (Int_t j=fDegreeOFF; j>1; j--)
1337 {
1338 // ignore fixed parameters
1339 if (fErrorsOFF[j] == 0)
1340 continue;
1341
1342 // this is the highest coefficient
1343 j0 = j;
1344 break;
1345 }
1346
1347 if (!fReduceDegree || j0==0 || TMath::Abs(fValuesOFF[j0]) > fErrorsOFF[j0])
1348 break;
1349
1350 // reduce the degree of the polynomial
1351 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; reduce the degree of the polynomial from "
1352 << fDegreeOFF << " to " << (j0-2) << endl;
1353 fDegreeOFF = j0 - 2;
1354
1355 funclist->Remove(fPolyOFF);
1356 //if (fPoly)
1357 delete fPolyOFF;
1358 fPolyOFF = NULL;
1359
1360 // delete the Minuit object in order to have independent starting
1361 // conditions for the next minimization
1362 //if (gMinuit)
1363 delete gMinuit;
1364 gMinuit = NULL;
1365 }
1366 //=========== end of loop for reducing the degree ==================
1367 // of the polynomial
1368
1369
1370 //--------------------------------------------------
1371 // get the error matrix of the fitted parameters
1372
1373 if (fIstatOFF >= 1)
1374 {
1375 // error matrix was calculated
1376 if (gMinuit)
1377 gMinuit->mnemat(&fEmatOFF[0][0], fNdimOFF);
1378
1379 // copy covariance matrix into a matrix which includes also the fixed
1380 // parameters
1381 TString name;
1382 Double_t bnd1, bnd2, val, err;
1383 Int_t jvarbl;
1384 Int_t kvarbl;
1385 for (Int_t j=0; j<=fDegreeOFF; j++)
1386 {
1387 if (gMinuit)
1388 gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
1389
1390 for (Int_t k=0; k<=fDegreeOFF; k++)
1391 {
1392 if (gMinuit)
1393 gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
1394
1395 fEmaOFF[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fEmatOFF[jvarbl-1][kvarbl-1];
1396 }
1397 }
1398 }
1399 else
1400 {
1401 // error matrix was not calculated, construct it
1402 *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; error matrix not defined"
1403 << endl;
1404 for (Int_t j=0; j<=fDegreeOFF; j++)
1405 {
1406 for (Int_t k=0; k<=fDegreeOFF; k++)
1407 fEmaOFF[j][k] = 0;
1408
1409 fEmaOFF[j][j] = fErrorsOFF[j]*fErrorsOFF[j];
1410 }
1411 }
1412
1413
1414
1415
1416 //--------------------------------------------------
1417 // calculate correlation matrix
1418 for (Int_t j=0; j<=fDegreeOFF; j++)
1419 {
1420 for (Int_t k=0; k<=fDegreeOFF; k++)
1421 {
1422 const Double_t sq = fEmaOFF[j][j]*fEmaOFF[k][k];
1423 fCorrOFF[j][k] =
1424 sq == 0 ? 0 : (fEmaOFF[j][k] / TMath::Sqrt(fEmaOFF[j][j]*fEmaOFF[k][k]));
1425 }
1426 }
1427
1428 //--------------------------------------------------
1429 // reset the errors of the points in the histogram
1430 for (Int_t i=1; i<=nbins; i++)
1431 {
1432 fHistOFF->SetBinError(i, saveError[i-1]);
1433 }
1434
1435 return kTRUE;
1436
1437
1438}
1439
1440
1441
1442
1443
1444
1445
1446
1447// --------------------------------------------------------------------------
1448//
1449// FitPolynomial
1450//
1451// - create a clone 'fHist' of the |alpha| distribution 'fHistOrig'
1452// - fit a polynomial of degree 'fDegree' to the alpha distribution
1453// 'fHist' in the region alphamin < |alpha| < alphamax
1454//
1455// in pathological cases the histogram is rebinned before fitting
1456// (this is done only if fRebin is kTRUE)
1457//
1458// if the highest coefficient of the polynomial is compatible with zero
1459// the fit is repeated with a polynomial of lower degree
1460// (this is done only if fReduceDegree is kTRUE)
1461//
1462//
1463
1464Bool_t MHFindSignificanceONOFF::FitPolynomial()
1465{
1466 //--------------------------------------------------
1467 // check the histogram :
1468 // - calculate initial values of the parameters
1469 // - check for bins with zero entries
1470 // - set minimum errors
1471 // - save the original errors
1472 // - set errors huge outside the fit range
1473 // (in 'fcnpoly' points with huge errors will be ignored)
1474
1475
1476 Double_t dummy = 1.e20;
1477
1478 Double_t mean;
1479 Double_t rms;
1480 Double_t nclose;
1481 Double_t nfar;
1482 Double_t a2init = 0.0;
1483 TArrayD saveError;
1484
1485 Int_t nbins;
1486 Int_t nrebin = 1;
1487
1488 //---------------- start while loop for rebinning -----------------
1489 while(1)
1490 {
1491
1492 fNzero = 0;
1493 fMbins = 0;
1494 fMlow = 0;
1495 fNbgtot = 0.0;
1496
1497 fAlphami = 10000.0;
1498 fAlphamm = 10000.0;
1499 fAlphama = -10000.0;
1500
1501 mean = 0.0;
1502 rms = 0.0;
1503 nclose = 0.0;
1504 nfar = 0.0;
1505
1506 nbins = fHist->GetNbinsX();
1507 saveError.Set(nbins);
1508
1509 for (Int_t i=1; i<=nbins; i++)
1510 {
1511 saveError[i-1] = fHist->GetBinError(i);
1512
1513 // bin should be completely contained in the fit range
1514 // (fAlphamin, fAlphamax)
1515 Double_t xlo = fHist->GetBinLowEdge(i);
1516 Double_t xup = fHist->GetBinLowEdge(i+1);
1517
1518 if ( xlo >= fAlphamin-fEps && xlo <= fAlphamax+fEps &&
1519 xup >= fAlphamin-fEps && xup <= fAlphamax+fEps )
1520 {
1521 fMbins++;
1522
1523 if ( xlo < fAlphami )
1524 fAlphami = xlo;
1525
1526 if ( xup > fAlphama )
1527 fAlphama = xup;
1528
1529 Double_t content = fHist->GetBinContent(i);
1530 fNbgtot += content;
1531
1532 mean += content;
1533 rms += content*content;
1534
1535 // count events in low-alpha and high-alpha region
1536 if ( xlo >= fAlphammm-fEps && xup >= fAlphammm-fEps)
1537 {
1538 nfar += content;
1539 if ( xlo < fAlphamm )
1540 fAlphamm = xlo;
1541 if ( xup < fAlphamm )
1542 fAlphamm = xup;
1543 }
1544 else
1545 {
1546 nclose += content;
1547 if ( xlo > fAlphamm )
1548 fAlphamm = xlo;
1549 if ( xup > fAlphamm )
1550 fAlphamm = xup;
1551 }
1552
1553 // count bins with zero entry
1554 if (content <= 0.0)
1555 fNzero++;
1556
1557 // set minimum error
1558 if (content < 9.0)
1559 {
1560 fMlow += 1;
1561 fHist->SetBinError(i, 3.0);
1562 }
1563
1564 //*fLog << "Take : i, content, error = " << i << ", "
1565 // << fHist->GetBinContent(i) << ", "
1566 // << fHist->GetBinError(i) << endl;
1567
1568 continue;
1569 }
1570 // bin is not completely contained in the fit range : set error huge
1571
1572 fHist->SetBinError(i, dummy);
1573
1574 //*fLog << "Omit : i, content, error = " << i << ", "
1575 // << fHist->GetBinContent(i) << ", " << fHist->GetBinError(i)
1576 // << endl;
1577
1578 }
1579
1580 // mean of entries/bin in the fit range
1581 if (fMbins > 0)
1582 {
1583 mean /= ((Double_t) fMbins);
1584 rms /= ((Double_t) fMbins);
1585 }
1586
1587 rms = sqrt( rms - mean*mean );
1588
1589 // if there are no events in the background region
1590 // there is no reason for rebinning
1591 // and this is the condition for assuming a constant background (= 0)
1592 if (mean <= 0.0)
1593 break;
1594
1595 Double_t helpmi = fAlphami*fAlphami*fAlphami;
1596 Double_t helpmm = fAlphamm*fAlphamm*fAlphamm;
1597 Double_t helpma = fAlphama*fAlphama*fAlphama;
1598 Double_t help = (helpma-helpmm) * (fAlphamm-fAlphami)
1599 - (helpmm-helpmi) * (fAlphama-fAlphamm);
1600 if (help != 0.0)
1601 a2init = ( (fAlphamm-fAlphami)*nfar - (fAlphama-fAlphamm)*nclose )
1602 * 1.5 * fHist->GetBinWidth(1) / help;
1603 else
1604 a2init = 0.0;
1605
1606
1607 //--------------------------------------------
1608 // rebin the histogram
1609 // - if a bin has no entries
1610 // - or if there are too many bins with too few entries
1611 // - or if the new bin width would exceed half the size of the
1612 // signal region
1613
1614 if ( !fRebin ||
1615 ( fNzero <= 0 && (Double_t)fMlow<0.05*(Double_t)fMbins ) ||
1616 (Double_t)(nrebin+1)/(Double_t)nrebin * fHist->GetBinWidth(1)
1617 > fAlphasig/2.0 )
1618 {
1619 //*fLog << "before break" << endl;
1620 break;
1621 }
1622
1623 nrebin += 1;
1624 TString histname = fHist->GetName();
1625 delete fHist;
1626 fHist = NULL;
1627
1628 *fLog << "MHFindSignificanceONOFF::FitPolynomial; rebin the |alpha| plot, grouping "
1629 << nrebin << " bins together" << endl;
1630
1631 // TH1::Rebin doesn't work properly
1632 //fHist = fHistOrig->Rebin(nrebin, "Rebinned");
1633 // use private routine RebinHistogram()
1634 fHist = new TH1F;
1635 fHist->Sumw2();
1636 fHist->SetNameTitle(histname, histname);
1637 fHist->UseCurrentStyle();
1638
1639 // do rebinning such that x0 remains a lower bin edge
1640 Double_t x0 = 0.0;
1641 if ( !RebinHistogram(x0, nrebin) )
1642 {
1643 *fLog << "MHFindSignificanceONOFF::FitPolynomial; RebinHistgram() failed"
1644 << endl;
1645 return kFALSE;
1646 }
1647
1648 fHist->SetXTitle("|alpha| [\\circ]");
1649 fHist->SetYTitle("Counts");
1650
1651 }
1652 //---------------- end of while loop for rebinning -----------------
1653
1654
1655 // if there are still too many bins with too few entries don't fit
1656 // and assume a constant background
1657
1658 fConstantBackg = kFALSE;
1659 if ( fNzero > 0 || (Double_t)fMlow>0.05*(Double_t)fMbins )
1660 {
1661 *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit not possible, fNzero, fMlow, fMbins = "
1662 << fNzero << ", " << fMlow << ", " << fMbins << endl;
1663 *fLog << " assume a constant background" << endl;
1664
1665 fConstantBackg = kTRUE;
1666 fDegree = 0;
1667
1668 TString funcname = "Poly";
1669 Double_t xmin = 0.0;
1670 Double_t xmax = 90.0;
1671
1672 TString formula = "[0]";
1673
1674 fPoly = new TF1(funcname, formula, xmin, xmax);
1675 TList *funclist = fHist->GetListOfFunctions();
1676 funclist->Add(fPoly);
1677
1678 //--------------------
1679 Int_t nparfree = 1;
1680 fChisq = 0.0;
1681 fNdf = fMbins - nparfree;
1682 fProb = 0.0;
1683 fIstat = 0;
1684
1685 fValues.Set(1);
1686 fErrors.Set(1);
1687
1688 Double_t val, err;
1689 val = mean;
1690 err = sqrt( mean / (Double_t)fMbins );
1691
1692 fPoly->SetParameter(0, val);
1693 fPoly->SetParError (0, err);
1694
1695 fValues[0] = val;
1696 fErrors[0] = err;
1697
1698 fEma[0][0] = err*err;
1699 fCorr[0][0] = 1.0;
1700 //--------------------
1701
1702 //--------------------------------------------------
1703 // reset the errors of the points in the histogram
1704 for (Int_t i=1; i<=nbins; i++)
1705 {
1706 fHist->SetBinError(i, saveError[i-1]);
1707 }
1708
1709
1710 return kTRUE;
1711 }
1712
1713
1714 //=========== start loop for reducing the degree ==================
1715 // of the polynomial
1716 while (1)
1717 {
1718 //--------------------------------------------------
1719 // prepare fit of a polynomial : (a0 + a1*x + a2*x**2 + a3*x**3 + ...)
1720
1721 TString funcname = "Poly";
1722 Double_t xmin = 0.0;
1723 Double_t xmax = 90.0;
1724
1725 TString formula = "[0]";
1726 TString bra1 = "+[";
1727 TString bra2 = "]";
1728 TString xpower = "*x";
1729 TString newpower = "*x";
1730 for (Int_t i=1; i<=fDegree; i++)
1731 {
1732 formula += bra1;
1733 formula += i;
1734 formula += bra2;
1735 formula += xpower;
1736
1737 xpower += newpower;
1738 }
1739
1740 //*fLog << "FitPolynomial : formula = " << formula << endl;
1741
1742 fPoly = new TF1(funcname, formula, xmin, xmax);
1743 TList *funclist = fHist->GetListOfFunctions();
1744 funclist->Add(fPoly);
1745
1746 //------------------------
1747 // attention : the dimensions must agree with those in CallMinuit()
1748 const UInt_t npar = fDegree+1;
1749
1750 TString parname[npar];
1751 TArrayD vinit(npar);
1752 TArrayD step(npar);
1753 TArrayD limlo(npar);
1754 TArrayD limup(npar);
1755 TArrayI fix(npar);
1756
1757 vinit[0] = mean;
1758 vinit[2] = a2init;
1759
1760 for (UInt_t j=0; j<npar; j++)
1761 {
1762 parname[j] = "p";
1763 parname[j] += j+1;
1764
1765 step[j] = vinit[j] != 0.0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
1766 }
1767
1768 // limit the first coefficient of the polynomial to positive values
1769 // because the background must not be negative
1770 limup[0] = fHist->GetEntries();
1771
1772 // use the subsequernt loop if you want to apply the
1773 // constraint : uneven derivatives (at alpha=0) = zero
1774 for (UInt_t j=1; j<npar; j+=2)
1775 {
1776 vinit[j] = 0;
1777 step[j] = 0;
1778 fix[j] = 1;
1779 }
1780
1781 //*fLog << "FitPolynomial : before CallMinuit()" << endl;
1782
1783 MMinuitInterface inter;
1784 const Bool_t rc = inter.CallMinuit(fcnpoly, parname, vinit, step,
1785 limlo, limup, fix, fHist, "Migrad",
1786 kFALSE);
1787
1788 //*fLog << "FitPolynomial : after CallMinuit()" << endl;
1789
1790 if (rc != 0)
1791 {
1792 // *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit failed"
1793 // << endl;
1794 // return kFALSE;
1795 }
1796
1797
1798 //-------------------
1799 // get status of minimization
1800 Double_t fmin = 0;
1801 Double_t fedm = 0;
1802 Double_t errdef = 0;
1803 Int_t npari = 0;
1804 Int_t nparx = 0;
1805
1806 if (gMinuit)
1807 gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fIstat);
1808
1809 *fLog << "MHFindSignificanceONOFF::FitPolynomial; fmin, fedm, errdef, npari, nparx, fIstat = "
1810 << fmin << ", " << fedm << ", " << errdef << ", " << npari
1811 << ", " << nparx << ", " << fIstat << endl;
1812
1813
1814 //-------------------
1815 // store the results
1816
1817 Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
1818 fChisq = fmin;
1819 fNdf = fMbins - nparfree;
1820 fProb = TMath::Prob(fChisq, fNdf);
1821
1822
1823 // get fitted parameter values and errors
1824 fValues.Set(npar);
1825 fErrors.Set(npar);
1826
1827 for (Int_t j=0; j<=fDegree; j++)
1828 {
1829 Double_t val, err;
1830 if (gMinuit)
1831 gMinuit->GetParameter(j, val, err);
1832
1833 fPoly->SetParameter(j, val);
1834 fPoly->SetParError(j, err);
1835
1836 fValues[j] = val;
1837 fErrors[j] = err;
1838 }
1839
1840
1841 //--------------------------------------------------
1842 // if the highest coefficient (j0) of the polynomial
1843 // is consistent with zero reduce the degree of the polynomial
1844
1845 Int_t j0 = 0;
1846 for (Int_t j=fDegree; j>1; j--)
1847 {
1848 // ignore fixed parameters
1849 if (fErrors[j] == 0)
1850 continue;
1851
1852 // this is the highest coefficient
1853 j0 = j;
1854 break;
1855 }
1856
1857 if (!fReduceDegree || j0==0 || TMath::Abs(fValues[j0]) > fErrors[j0])
1858 break;
1859
1860 // reduce the degree of the polynomial
1861 *fLog << "MHFindSignificanceONOFF::FitPolynomial; reduce the degree of the polynomial from "
1862 << fDegree << " to " << (j0-2) << endl;
1863 fDegree = j0 - 2;
1864
1865 funclist->Remove(fPoly);
1866 //if (fPoly)
1867 delete fPoly;
1868 fPoly = NULL;
1869
1870 // delete the Minuit object in order to have independent starting
1871 // conditions for the next minimization
1872 //if (gMinuit)
1873 delete gMinuit;
1874 gMinuit = NULL;
1875 }
1876 //=========== end of loop for reducing the degree ==================
1877 // of the polynomial
1878
1879
1880 //--------------------------------------------------
1881 // get the error matrix of the fitted parameters
1882
1883
1884 if (fIstat >= 1)
1885 {
1886 // error matrix was calculated
1887 if (gMinuit)
1888 gMinuit->mnemat(&fEmat[0][0], fNdim);
1889
1890 // copy covariance matrix into a matrix which includes also the fixed
1891 // parameters
1892 TString name;
1893 Double_t bnd1, bnd2, val, err;
1894 Int_t jvarbl;
1895 Int_t kvarbl;
1896 for (Int_t j=0; j<=fDegree; j++)
1897 {
1898 if (gMinuit)
1899 gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
1900
1901 for (Int_t k=0; k<=fDegree; k++)
1902 {
1903 if (gMinuit)
1904 gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
1905
1906 fEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fEmat[jvarbl-1][kvarbl-1];
1907 }
1908 }
1909 }
1910 else
1911 {
1912 // error matrix was not calculated, construct it
1913 *fLog << "MHFindSignificanceONOFF::FitPolynomial; error matrix not defined"
1914 << endl;
1915 for (Int_t j=0; j<=fDegree; j++)
1916 {
1917 for (Int_t k=0; k<=fDegree; k++)
1918 fEma[j][k] = 0;
1919
1920 fEma[j][j] = fErrors[j]*fErrors[j];
1921 }
1922 }
1923
1924
1925 //--------------------------------------------------
1926 // calculate correlation matrix
1927 for (Int_t j=0; j<=fDegree; j++)
1928 for (Int_t k=0; k<=fDegree; k++)
1929 {
1930 const Double_t sq = fEma[j][j]*fEma[k][k];
1931 fCorr[j][k] = sq==0 ? 0 : fEma[j][k] / TMath::Sqrt(fEma[j][j]*fEma[k][k]);
1932 }
1933
1934
1935 //--------------------------------------------------
1936 // reset the errors of the points in the histogram
1937 for (Int_t i=1; i<=nbins; i++)
1938 fHist->SetBinError(i, saveError[i-1]);
1939
1940
1941 return kTRUE;
1942}
1943
1944// --------------------------------------------------------------------------
1945//
1946// ReBinHistogramOFF
1947//
1948// rebin the histogram 'fHistOrig' by grouping 'nrebin' bins together
1949// put the result into the histogram 'fHist'
1950// the rebinning is made such that 'x0' remains a lower bound of a bin
1951//
1952
1953Bool_t MHFindSignificanceONOFF::RebinHistogramOFF(Double_t x0, Int_t nrebin)
1954{
1955 //-----------------------------------------
1956 // search bin i0 which has x0 as lower edge
1957
1958 Int_t i0 = -1;
1959 Int_t nbold = fHistOrigOFF->GetNbinsX();
1960 for (Int_t i=1; i<=nbold; i++)
1961 {
1962 if (TMath::Abs(fHistOrigOFF->GetBinLowEdge(i) - x0) < 1.e-4 )
1963 {
1964 i0 = i;
1965 break;
1966 }
1967 }
1968
1969 if (i0 == -1)
1970 {
1971 i0 = 1;
1972 *fLog << "MHFindsignificanceOFF::RebinOFF; no bin found with " << x0
1973 << " as lower edge, start rebinning with bin 1" << endl;
1974 }
1975
1976 Int_t istart = i0 - nrebin * ( (i0-1)/nrebin );
1977
1978 //-----------------------------------------
1979 // get new bin edges
1980
1981 const Int_t nbnew = (nbold-istart+1) / nrebin;
1982 const Double_t xmin = fHistOrigOFF->GetBinLowEdge(istart);
1983 const Double_t xmax = xmin + (Double_t)nbnew * nrebin * fHistOrigOFF->GetBinWidth(1);
1984 fHistOFF->SetBins(nbnew, xmin, xmax);
1985
1986 *fLog << "MHFindSignificanceONOFF::ReBinOFF; x0, i0, nbold, nbnew, xmin, xmax = "
1987 << x0 << ", " << i0 << ", " << nbold << ", " << nbnew << ", "
1988 << xmin << ", " << xmax << endl;
1989
1990 //-----------------------------------------
1991 // get new bin entries
1992
1993 for (Int_t i=1; i<=nbnew; i++)
1994 {
1995 Int_t j = nrebin*(i-1) + istart;
1996
1997 Double_t content = 0;
1998 Double_t error2 = 0;
1999 for (Int_t k=0; k<nrebin; k++)
2000 {
2001 content += fHistOrigOFF->GetBinContent(j+k);
2002 error2 += fHistOrigOFF->GetBinError(j+k) * fHistOrigOFF->GetBinError(j+k);
2003 }
2004 fHistOFF->SetBinContent(i, content);
2005 fHistOFF->SetBinError (i, sqrt(error2));
2006 }
2007 fHistOFF->SetEntries( fHistOrigOFF->GetEntries() );
2008
2009 return kTRUE;
2010}
2011
2012
2013
2014
2015// --------------------------------------------------------------------------
2016
2017
2018
2019//
2020// ReBinHistogram
2021//
2022// rebin the histogram 'fHistOrig' by grouping 'nrebin' bins together
2023// put the result into the histogram 'fHist'
2024// the rebinning is made such that 'x0' remains a lower bound of a bin
2025//
2026
2027Bool_t MHFindSignificanceONOFF::RebinHistogram(Double_t x0, Int_t nrebin)
2028{
2029 //-----------------------------------------
2030 // search bin i0 which has x0 as lower edge
2031
2032 Int_t i0 = -1;
2033 Int_t nbold = fHistOrig->GetNbinsX();
2034 for (Int_t i=1; i<=nbold; i++)
2035 {
2036 if (TMath::Abs(fHistOrig->GetBinLowEdge(i) - x0) < 1.e-4 )
2037 {
2038 i0 = i;
2039 break;
2040 }
2041 }
2042
2043 if (i0 == -1)
2044 {
2045 i0 = 1;
2046 *fLog << "MHFindsignificance::Rebin; no bin found with " << x0
2047 << " as lower edge, start rebinning with bin 1" << endl;
2048 }
2049
2050 Int_t istart = i0 - nrebin * ( (i0-1)/nrebin );
2051
2052 //-----------------------------------------
2053 // get new bin edges
2054
2055 const Int_t nbnew = (nbold-istart+1) / nrebin;
2056 const Double_t xmin = fHistOrig->GetBinLowEdge(istart);
2057 const Double_t xmax = xmin + (Double_t)nbnew * nrebin * fHistOrig->GetBinWidth(1);
2058 fHist->SetBins(nbnew, xmin, xmax);
2059
2060 *fLog << "MHFindSignificanceONOFF::ReBin; x0, i0, nbold, nbnew, xmin, xmax = "
2061 << x0 << ", " << i0 << ", " << nbold << ", " << nbnew << ", "
2062 << xmin << ", " << xmax << endl;
2063
2064 //-----------------------------------------
2065 // get new bin entries
2066
2067 for (Int_t i=1; i<=nbnew; i++)
2068 {
2069 Int_t j = nrebin*(i-1) + istart;
2070
2071 Double_t content = 0;
2072 Double_t error2 = 0;
2073 for (Int_t k=0; k<nrebin; k++)
2074 {
2075 content += fHistOrig->GetBinContent(j+k);
2076 error2 += fHistOrig->GetBinError(j+k) * fHistOrig->GetBinError(j+k);
2077 }
2078 fHist->SetBinContent(i, content);
2079 fHist->SetBinError (i, sqrt(error2));
2080 }
2081 fHist->SetEntries( fHistOrig->GetEntries() );
2082
2083 return kTRUE;
2084}
2085
2086// --------------------------------------------------------------------------
2087//
2088// FitGaussPoly
2089//
2090// fits a (Gauss + polynomial function) to the alpha distribution 'fhist'
2091//
2092//
2093Bool_t MHFindSignificanceONOFF::FitGaussPoly()
2094{
2095 *fLog << "Entry FitGaussPoly" << endl;
2096
2097 //--------------------------------------------------
2098 // check the histogram :
2099 // - calculate initial values of the parameters
2100 // - check for bins with zero entries
2101 // - set minimum errors
2102 // - save the original errors
2103 // - set errors huge outside the fit range
2104 // (in 'fcnpoly' points with huge errors will be ignored)
2105
2106
2107 Double_t dummy = 1.e20;
2108
2109 fGNzero = 0;
2110 fGMbins = 0;
2111
2112 //------------------------------------------
2113 // if a constant background has been assumed (due to low statistics)
2114 // fit only in the signal region
2115 if ( !fConstantBackg )
2116 {
2117 fAlphalow = 0.0;
2118 fAlphahig = fAlphamax;
2119 }
2120 else
2121 {
2122 fAlphalow = 0.0;
2123 fAlphahig = 2.0*fAlphasig>25.0 ? 25.0 : 2.0*fAlphasig;
2124 }
2125 //------------------------------------------
2126
2127
2128 fAlphalo = 10000.0;
2129 fAlphahi = -10000.0;
2130
2131
2132 Int_t nbins = fHist->GetNbinsX();
2133 TArrayD saveError(nbins);
2134
2135 for (Int_t i=1; i<=nbins; i++)
2136 {
2137 saveError[i-1] = fHist->GetBinError(i);
2138
2139 // bin should be completely contained in the fit range
2140 // (fAlphalow, fAlphahig)
2141 Double_t xlo = fHist->GetBinLowEdge(i);
2142 Double_t xup = fHist->GetBinLowEdge(i+1);
2143
2144 if ( xlo >= fAlphalow-fEps && xlo <= fAlphahig+fEps &&
2145 xup >= fAlphalow-fEps && xup <= fAlphahig+fEps )
2146 {
2147 fGMbins++;
2148
2149 if ( xlo < fAlphalo )
2150 fAlphalo = xlo;
2151
2152 if ( xup > fAlphahi )
2153 fAlphahi = xup;
2154
2155 Double_t content = fHist->GetBinContent(i);
2156
2157
2158 // count bins with zero entry
2159 if (content <= 0.0)
2160 fGNzero++;
2161
2162 // set minimum error
2163 if (content < 9.0)
2164 fHist->SetBinError(i, 3.0);
2165
2166 //*fLog << "Take : i, content, error = " << i << ", "
2167 // << fHist->GetBinContent(i) << ", "
2168 // << fHist->GetBinError(i) << endl;
2169
2170 continue;
2171 }
2172 // bin is not completely contained in the fit range : set error huge
2173
2174 fHist->SetBinError(i, dummy);
2175
2176 //*fLog << "Omit : i, content, error = " << i << ", "
2177 // << fHist->GetBinContent(i) << ", " << fHist->GetBinError(i)
2178 // << endl;
2179
2180 }
2181
2182
2183 // if a bin has no entries don't fit
2184 if (fGNzero > 0)
2185 {
2186 *fLog << "MHFindSignificanceONOFF::FitGaussPoly; out of " << fGMbins
2187 << " bins there are " << fGNzero
2188 << " bins with zero entry" << endl;
2189
2190 fGPoly = NULL;
2191 return kFALSE;
2192 }
2193
2194
2195 //--------------------------------------------------
2196 // prepare fit of a (polynomial+Gauss) :
2197 // (a0 + a1*x + a2*x**2 + a3*x**3 + ...) + A*exp( -0.5*((x-x0)/sigma)**2 )
2198
2199 TString funcname = "PolyGauss";
2200 Double_t xmin = 0.0;
2201 Double_t xmax = 90.0;
2202
2203 TString xpower = "*x";
2204 TString newpower = "*x";
2205
2206 TString formulaBackg = "[0]";
2207 for (Int_t i=1; i<=fDegree; i++)
2208 formulaBackg += Form("+[%d]*x^%d", i, i);
2209
2210 const TString formulaGauss =
2211 Form("[%d]/[%d]*exp(-0.5*((x-[%d])/[%d])^2)",
2212 fDegree+1, fDegree+3, fDegree+2, fDegree+3);
2213
2214 TString formula = formulaBackg;
2215 formula += "+";
2216 formula += formulaGauss;
2217
2218 *fLog << "FitGaussPoly : formulaBackg = " << formulaBackg << endl;
2219 *fLog << "FitGaussPoly : formulaGauss = " << formulaGauss << endl;
2220 *fLog << "FitGaussPoly : formula = " << formula << endl;
2221
2222 fGPoly = new TF1(funcname, formula, xmin, xmax);
2223 TList *funclist = fHist->GetListOfFunctions();
2224 funclist->Add(fGPoly);
2225
2226 fGBackg = new TF1("Backg", formulaBackg, xmin, xmax);
2227 //funclist->Add(fGBackg);
2228
2229 //------------------------
2230 // attention : the dimensions must agree with those in CallMinuit()
2231 Int_t npar = fDegree+1 + 3;
2232
2233 TString parname[npar];
2234 TArrayD vinit(npar);
2235 TArrayD step(npar);
2236 TArrayD limlo(npar);
2237 TArrayD limup(npar);
2238 TArrayI fix(npar);
2239
2240
2241 // take as initial values for the polynomial
2242 // the result from the polynomial fit
2243 for (Int_t j=0; j<=fDegree; j++)
2244 vinit[j] = fPoly->GetParameter(j);
2245
2246 Double_t sigma = 8;
2247 vinit[fDegree+1] = 2.0 * fNexONOFF * fHist->GetBinWidth(1) / TMath::Sqrt(TMath::Pi()*2);
2248 vinit[fDegree+2] = 0;
2249 vinit[fDegree+3] = sigma;
2250
2251 *fLog << "FitGaussPoly : starting value for Gauss-amplitude = "
2252 << vinit[fDegree+1] << endl;
2253
2254 for (Int_t j=0; j<npar; j++)
2255 {
2256 parname[j] = "p";
2257 parname[j] += j+1;
2258
2259 step[j] = vinit[j]!=0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
2260 }
2261
2262 // limit the first coefficient of the polynomial to positive values
2263 // because the background must not be negative
2264 limup[0] = fHist->GetEntries()*10;
2265
2266 // limit the sigma of the Gauss function
2267 limup[fDegree+3] = 20;
2268
2269
2270 // use the subsequernt loop if you want to apply the
2271 // constraint : uneven derivatives (at alpha=0) = zero
2272 for (Int_t j=1; j<=fDegree; j+=2)
2273 {
2274 vinit[j] = 0;
2275 step[j] = 0;
2276 fix[j] = 1;
2277 }
2278
2279 // fix position of Gauss function
2280 vinit[fDegree+2] = 0;
2281 step[fDegree+2] = 0;
2282 fix[fDegree+2] = 1;
2283
2284 // if a constant background has been assumed (due to low statistics)
2285 // fix the background
2286 if (fConstantBackg)
2287 {
2288 step[0] = 0;
2289 fix[0] = 1;
2290 }
2291
2292 MMinuitInterface inter;
2293 const Bool_t rc = inter.CallMinuit(fcnpolygauss, parname, vinit, step,
2294 limlo, limup, fix, fHist, "Migrad",
2295 kFALSE);
2296
2297 if (rc != 0)
2298 {
2299 // *fLog << "MHFindSignificanceONOFF::FitGaussPoly; (polynomial+Gauss) fit failed"
2300 // << endl;
2301 // return kFALSE;
2302 }
2303
2304
2305 //-------------------
2306 // get status of the minimization
2307 Double_t fmin;
2308 Double_t fedm;
2309 Double_t errdef;
2310 Int_t npari;
2311 Int_t nparx;
2312
2313 if (gMinuit)
2314 gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fGIstat);
2315
2316 *fLog << "MHFindSignificanceONOFF::FitGaussPoly; fmin, fedm, errdef, npari, nparx, fGIstat = "
2317 << fmin << ", " << fedm << ", " << errdef << ", " << npari
2318 << ", " << nparx << ", " << fGIstat << endl;
2319
2320
2321 //-------------------
2322 // store the results
2323
2324 Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
2325 fGChisq = fmin;
2326 fGNdf = fGMbins - nparfree;
2327 fGProb = TMath::Prob(fGChisq, fGNdf);
2328
2329
2330 // get fitted parameter values and errors
2331 fGValues.Set(npar);
2332 fGErrors.Set(npar);
2333
2334 for (Int_t j=0; j<npar; j++)
2335 {
2336 Double_t val, err;
2337 if (gMinuit)
2338 gMinuit->GetParameter(j, val, err);
2339
2340 fGPoly->SetParameter(j, val);
2341 fGPoly->SetParError(j, err);
2342
2343 fGValues[j] = val;
2344 fGErrors[j] = err;
2345
2346 if (j <=fDegree)
2347 {
2348 fGBackg->SetParameter(j, val);
2349 fGBackg->SetParError(j, err);
2350 }
2351 }
2352
2353 fSigmaGauss = fGValues[fDegree+3];
2354 fdSigmaGauss = fGErrors[fDegree+3];
2355 // fitted total number of excess events
2356 fNexGauss = fGValues[fDegree+1] * TMath::Sqrt(TMath::Pi()*2) /
2357 (fHist->GetBinWidth(1)*2 );
2358 fdNexGauss = fNexGauss * fGErrors[fDegree+1]/fGValues[fDegree+1];
2359
2360 //--------------------------------------------------
2361 // get the error matrix of the fitted parameters
2362
2363
2364 if (fGIstat >= 1)
2365 {
2366 // error matrix was calculated
2367 if (gMinuit)
2368 gMinuit->mnemat(&fGEmat[0][0], fGNdim);
2369
2370 // copy covariance matrix into a matrix which includes also the fixed
2371 // parameters
2372 TString name;
2373 Double_t bnd1, bnd2, val, err;
2374 Int_t jvarbl;
2375 Int_t kvarbl;
2376 for (Int_t j=0; j<npar; j++)
2377 {
2378 if (gMinuit)
2379 gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
2380
2381 for (Int_t k=0; k<npar; k++)
2382 {
2383 if (gMinuit)
2384 gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
2385
2386 fGEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fGEmat[jvarbl-1][kvarbl-1];
2387 }
2388 }
2389 }
2390 else
2391 {
2392 // error matrix was not calculated, construct it
2393 *fLog << "MHFindSignificanceONOFF::FitPolynomial; error matrix not defined"
2394 << endl;
2395 for (Int_t j=0; j<npar; j++)
2396 {
2397 for (Int_t k=0; k<npar; k++)
2398 fGEma[j][k] = 0;
2399
2400 fGEma[j][j] = fGErrors[j]*fGErrors[j];
2401 }
2402 }
2403
2404
2405 //--------------------------------------------------
2406 // calculate correlation matrix
2407 for (Int_t j=0; j<npar; j++)
2408 {
2409 for (Int_t k=0; k<npar; k++)
2410 {
2411 const Double_t sq = fGEma[j][j]*fGEma[k][k];
2412 fGCorr[j][k] = sq==0 ? 0 : fGEma[j][k] / sqrt( fGEma[j][j]*fGEma[k][k] );
2413 }
2414 }
2415
2416
2417 //--------------------------------------------------
2418 // reset the errors of the points in the histogram
2419 for (Int_t i=1; i<=nbins; i++)
2420 fHist->SetBinError(i, saveError[i-1]);
2421
2422 return kTRUE;
2423
2424}
2425
2426
2427
2428// --------------------------------------------------------------------------
2429//
2430// DetExcessONOFF
2431//
2432// using the result of the polynomial fit (fValuesOFF), DetExcessONOFF determines
2433//
2434// - the total number of events in the signal region (fNon)
2435// - the number of backgound events (fitted) in the signal region (fNoffSigFitted)
2436// - the number of OFF (normalized evetns) in the alpha OFF distribution (fNoffSig)
2437// - the number of excess events (fNexONOFF)
2438// - the number of excess events using the fitted OFF (fNexONOFFFitted)
2439// - the effective number of background events (fNoff), and fGamma :
2440// fNoffSigFitted = fGamma * fNoff; fdNoffSigFitted = fGamma * sqrt(fNoff);
2441//
2442// It assumed that the polynomial is defined as
2443// a0 + a1*x + a2*x**2 + a3*x**3 + ..
2444//
2445// and that the alpha distribution has the range 0 < alpha < 90 degrees
2446//
2447
2448Bool_t MHFindSignificanceONOFF::DetExcessONOFF()
2449{
2450 //*fLog << "MHFindSignificanceONOFF::DetExcessONOFF;" << endl;
2451
2452//--------------------------------------------
2453 // calculate the total number of events (fNon) in the signal region
2454
2455 fNon = 0.0;
2456 fdNon = 0.0;
2457
2458 Double_t alphaup = -1000.0;
2459 Double_t binwidth = fHist->GetBinWidth(1);
2460
2461 Int_t nbins = fHist->GetNbinsX();
2462 for (Int_t i=1; i<=nbins; i++)
2463 {
2464 Double_t xlo = fHist->GetBinLowEdge(i);
2465 Double_t xup = fHist->GetBinLowEdge(i+1);
2466
2467 // bin must be completely contained in the signal region
2468 if ( xlo <= (fAlphasig+fEps) && xup <= (fAlphasig+fEps) )
2469 {
2470 Double_t width = fabs(xup-xlo);
2471 if (fabs(width-binwidth) > fEps)
2472 {
2473 *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alpha plot has variable binning, which is not allowed"
2474 << endl;
2475 return kFALSE;
2476 }
2477
2478 if (xup > alphaup)
2479 alphaup = xup;
2480
2481 fNon += fHist->GetBinContent(i);
2482 fdNon += fHist->GetBinError(i) * fHist->GetBinError(i);
2483 }
2484 }
2485 fdNon = sqrt(fdNon);
2486
2487 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF:" << endl
2488 << "fNon = " << fNon << " ; fdNon = " << fdNon << endl;
2489
2490
2491 // tmp
2492 //if (fNon == 0)
2493 //{
2494 // cout << "ERROR... WITH COUNTED OF ON EVETNS: fAlphasig, fEps = "
2495// << fAlphasig << ", " << fEps << endl;
2496 // fHist-> DrawCopy();
2497 // gPad -> SaveAs("HistON.ps");
2498 //}
2499 // endtmp
2500
2501
2502 // the actual signal range is :
2503 if (alphaup == -1000.0)
2504 return kFALSE;
2505
2506 fAlphasi = alphaup;
2507
2508
2509 //*fLog << "fAlphasi, fNon, fdNon, binwidth, fDegree = " << fAlphasi << ", "
2510 // << fNon << ", " << fdNon << ", " << binwidth << ", "
2511 // << fDegree << endl;
2512
2513
2514 // Calculate the number of OFF events in the signal region
2515 // ___________________________________________________
2516
2517
2518 fNoffSig = 0.0;
2519 fdNoffSig = 0.0;
2520
2521 Double_t alphaupOFF = -1000.0;
2522 Double_t binwidthOFF = fHistOFF->GetBinWidth(1);
2523
2524 Int_t nbinsOFF = fHistOFF->GetNbinsX();
2525
2526 for (Int_t i=1; i<=nbinsOFF; i++)
2527 {
2528 Double_t xlo = fHistOFF->GetBinLowEdge(i);
2529 Double_t xup = fHistOFF->GetBinLowEdge(i+1);
2530
2531 // bin must be completely contained in the signal region
2532 if ( xlo <= (fAlphasig+fEps) && xup <= (fAlphasig+fEps) )
2533 {
2534 Double_t width = fabs(xup-xlo);
2535 if (fabs(width-binwidthOFF) > fEps)
2536 {
2537 *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alphaOFF plot has variable binning, which is not allowed"
2538 << endl;
2539 return kFALSE;
2540 }
2541
2542 if (xup > alphaupOFF)
2543 alphaup = xup;
2544
2545 fNoffSig += fHistOFF->GetBinContent(i);
2546 fdNoffSig += fHistOFF->GetBinError(i) * fHistOFF->GetBinError(i);
2547 }
2548 }
2549 fdNoffSig = sqrt(fdNoffSig);
2550
2551 // tmp
2552 //if (fNoffSig == 0)
2553 //{
2554 // cout << "ERROR... WITH COUNTED OF OFF EVETNS: fAlphasig, fEps = "
2555// << fAlphasig << ", " << fEps << endl;
2556 // fHistOFF-> DrawCopy();
2557 // gPad -> SaveAs("HistOFF.ps");
2558 //}
2559 //endtmp
2560
2561 // the actual signal range is :
2562 if (alphaup == -1000.0)
2563 return kFALSE;
2564
2565 fAlphasiOFF = alphaup;
2566
2567 if (fabs(fAlphasiOFF - fAlphasi) > fEps)
2568 {
2569 *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; fAlphasiOFF ("
2570 << fAlphasiOFF << ") is not equal to fAlphasi ("
2571 << fAlphasi << "), this is something that should not happen"
2572 << endl;
2573
2574 //return kFALSE; It might happen in pathological cases (very few OFF)
2575 // and I want to see the results, anyhow
2576 }
2577
2578
2579
2580 // Calculate the number of OFF events in the total OFF region
2581 // defined by fAlphaminOFF and fAlphamaxOFF
2582 // ___________________________________________________
2583
2584 fNoffTot = 0.0;
2585 fdNoffTot = 0.0;
2586
2587
2588
2589 for (Int_t i=1; i<=nbinsOFF; i++)
2590 {
2591 Double_t xlo = fHistOFF->GetBinLowEdge(i);
2592 Double_t xup = fHistOFF->GetBinLowEdge(i+1);
2593
2594 // bin must be completely contained in the signal region
2595 if ( xlo >= (fAlphaminOFF-fEps) && xup <= (fAlphamaxOFF+fEps) )
2596 {
2597 Double_t width = fabs(xup-xlo);
2598 if (fabs(width-binwidthOFF) > fEps)
2599 {
2600 *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alphaOFF plot has variable binning, which is not allowed"
2601 << endl;
2602 return kFALSE;
2603 }
2604
2605 fNoffTot += fHistOFF->GetBinContent(i);
2606 fdNoffTot += fHistOFF->GetBinError(i) * fHistOFF->GetBinError(i);
2607 }
2608 }
2609 fdNoffTot = sqrt(fdNoffTot);
2610
2611
2612
2613
2614
2615
2616
2617
2618
2619 //--------------------------------------------
2620 // calculate the number of OFF fitted events (fNoffSigFitted) in the signal region
2621 // and its error (fdNoffSigFitted)
2622
2623
2624
2625 if (fUseFittedQuantities)
2626 {
2627 //--------------------------------------------
2628 // calculate the number of OFF fitted events (fNoffSigFitted) in the signal region
2629 // and its error (fdNoffSigFitted)
2630
2631 Double_t fac = 1.0/binwidthOFF;
2632
2633 fNoffSigFitted = 0.0;
2634 Double_t altothejplus1 = fAlphasi; // Limit for signal found for ON data is used.
2635 for (Int_t j=0; j<=fDegreeOFF; j++)
2636 {
2637 fNoffSigFitted += fValuesOFF[j] * altothejplus1 / ((Double_t)(j+1));
2638 altothejplus1 *= fAlphasi;
2639 }
2640 fNoffSigFitted *= fac;
2641
2642 // derivative of fNoffSigFitted
2643 Double_t facj;
2644 Double_t fack;
2645
2646 Double_t sum = 0.0;
2647 altothejplus1 = fAlphasi;
2648 for (Int_t j=0; j<=fDegreeOFF; j++)
2649 {
2650 facj = altothejplus1 / ((Double_t)(j+1));
2651
2652 Double_t altothekplus1 = fAlphasi;
2653 for (Int_t k=0; k<=fDegreeOFF; k++)
2654 {
2655 fack = altothekplus1 / ((Double_t)(k+1));
2656 sum += facj * fack * fEmaOFF[j][k];
2657 altothekplus1 *= fAlphasi;
2658 }
2659 altothejplus1 *= fAlphasi;
2660 }
2661 sum *= fac*fac;
2662
2663 if (sum < 0.0)
2664 {
2665 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; error squared is negative"
2666 << endl;
2667 return kFALSE;
2668 }
2669
2670 fdNoffSigFitted = sqrt(sum);
2671
2672
2673 // We can now compare fNoffSig with fNoffSigFitted (and their errors)
2674 // NUmbers should agree within 10 % (errors within 20%)
2675
2676 if (fabs(fNoffSig - fNoffSigFitted) > 0.1 * fNoffSigFitted)
2677 {
2678 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; number of OFF events and Fitted number of OFF events in signal region do not agree (within 10 %)" << endl;
2679
2680 *fLog << "fNoffSig = " << fNoffSig << " ; fNoffSigFitted = " << fNoffSigFitted << endl;
2681
2682
2683 // return kFALSE; NOt yet...
2684 }
2685
2686 /*
2687 if (fabs(fdNoffSig - fdNoffSigFitted) > 0.2 * fdNoffSigFitted)
2688 {
2689 *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 %)"
2690 << endl;
2691
2692
2693 *fLog << "fdNoffSig = " << fdNoffSig << " ; fdNoffSigFitted = " << fdNoffSigFitted << endl;
2694
2695
2696 //return kFALSE; NOt yet...
2697 }
2698
2699 */
2700
2701
2702
2703 // Calculate the number of OFF events in the whole fit region (fAlphaminOFF-fAlphamaxOFF)
2704 // ___________________________________________________
2705
2706
2707
2708 fNoffTotFitted = 0.0;
2709
2710 altothejplus1 = fAlphamaxOFF; // Limit for OFF data fit
2711 for (Int_t j=0; j<=fDegreeOFF; j++)
2712 {
2713 fNoffTotFitted += fValuesOFF[j] * altothejplus1 / ((Double_t)(j+1));
2714 altothejplus1 *= fAlphamaxOFF;
2715 }
2716 fNoffTotFitted *= fac;
2717
2718 // derivative of fdNoffTotFitted
2719
2720
2721 sum = 0.0;
2722 altothejplus1 = fAlphamaxOFF;
2723 for (Int_t j=0; j<=fDegreeOFF; j++)
2724 {
2725 facj = altothejplus1 / ((Double_t)(j+1));
2726
2727 Double_t altothekplus1 = fAlphamaxOFF;
2728 for (Int_t k=0; k<=fDegreeOFF; k++)
2729 {
2730 fack = altothekplus1 / ((Double_t)(k+1));
2731
2732 sum += facj * fack * fEmaOFF[j][k];
2733 altothekplus1 *= fAlphamaxOFF;
2734 }
2735 altothejplus1 *= fAlphamaxOFF;
2736 }
2737 sum *= fac*fac;
2738
2739 if (sum < 0.0)
2740 {
2741 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; error squared is negative"
2742 << endl;
2743 return kFALSE;
2744 }
2745
2746 fdNoffTotFitted = sqrt(sum);
2747
2748
2749
2750 }
2751
2752else
2753 {
2754 fNoffSigFitted = 0.0;
2755 fdNoffSigFitted = 0.0;
2756 fNoffTotFitted = 0.0;
2757 fdNoffTotFitted = 0.0;
2758 }
2759
2760
2761
2762 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF; INFO ABOOUT COMPUTED OFF EVENTS..." << endl
2763 << "fNoffSig = " << fNoffSig << "; fdNoffSig = " << fdNoffSig << endl
2764 << "fNoffSigFitted = " << fNoffSigFitted
2765 << "; fdNoffSigFitted = " << fdNoffSigFitted << endl;
2766
2767
2768
2769
2770
2771 //--------------------------------------------
2772 // calculate the number of excess events in the signal region
2773
2774 fNexONOFF = fNon - fNoffSig*fNormFactor;
2775 fNexONOFFFitted = fNon - fNoffSigFitted*fNormFactor;
2776
2777
2778 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF;" << endl
2779 << "fNexONOFF (= fNon - fNoffSig*fNormFactor) = " << fNexONOFF << endl
2780 << "fNexONOFFFitted (= fNon - fNoffSigFitted*fNormFactor) = " << fNexONOFFFitted << endl;
2781
2782
2783 //--------------------------------------------
2784 // calculate the effective number of background events (fNoff) , and fGamma :
2785 // fNbg = fGamma * fNoff = fNormFactor* fNoffSigFitted;
2786 // dfNbg = fGamma * sqrt(fNoff) = fNormFactor * fdNoffSigFitted;
2787
2788 if (fNoffSigFitted < 0.0)
2789 {
2790 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF; number of fitted OFF events in signal region is negative, fNoffSigFitted, fdNoffSigFitted = "
2791 << fNoffSigFitted << ", " << fdNoffSigFitted << endl;
2792
2793 fGamma = 1.0;
2794 fNoff = 0.0;
2795 return kFALSE;
2796 }
2797
2798 if (fNoffSigFitted > 0.0)
2799 {
2800 fGamma = fNormFactor * fdNoffSigFitted*fdNoffSigFitted/fNoffSigFitted;
2801 fNoff = fNormFactor * fNoffSigFitted/fGamma;
2802 }
2803 else
2804 {
2805 fGamma = 1.0;
2806 fNoff = 0.0;
2807 }
2808
2809 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF: " << endl
2810 << "fGamma = " << fGamma << " ; fNoff = " << fNoff << endl;
2811
2812
2813 return kTRUE;
2814}
2815
2816
2817
2818
2819
2820// --------------------------------------------------------------------------
2821//
2822// SigmaLiMa
2823//
2824// calculates the significance according to Li & Ma
2825// ApJ 272 (1983) 317; formula 17
2826//
2827Bool_t MHFindSignificanceONOFF::SigmaLiMa(Double_t non, Double_t noff,
2828 Double_t gamma, Double_t *siglima)
2829{
2830 if (gamma <= 0.0 || non <= 0.0 || noff <= 0.0)
2831 {
2832 *siglima = 0.0;
2833 return kFALSE;
2834 }
2835
2836 Double_t help1 = non * log( (1.0+gamma)*non / (gamma*(non+noff)) );
2837 Double_t help2 = noff * log( (1.0+gamma)*noff / ( non+noff ) );
2838 *siglima = sqrt( 2.0 * (help1+help2) );
2839
2840 Double_t nex = non - gamma*noff;
2841 if (nex < 0.0)
2842 *siglima = - *siglima;
2843
2844 //*fLog << "MHFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = "
2845 // << non << ", " << noff << ", " << gamma << ", " << *siglima << endl;
2846
2847 return kTRUE;
2848}
2849
2850
2851
2852// calculates the significance according to Li & Ma
2853// ApJ 272 (1983) 317; formula 5
2854//
2855Bool_t MHFindSignificanceONOFF::SigmaLiMaForm5(Double_t non, Double_t noff,
2856 Double_t gamma, Double_t *siglima)
2857{
2858 if (gamma <= 0.0 || non <= 0.0 || noff <= 0.0)
2859 {
2860 *siglima = 0.0;
2861 return kFALSE;
2862 }
2863
2864 Double_t nex = non - (gamma*noff);
2865 Double_t tmp = non + (gamma*gamma)*noff;
2866 tmp = TMath::Sqrt(tmp);
2867
2868 *siglima = nex/tmp;
2869
2870 if (nex < 0.0)
2871 *siglima = - *siglima;
2872
2873 //*fLog << "MHFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = "
2874 // << non << ", " << noff << ", " << gamma << ", " << *siglima << endl;
2875
2876 return kTRUE;
2877}
2878
2879
2880
2881
2882// --------------------------------------------------------------------------
2883//
2884
2885// Following function computes a clone of fHistOFF and normalizes
2886// contents, errors and fPolyOFF (if exists) with the fNormFactor.
2887// This normalized OFF hist will be used when plotting OFF data
2888// together with ON data.
2889
2890Bool_t MHFindSignificanceONOFF::ComputeHistOFFNormalized()
2891{
2892
2893
2894 if (!fHist)
2895 {
2896 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fHist does not exist, normalization of HistOFF can not be performed properly..."
2897 << endl;
2898 return kFALSE;
2899
2900 }
2901
2902
2903
2904 if (!fHistOFF)
2905 {
2906 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fHistOFF does not exist, hence can not be normalized"
2907 << endl;
2908 return kFALSE;
2909
2910 }
2911
2912
2913 if (fNormFactor <= 0)
2914 {
2915 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fNormFactor is ZERO or NEGATIVE, it might have not been defined yet..."
2916 << endl;
2917 return kFALSE;
2918
2919 }
2920
2921
2922 Double_t BinWidthAlphaON = fHist -> GetBinWidth(1);
2923 Double_t BinWidthAlphaOFF = fHistOFF -> GetBinWidth(1);
2924 Double_t BinWidthRatioONOFF = BinWidthAlphaON/BinWidthAlphaOFF;
2925
2926
2927 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; INFO about alpha ON, OFF histo bins"
2928 << endl
2929 // << "fHist bin width = " << BinWidthAlphaON
2930 << " fHistOFF bin width = " << BinWidthAlphaOFF << endl;
2931
2932
2933
2934 TString fHistOFFNormalizedName = fHistOFF -> GetName();
2935 fHistOFFNormalizedName += (" (Normalized)");
2936 // fHistOFFNormalized = (TH1*) fHistOFF -> Clone();
2937 // fHistOFFNormalized -> SetNameTitle(fHistOFFNormalizedName, fHistOFFNormalizedName);
2938
2939
2940 Int_t nbinsOFFNormalized = 0;
2941 Int_t nbinsOFF = 0;
2942 Double_t xlow = 0.0;
2943 Double_t xup = 0.0;
2944 Double_t content = 0.0;
2945 Double_t error = 0.0;
2946 Double_t BinCenter = 0.0;
2947
2948
2949 // Bins for normalized OFF histo will be the ones of ON histo
2950
2951
2952 nbinsOFF = fHistOFF -> GetNbinsX();
2953 nbinsOFFNormalized = nbinsOFF;
2954 xlow = fHistOFF -> GetBinLowEdge(1);
2955 xup = fHistOFF -> GetBinLowEdge(nbinsOFFNormalized);
2956 xup = xup + BinWidthAlphaON;
2957
2958 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; Limits for fHistOFFNormalized: "
2959 << "nbins, xlow, xup: " << nbinsOFFNormalized << ", "
2960 << xlow << ", " << xup << endl;
2961
2962 fHistOFFNormalized = new TH1F (fHistOFFNormalizedName, fHistOFFNormalizedName,
2963 nbinsOFFNormalized, xlow, xup);
2964
2965
2966 // fHistOFFNormalized is filled with data from fHistOFF,
2967 // taken into account the possible bin difference
2968
2969
2970 for (Int_t i = 1; i <= nbinsOFF; i++)
2971 {
2972 BinCenter = fHistOFF -> GetBinCenter(i);
2973 fHistOFFNormalized -> Fill (BinCenter, fHistOFF -> GetBinContent(i));
2974 fHistOFFNormalized -> SetBinError(i, fHistOFF -> GetBinError(i));
2975 }
2976
2977
2978
2979
2980 for (Int_t i = 1; i <= nbinsOFFNormalized; i++)
2981 {
2982 content = fNormFactor * fHistOFFNormalized -> GetBinContent(i);
2983 error = fNormFactor * fHistOFFNormalized -> GetBinError(i);
2984
2985 fHistOFFNormalized -> SetBinContent (i, content);
2986 fHistOFFNormalized -> SetBinError (i, error);
2987 }
2988
2989
2990 // Number of entries is obtained from histOFF.
2991 // and set to histOFFNoramlized; otherwise, the number
2992 // of entries in histOFFNoramlized would be "nbins"
2993
2994 Double_t entries = fNormFactor * (fHistOFF -> GetEntries());
2995 fHistOFFNormalized -> SetEntries(entries);
2996
2997
2998
2999
3000 // If polynomial fit has been performed for fHistOFF,
3001 // it is defined a new polyfunction for fHistOFFNormalized,
3002 // which will be the polyfunction of fHistOFF normalized.
3003 // Function will be added to the function list of fHistOFFNormalized
3004
3005
3006
3007 if (fPolyOFF == NULL)
3008 {
3009 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fPolyOFF does not exist..."
3010 << endl;
3011 }
3012
3013 if (fPolyOFF)
3014 {
3015 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fPolyOFF exists... will be also normalized and included in the list of functions of fHistOFFNormalized"
3016 << endl;
3017
3018
3019
3020
3021 // Normalization of the function using fNormFactor and
3022 // BinWidthON/BinWidthOFF relation of alpha ON/OFF histograms
3023 // This makes possible to plot it together with ON alpha histo
3024
3025 TString FunctionName("PolyOFFNormalized");
3026
3027 Double_t xmin = fAlphaminOFF;
3028 Double_t xmax = fAlphamaxOFF;
3029
3030 TString formula = "[0]";
3031 TString bra1 = "+[";
3032 TString bra2 = "]";
3033 TString xpower = "*x";
3034 TString newpower = "*x";
3035 for (Int_t i=1; i<=fDegreeOFF; i++)
3036 {
3037 formula += bra1;
3038 formula += i;
3039 formula += bra2;
3040 formula += xpower;
3041
3042 xpower += newpower;
3043 }
3044
3045 *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; formula = " << formula << endl;
3046
3047
3048
3049 fPolyOFFNormalized = new TF1 (FunctionName, formula, xmin, xmax);
3050
3051
3052 Double_t Parameter = 0.0;
3053 Double_t ParameterError = 0.0;
3054
3055 *fLog << " MHFindsignificanceONOFF::ComputeHistOFFNormalized; Fit parameters info: " << endl;
3056 for (Int_t i = 0; i <= fDegreeOFF; i++)
3057 {
3058 Parameter = fNormFactor * BinWidthRatioONOFF * fValuesOFF[i];
3059 ParameterError = fNormFactor * BinWidthRatioONOFF * fErrorsOFF[i];
3060
3061 fPolyOFFNormalized -> SetParameter(i, Parameter);
3062 fPolyOFFNormalized -> SetParError(i,ParameterError);
3063
3064
3065
3066 // Parameters are shown :
3067
3068 *fLog << " fValuesOFF[" << i<< "] = " << fValuesOFF[i]
3069 << " ; Parameter for fPolyOFFNormalized = " << Parameter << endl;
3070
3071 }
3072
3073
3074 TList *funclist = fHistOFFNormalized->GetListOfFunctions();
3075
3076 // temporal...
3077 //*fLog << "INFO concerning list of functions of fHistOFFNormalized :" << endl
3078 // << "List before adding OFF Normal., after adding it and after removing fPolyOFF..."
3079 // << endl;
3080
3081 //funclist-> Print();
3082 funclist-> Add(fPolyOFFNormalized);
3083
3084 //funclist-> Print();
3085
3086
3087
3088 }
3089
3090 return kTRUE;
3091
3092}
3093
3094// --------------------------------------------------------------------------
3095//
3096
3097Bool_t MHFindSignificanceONOFF::DrawHistOFF()
3098{
3099 if (fHistOFF == NULL )
3100 {
3101 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; fHistOFF = NULL" << endl;
3102 return kFALSE;
3103 }
3104
3105// fPsFilename -> NewPage();
3106
3107 // PLOT DISABLE
3108 /*
3109 TCanvas* CanvasHistOFF = new TCanvas(fHistOFF->GetName(), fHistOFF->GetName(), 600, 600);
3110
3111 //gStyle->SetOptFit(1011);
3112
3113 gROOT->SetSelectedPad(NULL);
3114 gStyle->SetPadLeftMargin(0.1);
3115 gStyle -> SetOptStat(1);
3116
3117 CanvasHistOFF->cd();
3118
3119
3120 if (fHistOFF)
3121 {
3122 fHistOFF->DrawCopy();
3123 }
3124
3125 // TF1 *fpolyOFF = fHistOFF->GetFunction("PolyOFF");
3126 if (fPolyOFF == NULL)
3127 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; fpolyOFF = NULL" << endl;
3128
3129 if (fPolyOFF)
3130 {
3131 // 2, 1 is red and solid
3132 fPolyOFF->SetLineColor(2);
3133 fPolyOFF->SetLineStyle(1);
3134 fPolyOFF->SetLineWidth(2);
3135 fPolyOFF->DrawCopy("same");
3136 }
3137
3138 CanvasHistOFF -> Update();
3139
3140
3141 */
3142
3143 return kTRUE;
3144}
3145
3146
3147// --------------------------------------------------------------------------
3148//
3149
3150Bool_t MHFindSignificanceONOFF::DrawHistOFFNormalized()
3151{
3152 if (fHistOFFNormalized == NULL )
3153 {
3154 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; fHistOFFNormalized = NULL" << endl;
3155 return kFALSE;
3156 }
3157
3158 // fPsFilename -> NewPage();
3159
3160 // PLOT DISABLE TO PERFORM GRID ANALYSIS
3161 /*
3162 TCanvas* CanvasHistOFFNormalized = new TCanvas(fHistOFFNormalized->GetName(),
3163 fHistOFFNormalized->GetName(), 600, 600);
3164
3165 //gStyle->SetOptFit(1011);
3166
3167 // gROOT->SetSelectedPad(NULL);
3168 gStyle->SetPadLeftMargin(0.1);
3169 gStyle -> SetOptStat(1);
3170
3171 CanvasHistOFFNormalized->cd();
3172
3173
3174 if (fHistOFFNormalized)
3175 {
3176 fHistOFFNormalized->DrawCopy();
3177 }
3178
3179 // TF1 *fpolyOFFNormalized = fHistOFFNormalized->GetFunction("PolyOFFNormalized");
3180 if (fPolyOFFNormalized == NULL)
3181 *fLog << "MHFindSignificanceONOFF::DrawHistOFF; fPolyOFFNormalized = NULL" << endl;
3182
3183 if (fPolyOFFNormalized)
3184 {
3185 // 2, 1 is red and solid
3186 fPolyOFFNormalized->SetLineColor(2);
3187 fPolyOFFNormalized->SetLineStyle(1);
3188 fPolyOFFNormalized->SetLineWidth(2);
3189 fPolyOFFNormalized->DrawCopy("same");
3190 }
3191
3192 CanvasHistOFFNormalized -> Update();
3193
3194
3195 */
3196
3197 return kTRUE;
3198
3199}
3200
3201
3202// --------------------------------------------------------------------------
3203//
3204
3205Bool_t MHFindSignificanceONOFF::DrawFit(const Option_t *opt)
3206{
3207 if (fHistOFFNormalized == NULL || fHist == NULL)
3208 *fLog << "MHFindSignificanceONOFF::DrawFit; fHistOFFNormalized = NULL or fHist == NULL" << endl;
3209
3210 //fPsFilename -> NewPage();
3211
3212
3213 // PLOT DISABLE TO PERFORM GRID ANALYSIS
3214 // I DO SAVE PS FILE, BUT CANVAS IS DELETED AFTERWARDS
3215
3216 fCanvas = new TCanvas("Alpha", "Alpha plot", 600, 600);
3217 fCanvas -> SetFillColor(10);
3218
3219
3220
3221
3222 //gStyle->SetOptFit(1011);
3223
3224 gROOT->SetSelectedPad(NULL);
3225 gStyle -> SetFrameFillColor(10);
3226 gStyle->SetPadLeftMargin(0.15);
3227 gStyle -> SetOptStat(1);
3228
3229 fCanvas->cd();
3230
3231 if (fHist)
3232 {
3233 fHist -> SetTitle("Alpha Plot");
3234 fHist-> SetTitleOffset(1.5, "Y");
3235 fHist-> DrawCopy();
3236
3237 }
3238
3239
3240 if (fHistOFFNormalized)
3241 {
3242 TF1 *fpoly = fHistOFFNormalized->GetFunction("PolyOFFNormalized");
3243 if (fpoly == NULL)
3244 *fLog << "MHFindSignificanceONOFF::DrawFit; fPolyOFFNormalized = NULL" << endl;
3245
3246 if (fpoly)
3247 {
3248 // 2, 1 is red and solid
3249 fpoly->SetLineColor(2);
3250 fpoly->SetLineStyle(1);
3251 fpoly->SetLineWidth(2);
3252 fpoly->DrawCopy("same");
3253 }
3254 }
3255
3256 if (fFitGauss)
3257 {
3258 TF1 *fpolygauss = fHist->GetFunction("PolyGauss");
3259 if (fpolygauss == NULL)
3260 *fLog << "MHFindSignificanceONOFF::DrawFit; fpolygauss = NULL" << endl;
3261
3262 if (fpolygauss)
3263 {
3264 // 4, 1 is blue and solid
3265 fpolygauss->SetLineColor(4);
3266 fpolygauss->SetLineStyle(1);
3267 fpolygauss->SetLineWidth(4);
3268 fpolygauss->DrawCopy("same");
3269 }
3270
3271 TF1 *fbackg = fHist->GetFunction("Backg");
3272 if (fbackg == NULL)
3273 *fLog << "MHFindSignificanceONOFF::DrawFit; fbackg = NULL" << endl;
3274
3275 if (fbackg)
3276 {
3277 // 10, 1 is white and solid
3278 fbackg->SetLineColor(10);
3279 fbackg->SetLineStyle(1);
3280 fbackg->SetLineWidth(4);
3281 // fbackg->DrawCopy("same"); I do not want to draw it... already too many things.
3282 }
3283 }
3284
3285
3286 //-------------------------------
3287 // print results onto the figure
3288
3289
3290
3291 TPaveText *pt = new TPaveText(0.30, 0.35, 0.70, 0.90, "NDC");
3292 char tx[100];
3293
3294 sprintf(tx, "Results of polynomial fit to OFF (order %2d) :", fDegreeOFF);
3295 TText *t1 = pt->AddText(tx);
3296 t1->SetTextSize(0.03);
3297 t1->SetTextColor(2);
3298
3299 sprintf(tx, " (%6.2f< |alpha| <%6.2f [\\circ])", fAlphaminOFF, fAlphamaxOFF);
3300 pt->AddText(tx);
3301
3302 sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f",
3303 fChisqOFF, fNdfOFF, fProbOFF);
3304 pt->AddText(tx);
3305
3306 sprintf(tx, " OFF events (fit)= %8.1f #pm %8.1f",
3307 fNoffTotFitted, fdNoffTotFitted);
3308 pt->AddText(tx);
3309
3310 sprintf(tx, " OFF events (meas) = %8.1f #pm %8.1f", fNoffTot, fdNoffTot);
3311 pt->AddText(tx);
3312
3313 sprintf(tx, " OFF Normalization Factor (= Non/Noff) = %4.4f", fNormFactor);
3314 pt->AddText(tx);
3315
3316
3317
3318
3319 //sprintf(tx, " ");
3320 //pt->AddText(tx);
3321
3322 //--------------
3323 sprintf(tx, "Results for |alpha|< %6.2f [\\circ] :", fAlphasi);
3324 TText *t6 = pt->AddText(tx);
3325 t6->SetTextSize(0.03);
3326 t6->SetTextColor(8);
3327
3328 sprintf(tx, " Non = %8.1f #pm %8.1f", fNon, fdNon);
3329 pt->AddText(tx);
3330
3331
3332 if(fUseFittedQuantities)
3333 {
3334 //// **************************************************
3335 ///// PRINT INFORMATION ABOUT FITTED QUANTITIES /////////
3336
3337
3338
3339 Double_t NoffFitNormalized = fNoffSigFitted * fNormFactor;
3340 Double_t ErrorNoffFitNormalized = fdNoffSigFitted * fNormFactor;
3341 Double_t SignificanceUsed = GetSignificance();
3342
3343 sprintf(tx, " Noff Fitted (Normalized) = %8.1f #pm %8.1f",
3344 NoffFitNormalized, ErrorNoffFitNormalized);
3345 pt->AddText(tx);
3346
3347
3348 sprintf(tx, " Nex (ON - OFF Fitted) = %8.1f #pm %8.1f",
3349 fNexONOFFFitted, fdNexONOFFFitted);
3350 pt->AddText(tx);
3351
3352
3353 sprintf(tx, " Gamma = %4.4f, Effective Noff (i.e. fNoff) = %6.1f",
3354 fGamma, fNoff);
3355 pt->AddText(tx);
3356
3357
3358 Double_t ratio = fNoffSigFitted>0.0 ? fNexONOFFFitted/(fNoffSigFitted*fNormFactor) : 0.0;
3359 sprintf(tx, " Significance = %6.2f, Nex/(Nbg*NormFactor) = %6.2f",
3360 SignificanceUsed, ratio);
3361 pt->AddText(tx);
3362
3363
3364 }
3365
3366 else
3367 {
3368 //// **************************************************
3369 ///// PRINT INFORMATION ABOUT MEASURED QUANTITIES /////////
3370
3371
3372 Double_t NoffNormalized = fNoffSig * fNormFactor;
3373 Double_t ErrorNoffNormalized = fdNoffSig * fNormFactor;
3374 Double_t SignificanceUsed = GetSignificance();
3375
3376 sprintf(tx, " Noff measured (Normalized) = %8.1f #pm %8.1f",
3377 NoffNormalized, ErrorNoffNormalized);
3378 pt->AddText(tx);
3379
3380 sprintf(tx, " Nex (ON - OFF measured) = %8.1f #pm %8.1f",
3381 fNexONOFF, fdNexONOFF);
3382 pt->AddText(tx);
3383
3384 Double_t ratio = fNoffSig>0.0 ? fNexONOFF/(fNoffSig*fNormFactor) : 0.0;
3385 sprintf(tx, " Significance = %6.2f, Nex/(Nbg*NormFactor) = %6.2f",
3386 SignificanceUsed, ratio);
3387 pt->AddText(tx);
3388
3389 }
3390
3391 /*
3392 // Temporally I will also show ALL SIGMALIMA COMPUTED.
3393
3394 sprintf(tx,
3395 " fSigLiMa1 = %6.2f, fSigLiMa2 = %6.2f, fSigLiMa3 = %6.2f",
3396 fSigLiMa,fSigLiMa2, fSigLiMa3);
3397 pt->AddText(tx);
3398 */
3399
3400
3401 //--------------
3402 if (fFitGauss)
3403 {
3404 sprintf(tx, "Results of (polynomial+Gauss) fit :");
3405 TText *t7 = pt->AddText(tx);
3406 t7->SetTextSize(0.03);
3407 t7->SetTextColor(4);
3408
3409 sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f",
3410 fGChisq, fGNdf, fGProb);
3411 pt->AddText(tx);
3412
3413 sprintf(tx, " Sigma of Gauss = %8.1f #pm %8.1f [\\circ]",
3414 fSigmaGauss, fdSigmaGauss);
3415 pt->AddText(tx);
3416
3417 sprintf(tx, " total no.of excess events = %8.1f #pm %8.1f",
3418 fNexGauss, fdNexGauss);
3419 pt->AddText(tx);
3420 }
3421 //--------------
3422
3423 pt->SetFillStyle(0);
3424 pt->SetBorderSize(0);
3425 pt->SetTextAlign(12);
3426
3427
3428 if(fPrintResultsOntoAlphaPlot)
3429 {
3430 pt->Draw();
3431 }
3432 fCanvas->Modified();
3433 fCanvas->Update();
3434
3435 // fPsFilename -> NewPage();
3436
3437
3438
3439 if (fSavePlots)
3440 {
3441 // ********************************************
3442 // TMP solution while the TPostScript thing is not working.
3443 // PsFileName for storing these histograms is derived
3444 // from fPsFilenameString.
3445
3446
3447 cout << "Alpha plot with ON-OFF data will be saved in PostScript file " ;
3448
3449
3450 if (!fPsFilenameString.IsNull())
3451 {
3452 TString filename = (fPsFilenameString);
3453 // Train or Test Sample is specified outside
3454 // class MHFindSignificanceONOFF, and included in
3455 // fPsFilenameString
3456
3457 filename += ("AlphaPlotAfterSupercuts.ps");
3458 cout << filename << endl;
3459 fCanvas -> SaveAs(filename);
3460 }
3461
3462 // END OF TEMPORAL SOLUTION
3463 // ********************************************
3464
3465 }
3466
3467 // Canvvas deleted to allow for GRID analysis
3468
3469 delete fCanvas;
3470
3471
3472 return kTRUE;
3473}
3474
3475
3476// --------------------------------------------------------------------------
3477//
3478// Print the results of the polynomial fit to the alpha OFF distribution
3479//
3480//
3481void MHFindSignificanceONOFF::PrintPolyOFF(Option_t *o)
3482{
3483 *fLog << "---------------------------" << endl;
3484 *fLog << "MHFindSignificanceONOFF::PrintPolyOFF :" << endl;
3485
3486 *fLog << "fAlphaminOFF, fAlphamaxOFF, fDegreeOFF "
3487 << fAlphaminOFF << ", " << fAlphamaxOFF << ", " << fDegreeOFF << endl;
3488
3489 *fLog << "fMbinsOFF, fNzeroOFF, fIstatOFF = " << fMbinsOFF << ", "
3490 << fNzeroOFF << ", " << fIstatOFF << endl;
3491
3492 *fLog << "fChisqOFF, fNdfOFF, fProbOFF = " << fChisqOFF << ", "
3493 << fNdfOFF << ", " << fProbOFF << endl;
3494
3495 *fLog << "fNon; fNoffSigFitted, fdNoffSigFitted; fNoffSig, fdNoffSig = "
3496 << fNon << "; " << fNoffSigFitted << ", " << fdNoffSigFitted
3497 << "; " << fNoffSig << ", " << fdNoffSig << endl;
3498
3499 Double_t sigtoback = fNoffSigFitted >0.0 ? fNexONOFFFitted/(fNoffSigFitted*fNormFactor) : 0.0;
3500
3501
3502 *fLog << "fNexONOFFFitted, fdNexONOFFFitted, fGamma, fNoff, fSigLiMa, sigtoback = "
3503 << fNexONOFFFitted << ", " << fdNexONOFFFitted << ", "
3504 << fGamma << ", " << fNoff
3505 << ", " << fSigLiMa << ", "
3506 << sigtoback << endl;
3507
3508
3509 Double_t sigtoback2 = fNoffSig >0.0 ? fNexONOFF/(fNoffSig*fNormFactor) : 0.0;
3510
3511 *fLog << "fNexONOFF, fdNexONOFF, fNormFactor, fNoffSig, fSigLiMa2, sigtoback2 = "
3512 << fNexONOFF << ", " << fdNexONOFF << ", "
3513 << fNormFactor << ", " << fNoffSig
3514 << ", " << fSigLiMa2 << ", "
3515 << sigtoback2 << endl;
3516
3517 *fLog << "---------------------------" << endl;
3518}
3519
3520// --------------------------------------------------------------------------
3521//
3522// Print the results of the polynomial fit to the alpha distribution
3523//
3524//
3525void MHFindSignificanceONOFF::PrintPoly(Option_t *o)
3526{
3527 *fLog << "---------------------------" << endl;
3528 *fLog << "MHFindSignificanceONOFF::PrintPoly :" << endl;
3529
3530 *fLog << "fAlphami, fAlphama, fDegree, fAlphasi = "
3531 << fAlphami << ", " << fAlphama << ", " << fDegree << ", "
3532 << fAlphasi << endl;
3533
3534 *fLog << "fMbins, fNzero, fIstat = " << fMbins << ", "
3535 << fNzero << ", " << fIstat << endl;
3536
3537 *fLog << "fChisq, fNdf, fProb = " << fChisq << ", "
3538 << fNdf << ", " << fProb << endl;
3539
3540 *fLog << "fNon, fNbg, fdNbg, fNbgtot, fNbgtotFitted, fdNbgtotFitted = "
3541 << fNon << ", " << fNbg << ", " << fdNbg << ", " << fNbgtot
3542 << ", " << fNbgtotFitted << ", " << fdNbgtotFitted << endl;
3543
3544 Double_t sigtoback = fNbg>0.0 ? fNex/fNbg : 0.0;
3545 *fLog << "fNex, fdNex, fGamma, fNoff, fSigLiMa, sigtoback = "
3546 << fNex << ", " << fdNex << ", " << fGamma << ", " << fNoff
3547 << ", " << fSigLiMa << ", " << sigtoback << endl;
3548
3549 //------------------------------------
3550 // get errors
3551
3552 /*
3553 Double_t eplus;
3554 Double_t eminus;
3555 Double_t eparab;
3556 Double_t gcc;
3557 Double_t errdiag;
3558
3559
3560 if ( !fConstantBackg )
3561 {
3562 *fLog << "parameter value error eplus eminus eparab errdiag gcc"
3563 << endl;
3564
3565 for (Int_t j=0; j<=fDegree; j++)
3566 {
3567 if (gMinuit)
3568 gMinuit->mnerrs(j, eplus, eminus, eparab, gcc);
3569 errdiag = sqrt(fEma[j][j]);
3570 *fLog << j << " " << fValues[j] << " " << fErrors[j] << " "
3571 << eplus << " " << eminus << " " << eparab << " "
3572 << errdiag << " " << gcc << endl;
3573 }
3574 }
3575 else
3576 {
3577 *fLog << "parameter value error errdiag "
3578 << endl;
3579
3580 for (Int_t j=0; j<=fDegree; j++)
3581 {
3582 errdiag = sqrt(fEma[j][j]);
3583 *fLog << j << " " << fValues[j] << " " << fErrors[j] << " "
3584 << errdiag << endl;
3585 }
3586 }
3587 */
3588
3589 //----------------------------------------
3590 /*
3591 *fLog << "Covariance matrix :" << endl;
3592 for (Int_t j=0; j<=fDegree; j++)
3593 {
3594 *fLog << "j = " << j << " : ";
3595 for (Int_t k=0; k<=fDegree; k++)
3596 {
3597 *fLog << fEma[j][k] << " ";
3598 }
3599 *fLog << endl;
3600 }
3601
3602 *fLog << "Correlation matrix :" << endl;
3603 for (Int_t j=0; j<=fDegree; j++)
3604 {
3605 *fLog << "j = " << j << " : ";
3606 for (Int_t k=0; k<=fDegree; k++)
3607 {
3608 *fLog << fCorr[j][k] << " ";
3609 }
3610 *fLog << endl;
3611 }
3612 */
3613
3614 *fLog << "---------------------------" << endl;
3615}
3616
3617// --------------------------------------------------------------------------
3618//
3619// Print the results of the (polynomial+Gauss) fit to the alpha distribution
3620//
3621//
3622void MHFindSignificanceONOFF::PrintPolyGauss(Option_t *o)
3623{
3624 *fLog << "---------------------------" << endl;
3625 *fLog << "MHFindSignificanceONOFF::PrintPolyGauss :" << endl;
3626
3627 *fLog << "fAlphalo, fAlphahi = "
3628 << fAlphalo << ", " << fAlphahi << endl;
3629
3630 *fLog << "fGMbins, fGNzero, fGIstat = " << fGMbins << ", "
3631 << fGNzero << ", " << fGIstat << endl;
3632
3633 *fLog << "fGChisq, fGNdf, fGProb = " << fGChisq << ", "
3634 << fGNdf << ", " << fGProb << endl;
3635
3636
3637 //------------------------------------
3638 // get errors
3639
3640 Double_t eplus;
3641 Double_t eminus;
3642 Double_t eparab;
3643 Double_t gcc;
3644 Double_t errdiag;
3645
3646 *fLog << "parameter value error eplus eminus eparab errdiag gcc"
3647 << endl;
3648 for (Int_t j=0; j<=(fDegree+3); j++)
3649 {
3650 if (gMinuit)
3651 gMinuit->mnerrs(j, eplus, eminus, eparab, gcc);
3652 errdiag = sqrt(fGEma[j][j]);
3653 *fLog << j << " " << fGValues[j] << " " << fGErrors[j] << " "
3654 << eplus << " " << eminus << " " << eparab << " "
3655 << errdiag << " " << gcc << endl;
3656 }
3657
3658
3659 *fLog << "Covariance matrix :" << endl;
3660 for (Int_t j=0; j<=(fDegree+3); j++)
3661 {
3662 *fLog << "j = " << j << " : ";
3663 for (Int_t k=0; k<=(fDegree+3); k++)
3664 {
3665 *fLog << fGEma[j][k] << " ";
3666 }
3667 *fLog << endl;
3668 }
3669
3670 *fLog << "Correlation matrix :" << endl;
3671 for (Int_t j=0; j<=(fDegree+3); j++)
3672 {
3673 *fLog << "j = " << j << " : ";
3674 for (Int_t k=0; k<=(fDegree+3); k++)
3675 {
3676 *fLog << fGCorr[j][k] << " ";
3677 }
3678 *fLog << endl;
3679 }
3680
3681 *fLog << "---------------------------" << endl;
3682}
3683
3684//============================================================================
3685
3686
3687
3688
3689
3690
3691
3692
3693
3694
3695
3696
3697
3698
3699
3700
3701
3702
3703
3704
3705
3706
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718
3719
3720
Note: See TracBrowser for help on using the repository browser.