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

Last change on this file since 6724 was 6354, checked in by mazin, 20 years ago
*** empty log message ***
File size: 117.6 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): David Paneque 12/2003 <mailto:dpaneque@mppmu.mpg.de>
19!
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27// //
28// MFindSupercutsONOFFThetaLoop //
29// //
30// Class for optimizing the parameters of the supercuts //
31// Using ON and OFF data runing over data which are subdivided in
32// several Theta ranges. //
33// //
34// //
35/////////////////////////////////////////////////////////////////////////////
36
37
38#include "MFindSupercutsONOFFThetaLoop.h"
39
40
41#include <math.h> // fabs
42
43#include <TFile.h>
44#include <TArrayD.h>
45#include <TMinuit.h>
46#include <TCanvas.h>
47#include <TStopwatch.h>
48#include <TVirtualFitter.h>
49
50#include <iostream.h>
51#include <fstream.h>
52#include <stdlib.h>
53#include <stdio.h>
54#include <string.h>
55
56#include "MFindSupercutsONOFF.h"
57#include "MBinning.h"
58#include "MContinue.h"
59#include "MSupercuts.h"
60#include "MSupercutsCalcONOFF.h"
61#include "MDataElement.h"
62#include "MDataMember.h"
63
64
65#include <TPostScript.h>
66
67#include "MEvtLoop.h"
68#include "MFCT1SelFinal.h"
69#include "MF.h"
70#include "MFEventSelector.h"
71#include "MFEventSelector2.h"
72#include "MFillH.h"
73//#include "MGeomCamCT1Daniel.h"
74//#include "MGeomCamCT1.h"
75#include "MGeomCamMagic.h" // only for magic
76#include "MFRandomSplit.h"
77#include "MH3.h"
78#include "MHCT1Supercuts.h"
79#include "MHFindSignificance.h" // To be removed at some point...
80#include "MHFindSignificanceONOFF.h"
81#include "MHMatrix.h"
82#include "MHOnSubtraction.h"
83#include "MDataValue.h"
84// #include "MDataString.h"
85
86#include "MLog.h"
87#include "MLogManip.h"
88#include "MMatrixLoop.h"
89#include "MMinuitInterface.h"
90#include "MParList.h"
91#include "MProgressBar.h"
92#include "MReadMarsFile.h"
93#include "MReadTree.h"
94#include "MTaskList.h"
95
96
97ClassImp(MFindSupercutsONOFF);
98
99using namespace std;
100
101// --------------------------------------------------------------------------
102//
103// Default constructor.
104//
105MFindSupercutsONOFFThetaLoop::MFindSupercutsONOFFThetaLoop(const char *name, const char *title)
106{
107 fName = name ? name : "MFindSupercutsONOFF";
108 fTitle = title ? title : "Optimizer of the supercuts";
109
110 //---------------------------
111
112 // *********************************
113 // ** NOT TRUE ANY LONGER... ** //
114 // Cuts (low and up) in variable ThetaOrig.fVal
115 // The default is not cut, i.e. all values (0-1) are taken
116 // fThetaOrig.fVal is measured in Radians; thus 1 = 57 degrees.
117
118 // fThetaMin = 0.0;
119 // fThetaMax = 1.0;
120 // *********************************
121
122
123 // For MAGIC new (July 2004) Theta is measured in degrees
124
125 fThetaMin = 0.0;
126 fThetaMax = 90.0;
127
128
129
130 fAlphaSig = 20; // By default, signal is expected in alpha<20 for CT1
131
132 // By default, bkg region is set to alpha range 30-90 for CT1
133 fAlphaBkgMin = 30;
134 fAlphaBkgMax = 90;
135
136 fActualCosThetaBinCenter = 0;
137
138 fNormFactorTrainHist = NULL;
139 //fNormFactorTrainHist -> SetDirectory(0);
140 fNormFactorTestHist = NULL;
141 //fNormFactorTestHist -> SetDirectory(0);
142 fSigmaLiMaTrainHist = NULL;
143 //fSigmaLiMaTrainHist -> SetDirectory(0);
144 fSigmaLiMaTestHist = NULL;
145 //fSigmaLiMaTestHist -> SetDirectory(0);
146 fNexTrainHist = NULL;
147 //fSigmaLiMaTestHist -> SetDirectory(0);
148 fNexTestHist = NULL;
149 //fNexTestHist -> SetDirectory(0);
150
151
152 fSuccessfulThetaBinsHist = NULL;
153
154 /*
155 fNexVSAlphaSigTrain = NULL;
156 fNexVSAlphaSigTest = NULL;
157
158 fSignificanceVSAlphaSigTrain = NULL;
159 fSignificanceVSAlphaSigTest = NULL;
160 */
161
162
163 fNEvtsInTrainMatrixONHist = NULL;
164 fNEvtsInTestMatrixONHist = NULL;
165 fNEvtsInTrainMatrixOFFHist = NULL;
166 fNEvtsInTestMatrixOFFHist = NULL;
167
168
169 fTrainMatrixONFilenameVector = NULL;
170 fTestMatrixONFilenameVector= NULL;
171 fTrainMatrixOFFFilenameVector= NULL;
172 fTestMatrixOFFFilenameVector= NULL;
173
174
175 fOptSCParamFilenameVector = NULL;
176
177 fThetaRangeStringVector = NULL;
178
179 fPsFilename = NULL;
180
181
182
183 fTuneNormFactor = kTRUE; // Norm factors will be corrected using the total amount of OFF events before cuts and the estimated excess events
184
185 // fNormFactorTrain = fNormFactorTrain - Ngammas/EventsInTrainMatrixOFF
186
187 fGammaEfficiency = 0.5; // Fraction of gammas that remain after cuts
188 // Quantity that will have to be determined with MC, yet for the
189 // time being I set it to 0.5
190
191 // boolean variables are set to kFALSE as default
192
193 fReadMatricesFromFile = kFALSE;
194 fOptimizeParameters = kFALSE;
195 fTestParameters = kFALSE;
196
197 // Fraction of Train/Test ON/OFF samples is set to 0.5 as default
198
199 fWhichFractionTrain = 0.5; // number <= 1; specifying fraction of ON Train events
200 fWhichFractionTest = 0.5; // number <= 1; specifying fraction of ON Test events
201
202 fWhichFractionTrainOFF = 0.5; // number <= 1; specifying fraction of OFF Train events
203 fWhichFractionTestOFF = 0.5; // number <= 1; specifying fraction of OFF Test events
204
205
206
207 // use quantities computed from the fits
208 // The variable allows the user to NOT use these quantities when there is not
209 // enough statistics and fit not always is possible
210 fUseFittedQuantities = kTRUE;
211
212
213
214 // Boolean variable that controls wether the optimization of the
215 // parameters (MMinuitInterface::CallMinuit(..) in function FindParams(..))
216 // takes place or not. kTRUE will skip such optimization.
217 // This variable is useful to test the optmized parameters (previously found
218 // and stored in root file) on the TRAIN sample.
219
220 fSkipOptimization = kFALSE;
221
222 // Boolean variable that allows the user to write the initial parameters
223 // into the root file that will be used to store the optimum cuts.
224 // If fUseInitialSCParams = kTRUE , parameters are written.
225 // In this way, the initial SC parameters can be applied on the data (train/test)
226
227 // The initial parameters are ONLY written to the root file if
228 // there is NO SC params optimization, i.e., if variable
229 // fSkipOptimization = kTRUE;
230
231 // The default value is obviously kFALSE.
232
233 fUseInitialSCParams = kTRUE;
234
235
236
237 // fInitSCPar and fInitSCParSteps initialized to zero
238
239 fInitSCPar.Set(0);
240
241 fInitSCParSteps.Set(0);
242
243}
244
245// --------------------------------------------------------------------------
246//
247// Default destructor.
248//
249MFindSupercutsONOFFThetaLoop::~MFindSupercutsONOFFThetaLoop()
250{
251
252 *fLog << "Destructor of MFindSupercutsONOFFThetaLoop is called"
253 << endl;
254
255 fPsFilename = NULL;
256
257 delete fNormFactorTrainHist;
258 delete fNormFactorTestHist;
259 delete fSigmaLiMaTrainHist;
260 delete fSigmaLiMaTestHist;
261 delete fNexTrainHist;
262 delete fNexTestHist;
263
264
265
266 delete fSuccessfulThetaBinsHist;
267
268
269 delete fNEvtsInTrainMatrixONHist;
270 delete fNEvtsInTestMatrixONHist;
271 delete fNEvtsInTrainMatrixOFFHist;
272 delete fNEvtsInTestMatrixOFFHist;
273
274 delete [] fTrainMatrixONFilenameVector;
275 delete [] fTestMatrixONFilenameVector;
276 delete [] fTrainMatrixOFFFilenameVector;
277 delete [] fTestMatrixOFFFilenameVector;
278
279 delete [] fOptSCParamFilenameVector;
280
281 delete [] fThetaRangeStringVector;
282
283
284 *fLog << "Destructor of MFindSupercutsONOFFThetaLoop finished "
285 << "successfully."
286 << endl;
287
288}
289
290
291// --------------------------------------------------------------------------
292//
293// Function that sets the name of the PostScript file where alpha distributions
294// for the different Theta bins will be stored.
295// It also initializes
296
297
298
299void MFindSupercutsONOFFThetaLoop::SetPostScriptFile (TPostScript* PsFile)
300{
301 fPsFilename = PsFile;
302
303 *fLog << "MFindSupercutsONOFF : Results (alpha distributions with excess and significances) will be stored in PostScript file "
304 << fPsFilename -> GetName() << endl;
305
306}
307
308
309Bool_t MFindSupercutsONOFFThetaLoop::SetGammaEfficiency (Double_t gammaeff)
310{
311 if (gammaeff <= 0.0 || gammaeff > 1)
312 {
313 *fLog << " ****** ERROR *********" << endl
314 << "MFindSupercutsONOFFThetaLoop :: SetGammaEfficiency; " << endl
315 << "The requested Gamma efficiency " << gammaeff << " can not occur" << endl
316 << "Member data fGammaEfficiency can not be set properly." << endl
317 << "Default value for fGammaEfficiency, which is "
318 << fGammaEfficiency << " will remain." << endl;
319
320
321 return kFALSE;
322 }
323
324 fGammaEfficiency = gammaeff;
325
326 return kTRUE;
327
328}
329
330
331Bool_t MFindSupercutsONOFFThetaLoop::ReadSCParamsFromAsciiFile(const char* filename,
332 Int_t Nparams)
333{
334 const Int_t NColumns = 2;
335
336 Double_t tmp_vector[NColumns];
337 Int_t counter = 0;
338 Double_t tmpdouble = 0.0;
339
340
341 // File is read and data stored in data_vector
342
343 ifstream infile (filename, ios::in);
344
345 if ( !infile )
346 {
347 cout << "Input file could not be opened" << endl;
348 return kFALSE;
349 // exit (1);
350 }
351
352 counter = 0;
353 while (infile >> tmp_vector [0])
354 {
355 for (Int_t i = 1; i < NColumns; i++)
356 {
357 infile >> tmp_vector[i];
358 }
359 counter++;
360 }
361
362 infile.close();
363
364
365 // Length of TArray vectors is set to the length of the
366 // vectors defined in ascii file.
367 // If such length does not match with the specified
368 // by variable Nparams, the function complains and
369 // aborts.
370
371 if (counter != Nparams)
372 {
373 *fLog << "MFindSupercutsONOFFThetaLoop::ReadSCParamsFromAsciiFile" << endl
374 << "Length of vectors in ascii file " << filename << " is " << counter
375 << " , which does not match with the expected length " << Nparams << endl
376 << " fInitSCPar and fInitSCParSteps are NOT retrieved from ascii file" << endl;
377
378 return kFALSE;
379 }
380
381 fInitSCPar.Set(counter);
382 fInitSCParSteps.Set(counter);
383
384
385 // Vectors are filled
386
387 ifstream infile2 (filename, ios::in);
388
389 if ( !infile2 )
390 {
391 cout << "Input file could not be opened" << endl;
392 return kFALSE;
393 // exit (1);
394 }
395
396 for (Int_t i = 0; i < fInitSCPar.GetSize(); i++)
397 {
398
399 infile2 >> tmpdouble;
400 fInitSCPar[i] = tmpdouble;
401
402 infile2 >> tmpdouble;
403 fInitSCParSteps[i] = tmpdouble;
404
405 }
406
407
408 // Print the vectors read
409
410 *fLog << "***** INITIAL SC PARAMETERS AND STEPS TAKEN FROM ASCII FILE ********* " << endl;
411 for (Int_t i = 0; i < fInitSCPar.GetSize(); i++)
412 {
413 cout << fInitSCPar[i] << setw(20) << fInitSCParSteps[i] << endl;
414 }
415
416 //
417
418 // File close
419 infile2.close();
420
421 return kTRUE;
422}
423
424
425Bool_t MFindSupercutsONOFFThetaLoop::SetInitSCPar (TArrayD &d)
426{
427 *fLog << "MFindSupercutsONOFFThetaLoop; Initial SC parameters are "
428 << "set using function SetInitSCPar (TArrayD &d)" << endl;
429
430 fInitSCPar = d;
431
432 return kTRUE;
433}
434
435
436Bool_t MFindSupercutsONOFFThetaLoop::SetInitSCParSteps (TArrayD &d)
437{
438
439 *fLog << "MFindSupercutsONOFFThetaLoop; Initial SC parameters STEPS are "
440 << "set using function SetInitSCParSteps (TArrayD &d)" << endl;
441
442 fInitSCParSteps = d;
443
444 return kTRUE;
445}
446
447
448
449Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaSig(Double_t alphasig)
450{
451 // check that alpha is within the limits 0-90
452 if (alphasig <= 0 || alphasig > 90)
453 {
454 *fLog << "MFindSupercutsONOFFThetaLoop ::SetAlphaSig; "
455 << "value " << alphasig << " is not within the the "
456 << "logical limits of alpha; 0-90" << endl;
457 return kFALSE;
458 }
459
460
461 fAlphaSig = alphasig;
462
463 return kTRUE;
464}
465
466Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaBkgMin(Double_t alpha)
467{
468 // check that alpha is within the limits 0-90
469 if (alpha <= 0 || alpha >= 90)
470 {
471 *fLog << "MFindSupercutsONOFFThetaLoop::SetAlphaBkgMin; "
472 << "value " << alpha << " is not within the the "
473 << "logical limits of alpha; 0-90" << endl;
474 return kFALSE;
475 }
476
477
478 fAlphaBkgMin = alpha;
479
480 return kTRUE;
481}
482
483
484Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaBkgMax(Double_t alpha)
485{
486 // check that alpha is within the limits 0-90
487 if (alpha <= 0 || alpha > 90.001)
488 {
489 *fLog << "MFindSupercutsONOFFThetaLoop ::SetAlphaBkgMax; "
490 << "value " << alpha << " is not within the the "
491 << "logical limits of alpha; 0-90" << endl;
492 return kFALSE;
493 }
494
495
496 fAlphaBkgMax = alpha;
497
498 return kTRUE;
499}
500
501
502// Function that checks that the values of the member data
503// fAlphaSig, fAlphaBkgMin and fAlphaBkgMax make sense
504// (ie, fAlphaSig < fAlphaBkgMin < fAlphaBkgMax)
505
506Bool_t MFindSupercutsONOFFThetaLoop::CheckAlphaSigBkg()
507{
508
509 if (fAlphaSig > fAlphaBkgMin)
510 {
511 *fLog << "MFindSupercutsONOFFThetaLoop ::CheckAlphaSigBkg(); " << endl
512 << "fAlphaSig > fAlphaBkgMin, which should not occur..." << endl
513 << "fAlphaSig = " << fAlphaSig << ", fAlphaBkgMin = " << fAlphaBkgMin
514 << endl;
515
516 return kFALSE;
517 }
518
519 if (fAlphaBkgMax < fAlphaBkgMin)
520 {
521 *fLog << "MFindSupercutsONOFFThetaLoop::CheckAlphaSigBkg(); " << endl
522 << "fAlphaBkgMin > fAlphaBkgMax, which should not occur..." << endl
523 << "fAlphaBkgMin = " << fAlphaBkgMin << ", fAlphaBkgMax = " << fAlphaBkgMax
524 << endl;
525
526 return kFALSE;
527 }
528
529 return kTRUE;
530
531}
532
533
534
535
536// Function that sets the value of fCosThetaRangeVector and
537// fCosThetaBinCenterVector
538
539Bool_t MFindSupercutsONOFFThetaLoop::SetCosThetaRangeVector(
540 const TArrayD &d)
541{
542 if (d.GetSize() < 2)
543 {
544 *fLog << err
545 << "MFindSupercutsONOFFThetaLoop::SetCosThetaRangeVector: "
546 << endl
547 << "Size of d is smaller than 2... "
548 << "fCosThetaRangeVector can not be defined properly"
549 << endl;
550 return kFALSE;
551 }
552
553 fCosThetaRangeVector.Set(d.GetSize());
554 fCosThetaRangeVector = d;
555
556
557
558
559 // fCosThetaBinCenterVector is defined
560
561 Int_t dim = fCosThetaRangeVector.GetSize() - 1;
562 fCosThetaBinCenterVector.Set(dim);
563
564 for (Int_t i = 0; i < dim; i++)
565 {
566 fCosThetaBinCenterVector[i] =
567 (fCosThetaRangeVector[i+1]+fCosThetaRangeVector[i])/2;
568 }
569
570 return kTRUE;
571}
572
573// Function that sets the range of Theta (fThetaMin, fThetaMax)
574// where data will be used in the iteration
575// thetabin. This theta range is used when geting data from
576// file and filling matrices. If matrices are already built
577// and just read from root files, this theta range is useless.
578
579
580// **** IMPORTANT *********
581// thetabin is an integer value that goes from 1 to
582// fCosThetaBinCenterVector.GetSize()
583// ************************
584
585Bool_t MFindSupercutsONOFFThetaLoop::SetThetaRange (Int_t thetabin)
586{
587 // Check that iteration thetabin is smaller than
588 // components of fCosThetaBinCenterVector
589
590 if (fCosThetaRangeVector.GetSize() < 2 ||
591 fCosThetaBinCenterVector.GetSize() < 1)
592
593 {
594 *fLog << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl
595 << "Dimension of vectors fCosThetaRangeVector or/and "
596 << "fCosThetaBinCenterVector are smaller than 2 and 1 "
597 << "respectively. Range for theta can not be defined." << endl;
598
599 return kFALSE;
600 }
601
602 if (thetabin < 1 ||
603 thetabin > fCosThetaBinCenterVector.GetSize())
604 {
605 *fLog << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl
606 << "thetabin value is " << thetabin
607 << " Is is outside the possible values given by "
608 << "vector fCosThetaRangeVector. "
609 << "Range for theta can not be defined." << endl;
610
611 return kFALSE;
612 }
613
614
615 fThetaMin = fCosThetaRangeVector[thetabin-1];
616 fThetaMax = fCosThetaRangeVector[thetabin];
617
618
619 // Now costheta values have to be converted to Radians
620 // Function TMath::ACos(x) gives angle in Radians
621
622
623 fThetaMin = TMath::ACos(fThetaMin);
624 fThetaMax = TMath::ACos(fThetaMax);
625
626
627 // For new MAGIC (July 2004), angles are given in degrees.
628 // So I convert angles into degrees by multiplying by a the
629 // conversion rad/degree
630
631 fThetaMin = fThetaMin * TMath::RadToDeg();
632 fThetaMax = fThetaMax * TMath::RadToDeg();
633
634
635
636 // If fThetaMin > fThetaMax means that fCosThetaRangeVector
637 // is defined with increasing values of CosTheta, which
638 // means decreasing values of Theta
639
640 // In such case, values of fThetaMin and fThetaMax are swapped
641
642 if (fThetaMin > fThetaMax)
643 {
644 Double_t tmp = fThetaMin;
645 fThetaMin = fThetaMax;
646 fThetaMax = tmp;
647 }
648
649 // Info for user
650 *fLog << "-------------------------------------------------"<< endl
651 << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl
652 << "Theta range for this iteration is set to "
653 << fThetaMin << "-" << fThetaMax << " degrees." << endl
654 << "-------------------------------------------------"<< endl;
655
656 return kTRUE;
657
658
659}
660
661
662Bool_t MFindSupercutsONOFFThetaLoop::SetNormFactorTrainHist()
663{
664 if (fCosThetaRangeVector.GetSize() <= 1)
665 {
666 *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNormFactorTrainHist: "
667 << endl
668 << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
669 << endl;
670 return kFALSE;
671 }
672
673
674 Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
675 Double_t* BinsVector = fCosThetaRangeVector.GetArray();
676
677 fNormFactorTrainHist = new TH1F ("NormFactorTrainHist",
678 "NormFactorTrainHist",
679 HistoNbins,BinsVector );
680
681 return kTRUE;
682
683}
684
685
686Bool_t MFindSupercutsONOFFThetaLoop::SetNormFactorTestHist()
687{
688 if (fCosThetaRangeVector.GetSize() <= 1)
689 {
690 *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNormFactorTestHist: "
691 << endl
692 << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
693 << endl;
694 return kFALSE;
695 }
696
697
698 Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
699 Double_t* BinsVector = fCosThetaRangeVector.GetArray();
700
701 fNormFactorTestHist = new TH1F ("NormFactorTestHist",
702 "NormFactorTestHist",
703 HistoNbins, BinsVector);
704
705 return kTRUE;
706}
707
708
709Bool_t MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTrainHist()
710{
711 if (fCosThetaRangeVector.GetSize() <= 1)
712 {
713 *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTrainHist: "
714 << endl
715 << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
716 << endl;
717 return kFALSE;
718 }
719
720
721 // At some point, definition of bins should be able to
722 // contain variable bins...
723
724 Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
725
726 Double_t* BinsVector = fCosThetaRangeVector.GetArray();
727
728 fSigmaLiMaTrainHist = new TH1F ("SigmaLiMaTrainHist",
729 "SigmaLiMaTrainHist",
730 HistoNbins, BinsVector);
731
732 return kTRUE;
733}
734
735
736
737Bool_t MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTestHist()
738{
739 if (fCosThetaRangeVector.GetSize() <= 1)
740 {
741 *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTestHist: "
742 << endl
743 << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
744 << endl;
745 return kFALSE;
746 }
747
748
749
750 Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
751
752 Double_t* BinsVector = fCosThetaRangeVector.GetArray();
753 fSigmaLiMaTestHist = new TH1F ("SigmaLiMaTestHist",
754 "SigmaLiMaTestHist",
755 HistoNbins, BinsVector);
756
757 return kTRUE;
758
759}
760
761
762Bool_t MFindSupercutsONOFFThetaLoop::SetNexTrainHist()
763{
764 if (fCosThetaRangeVector.GetSize() <= 1)
765 {
766 *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNexTrainHist: "
767 << endl
768 << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
769 << endl;
770 return kFALSE;
771 }
772
773
774 Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
775
776 Double_t* BinsVector = fCosThetaRangeVector.GetArray();
777
778 fNexTrainHist = new TH1F ("NexTrainHist",
779 "NexTrainHist",
780 HistoNbins, BinsVector);
781
782 return kTRUE;
783
784}
785
786
787Bool_t MFindSupercutsONOFFThetaLoop::SetNexTestHist()
788{
789 if (fCosThetaRangeVector.GetSize() <= 1)
790 {
791 *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNexTestHist: "
792 << endl
793 << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
794 << endl;
795 return kFALSE;
796 }
797
798
799 Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
800 Double_t* BinsVector = fCosThetaRangeVector.GetArray();
801
802 fNexTestHist = new TH1F ("NexTestHist",
803 "NexTestHist",
804 HistoNbins, BinsVector);
805
806 return kTRUE;
807
808}
809
810Bool_t MFindSupercutsONOFFThetaLoop::SetSuccessfulThetaBinsHist()
811{
812 if (fCosThetaRangeVector.GetSize() <= 1)
813 {
814 *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSuccessfulThetaBinsHist(): "
815 << endl
816 << "Size of fCosThetaRangeVector is <=1 ... "
817 << "fCosThetaRangeVector must exist and have >= 2 components (defining a single bin)"
818 << endl;
819 return kFALSE;
820 }
821
822
823 Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
824 Double_t* BinsVector = fCosThetaRangeVector.GetArray();
825
826 fSuccessfulThetaBinsHist = new TH1F ("SuccessfulThetaBinsHist",
827 "SuccessfulThetaBinsHist",
828 HistoNbins, BinsVector);
829
830 return kTRUE;
831
832}
833
834
835
836
837Bool_t MFindSupercutsONOFFThetaLoop::SetNexSigmaLiMaNormFactorNEvtsTrainTestHist()
838{
839 if (fCosThetaRangeVector.GetSize() <= 1)
840 {
841 *fLog << err
842 << "MFindSupercutsONOFFThetaLoop::SetNexSigmaLiMaNormFactorTrainTestHist: "
843 << endl
844 << "Size of fCosThetaRangeVector is <=1 ... "
845 << "fCosThetaRangeVector must exist and have, at least, "
846 << "2 components defining a single bin"
847 << endl;
848 return kFALSE;
849 }
850
851
852 Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
853
854 Double_t* BinsVector = fCosThetaRangeVector.GetArray();
855
856 char* x_axis_title = {"Cos(Theta)"};
857
858 fNormFactorTrainHist = new TH1F ("NormFactorTrainHist",
859 "NormFactorTrainHist",
860 HistoNbins, BinsVector);
861
862
863 fNormFactorTrainHist -> SetXTitle(x_axis_title);
864 fNormFactorTrainHist -> SetYTitle("Normalization Factor (Non/Noff)");
865
866
867 fNormFactorTestHist = new TH1F ("NormFactorTestHist",
868 "NormFactorTestHist",
869 HistoNbins, BinsVector);
870
871 fNormFactorTestHist -> SetXTitle(x_axis_title);
872 fNormFactorTestHist -> SetYTitle("Normalization Factor (Non/Noff)");
873
874 fSigmaLiMaTrainHist = new TH1F ("SigmaLiMaTrainHist",
875 "SigmaLiMaTrainHist",
876 HistoNbins, BinsVector);
877
878 fSigmaLiMaTrainHist -> SetXTitle(x_axis_title);
879 fSigmaLiMaTrainHist -> SetYTitle("Significance");
880
881
882 fSigmaLiMaTestHist = new TH1F ("SigmaLiMaTestHist",
883 "SigmaLiMaTestHist",
884 HistoNbins, BinsVector);
885
886 fSigmaLiMaTestHist -> SetXTitle(x_axis_title);
887 fSigmaLiMaTestHist -> SetYTitle("Significance");
888
889
890
891 fNexTrainHist = new TH1F ("NexTrainHist",
892 "NexTrainHist",
893 HistoNbins, BinsVector);
894
895 fNexTrainHist -> SetXTitle(x_axis_title);
896 fNexTrainHist -> SetYTitle("Excess Events");
897
898
899 fNexTestHist = new TH1F ("NexTestHist",
900 "NexTestHist",
901 HistoNbins, BinsVector);
902
903 fNexTestHist -> SetXTitle(x_axis_title);
904 fNexTestHist -> SetYTitle("Excess Events");
905
906
907
908 fNEvtsInTrainMatrixONHist = new TH1F("NEvtsInTrainMatrixONHist",
909 "NEvtsInTrainMatrixONHist",
910 HistoNbins, BinsVector);
911
912
913 fNEvtsInTrainMatrixONHist -> SetXTitle(x_axis_title);
914 fNEvtsInTrainMatrixONHist -> SetYTitle("Number of ON Events");
915
916
917 fNEvtsInTestMatrixONHist = new TH1F("NEvtsInTestMatrixONHist",
918 "NEvtsInTestMatrixONHist",
919 HistoNbins, BinsVector);
920
921 fNEvtsInTestMatrixONHist -> SetXTitle(x_axis_title);
922 fNEvtsInTestMatrixONHist -> SetYTitle("Number of ON Events");
923
924
925 fNEvtsInTrainMatrixOFFHist = new TH1F("NEvtsInTrainMatrixOFFHist",
926 "NEvtsInTrainMatrixOFFHist",
927 HistoNbins, BinsVector);
928
929 fNEvtsInTrainMatrixOFFHist -> SetXTitle(x_axis_title);
930 fNEvtsInTrainMatrixOFFHist -> SetYTitle("Number of OFF Events");
931
932
933 fNEvtsInTestMatrixOFFHist = new TH1F("NEvtsInTestMatrixOFFHist",
934 "NEvtsInTestMatrixOFFHist",
935 HistoNbins, BinsVector);
936
937 fNEvtsInTestMatrixOFFHist -> SetXTitle(x_axis_title);
938 fNEvtsInTestMatrixOFFHist -> SetYTitle("Number of OFF Events");
939
940
941
942
943
944 return kTRUE;
945
946}
947
948void MFindSupercutsONOFFThetaLoop::WriteSuccessfulThetaBinsHistToFile()
949{
950 *fLog << "MFindSupercutsONOFFThetaLoop : Writing histogram "
951 << " 'SuccessfulThetaBinsHist'"
952 << " (if contents != 0) into root file "
953 << fAlphaDistributionsRootFilename << endl;
954
955 TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE",
956 "Alpha Distributions for several Theta bins");
957
958
959 if (fSuccessfulThetaBinsHist)
960 {
961 if (fSuccessfulThetaBinsHist->GetEntries() > 0.5)
962 fSuccessfulThetaBinsHist -> Write();
963 }
964
965
966 rootfile.Close();
967
968}
969
970void MFindSupercutsONOFFThetaLoop::WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile()
971{
972
973 *fLog << "MFindSupercutsONOFFThetaLoop : Writing histograms 'NexSigmaLiMaNormFactorNEvtsTrainTest'"
974 << " (if they have contents != 0) into root file "
975 << fAlphaDistributionsRootFilename << endl;
976
977 TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE",
978 "Alpha Distributions for several Theta bins");
979
980
981 if (fNormFactorTrainHist)
982 {
983 if (fNormFactorTrainHist->GetEntries() > 0.5)
984 fNormFactorTrainHist -> Write();
985 }
986
987 if (fNormFactorTestHist)
988 {
989 if (fNormFactorTestHist->GetEntries() > 0.5)
990 fNormFactorTestHist -> Write();
991 }
992
993 if (fSigmaLiMaTrainHist)
994 {
995 if (fSigmaLiMaTrainHist->GetEntries() > 0.5)
996 fSigmaLiMaTrainHist -> Write();
997 }
998
999 if (fSigmaLiMaTestHist)
1000 {
1001 if (fSigmaLiMaTestHist->GetEntries() > 0.5)
1002 fSigmaLiMaTestHist -> Write();
1003 }
1004
1005 if (fNexTrainHist)
1006 {
1007 if (fNexTrainHist->GetEntries() > 0.5)
1008 fNexTrainHist -> Write();
1009 }
1010
1011 if (fNexTestHist)
1012 {
1013 if (fNexTestHist->GetEntries() > 0.5)
1014 fNexTestHist -> Write();
1015 }
1016
1017
1018 if (fNEvtsInTrainMatrixONHist)
1019 {
1020 if (fNEvtsInTrainMatrixONHist->GetEntries() > 0.5)
1021 fNEvtsInTrainMatrixONHist -> Write();
1022 }
1023
1024 if (fNEvtsInTestMatrixONHist)
1025 {
1026 if (fNEvtsInTestMatrixONHist->GetEntries() > 0.5)
1027 fNEvtsInTestMatrixONHist -> Write();
1028 }
1029
1030 if (fNEvtsInTrainMatrixOFFHist)
1031 {
1032 if (fNEvtsInTrainMatrixOFFHist->GetEntries() > 0.5)
1033 fNEvtsInTrainMatrixOFFHist -> Write();
1034 }
1035
1036
1037 if (fNEvtsInTestMatrixOFFHist)
1038 {
1039 if (fNEvtsInTestMatrixOFFHist->GetEntries() > 0.5)
1040 fNEvtsInTestMatrixOFFHist -> Write();
1041 }
1042
1043
1044 rootfile.Close();
1045
1046}
1047
1048void MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents(
1049 Double_t fontrain,
1050 Double_t fontest,
1051 Double_t fofftrain,
1052 Double_t fofftest)
1053{
1054
1055 // Check that fractions set by user are within the expected limits
1056
1057 if (fontrain < 0.00 || fontrain > 1.00 ||
1058 fontest < 0.00 || fontest > 1.00 ||
1059 fofftrain < 0.00 || fofftrain > 1.00 ||
1060 fofftest < 0.00 || fontest > 1.00)
1061 {
1062 *fLog << err
1063 << "MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents: "
1064 << endl
1065 << "Fractions of ON/OFF data for TRAIN/TEST sample which were specified by user " << endl
1066 << "are outside the natural limits 0.0-1.0. " << endl
1067 << "Set reasonable fraction numbers for TRAIN/TEST please. The program will be aborted."
1068 << endl;
1069
1070 exit(1);
1071
1072
1073 }
1074
1075
1076 if((fontrain+fontest) > 1.00 ||
1077 (fofftrain+fofftest) >1.00)
1078
1079 {
1080 *fLog << err
1081 << "MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents: "
1082 << endl
1083 << "The sum of the fractions of TRAIN/TEST for ON or OFF data which were specified by user " << endl
1084 << "are larger than 1.0. " << endl
1085 << "Set reasonable fraction numbers for TRAIN/TEST please. The program will be aborted."
1086 << endl;
1087
1088 exit(1);
1089
1090
1091 }
1092
1093
1094 fWhichFractionTrain = fontrain;
1095 fWhichFractionTest = fontest;
1096 fWhichFractionTrainOFF = fofftrain;
1097 fWhichFractionTestOFF = fofftest;
1098
1099
1100}
1101
1102void MFindSupercutsONOFFThetaLoop::SetReadMatricesFromFile (Bool_t b)
1103{
1104
1105 *fLog << "---------------------------------------------------" << endl
1106 << "MFindSupercutsONOFFThetaLoop: " << endl;
1107
1108 if (b)
1109 { *fLog << "Matrices are read from root files." << endl;}
1110 else
1111 { *fLog << "Matrices are produced and stored into root files."
1112 << endl;
1113 }
1114
1115
1116 fReadMatricesFromFile = b;
1117}
1118
1119Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaDistributionsRootFilename(
1120 const TString &name)
1121{
1122
1123 if (fPathForFiles.IsNull())
1124 {
1125 *fLog << err
1126 << "MFindSupercutsONOFFThetaLoop::SetAlphaDistributionsRootFilename; "
1127 << endl
1128 << "fPathForFiles is not defined yet, hence "
1129 << "name for rootfile containing alpha distributions "
1130 << "can not be defined properly" << endl;
1131
1132 return kFALSE;
1133 }
1134
1135
1136 fAlphaDistributionsRootFilename = (fPathForFiles);
1137 fAlphaDistributionsRootFilename += (name);
1138
1139
1140 *fLog << "MFindSupercutsONOFFThetaLoop: fAlphaDistributionsRootFilename = "
1141 << fAlphaDistributionsRootFilename << endl;
1142
1143
1144 return kTRUE;
1145
1146
1147}
1148
1149// Function to set Names of root files containing matrices and
1150// optimizedparameters and also fThetaRangeStringVector vector
1151
1152Bool_t MFindSupercutsONOFFThetaLoop::SetALLNames()
1153{
1154 // Names for files only possible if CosThetaRangeVector
1155 // is defined properly
1156 if (fCosThetaRangeVector.GetSize() < 2)
1157 {
1158 *fLog << err
1159 << "MFindSupercutsONOFFThetaLoop::SetALLNames; "
1160 << endl
1161 << "Size of fCosThetaRangeVector is smaller than 2... "
1162 << "fCosThetaRangeVector is not defined properly, "
1163 << "and thus, names for rootfiles containing data matrices "
1164 << "and optimized parameters for the different theta "
1165 << "ranges can not defined properly." << endl;
1166
1167 return kFALSE;
1168 }
1169
1170 if (fPathForFiles.IsNull())
1171 {
1172 *fLog << err
1173 << "MFindSupercutsONOFFThetaLoop::Setnames; "
1174 << endl
1175 << "fPathForFiles is not defined yet, hence "
1176 << "names for rootfiles containing data matrices "
1177 << "and optimized parameters for the different theta "
1178 << "ranges can not defined properly." << endl;
1179
1180 return kFALSE;
1181 }
1182
1183
1184 // Name of rootfiles for Supercuts parameters follow the
1185 // follwing structure
1186
1187 // Path + "OptSCParametersONOFF" + "ThetaRange"
1188 // + int {(fThetaMin_fThetaMax)} + ".root"
1189
1190
1191 // Name of rootfiles for Matrices follow the follwing structure
1192
1193 // Path + "ON(OFF)Data" + "Train(Test)" + "Matrix" + Fraction
1194 // + "ThetaRange" + int {(fThetaMin_fThetaMax)} + ".root"
1195
1196
1197 Int_t tmpdim = fCosThetaRangeVector.GetSize() - 1;
1198 const Int_t VectorDimension = tmpdim;
1199
1200 fThetaRangeStringVector = new TString[VectorDimension];
1201
1202 fOptSCParamFilenameVector = new TString[VectorDimension];
1203
1204 fTrainMatrixONFilenameVector = new TString[VectorDimension];
1205 fTestMatrixONFilenameVector = new TString[VectorDimension];
1206 fTrainMatrixOFFFilenameVector = new TString[VectorDimension];
1207 fTestMatrixOFFFilenameVector = new TString[VectorDimension];
1208
1209
1210 // Int_t ThetaMinTimes1000 = 0;
1211 // Int_t ThetaMaxTimes1000 = 0;
1212
1213 // For MAGIC new (July 2004) Theta is given in deg
1214
1215 Int_t ThetaMinInt = 0;
1216 Int_t ThetaMaxInt = 0;
1217
1218
1219 Int_t FractionONTrainTimes100 = 0;
1220 Int_t FractionONTestTimes100 = 0;
1221 Int_t FractionOFFTrainTimes100 = 0;
1222 Int_t FractionOFFTestTimes100 = 0;
1223
1224 char ThetaRangeString[100];
1225
1226 for (Int_t i = 0; i < VectorDimension; i++)
1227 {
1228 if (!SetThetaRange(i+1))
1229 {
1230 *fLog << "MFindSupercutsONOFFThetaLoop::SetALLNames; "
1231 << endl
1232 << "Values for fThetaMin and fThetaMax could NOT "
1233 << "be computed for theta bin "
1234 << (i+1) << "; SetNames can not go on successfully."
1235 << endl;
1236
1237 return kFALSE;
1238 }
1239
1240
1241 //ThetaMinTimes1000 = int(fThetaMin*1000);
1242 //ThetaMaxTimes1000 = int(fThetaMax*1000);
1243 //sprintf(ThetaRangeString, "ThetaRange%d_%dmRad",
1244 // ThetaMinTimes1000, ThetaMaxTimes1000);
1245
1246 // For New MAGIC (July 2004)
1247 // Theta is given in deg, no need to multiply by 1000
1248
1249 // 0.5 is added so that the rounding to integer is correct
1250 ThetaMinInt = int(fThetaMin + 0.5);
1251 ThetaMaxInt = int(fThetaMax + 0.5);
1252 sprintf(ThetaRangeString, "ThetaRange%d_%d_Degrees",
1253 ThetaMinInt, ThetaMaxInt);
1254
1255
1256
1257 fThetaRangeStringVector[i] = (ThetaRangeString);
1258
1259 // tmp
1260 *fLog << "ThetaRangeString = " << fThetaRangeStringVector[i] << endl;
1261 // endtmp
1262
1263 // 0.5 is added so that the rounding to integer is correct
1264
1265 FractionONTrainTimes100 = int(fWhichFractionTrain*100 + 0.5);
1266 FractionONTestTimes100 = int(fWhichFractionTest*100 + 0.5);
1267 FractionOFFTrainTimes100 = int(fWhichFractionTrainOFF*100 + 0.5);
1268 FractionOFFTestTimes100 = int(fWhichFractionTestOFF*100 + 0.5);
1269
1270
1271
1272
1273
1274 // File for SC parameters
1275 fOptSCParamFilenameVector[i] = (fPathForFiles);
1276 fOptSCParamFilenameVector[i] += ("OptSCParametersONOFF");
1277 fOptSCParamFilenameVector[i] += (fThetaRangeStringVector[i]);
1278 fOptSCParamFilenameVector[i] += (".root");
1279
1280
1281
1282 // File for Data Matrices
1283
1284 fTrainMatrixONFilenameVector[i] = (fPathForFiles);
1285 fTrainMatrixONFilenameVector[i] += ("ONDataTrainMatrixFraction");
1286 fTrainMatrixONFilenameVector[i] += (FractionONTrainTimes100);
1287 fTrainMatrixONFilenameVector[i] += (fThetaRangeStringVector[i]);
1288 fTrainMatrixONFilenameVector[i] += (".root");
1289
1290 fTestMatrixONFilenameVector[i] = (fPathForFiles);
1291 fTestMatrixONFilenameVector[i] += ("ONDataTestMatrixFraction");
1292 fTestMatrixONFilenameVector[i] += (FractionONTestTimes100);
1293 fTestMatrixONFilenameVector[i] += (fThetaRangeStringVector[i]);
1294 fTestMatrixONFilenameVector[i] += (".root");
1295
1296
1297 fTrainMatrixOFFFilenameVector[i] = (fPathForFiles);
1298 fTrainMatrixOFFFilenameVector[i] += ("OFFDataTrainMatrixFraction");
1299 fTrainMatrixOFFFilenameVector[i] += (FractionOFFTrainTimes100);
1300 fTrainMatrixOFFFilenameVector[i] += (fThetaRangeStringVector[i]);
1301 fTrainMatrixOFFFilenameVector[i] += (".root");
1302
1303 fTestMatrixOFFFilenameVector[i] = (fPathForFiles);
1304 fTestMatrixOFFFilenameVector[i] += ("OFFDataTestMatrixFraction");
1305 fTestMatrixOFFFilenameVector[i] += (FractionOFFTestTimes100);
1306 fTestMatrixOFFFilenameVector[i] += (fThetaRangeStringVector[i]);
1307 fTestMatrixOFFFilenameVector[i] += (".root");
1308
1309 }
1310 // Info concerning names is printed
1311
1312 *fLog << "--------------------------------------------------" << endl
1313 << "MFindSupercutsONOFFThetaLoop::Setnames; "
1314 << endl
1315 << "Names for root files containing Opt SC Param. and Data "
1316 << "matrices for the different theta bins are the "
1317 << "following ones: " << endl
1318 << "MeanCosTheta, OptSCParm, ONDataTrainMatrix, ..." << endl;
1319
1320 for (Int_t i = 0; i < VectorDimension; i++)
1321 {
1322 *fLog << fCosThetaBinCenterVector[i] << ", " << endl
1323 << fOptSCParamFilenameVector[i] << ", "<< endl
1324 << fTrainMatrixONFilenameVector[i] << ", " << endl
1325 << fTestMatrixONFilenameVector[i] << ", " << endl
1326 << fTrainMatrixOFFFilenameVector[i] << ", " << endl
1327 << fTestMatrixOFFFilenameVector[i] << ", " << endl
1328 << endl << endl;
1329 }
1330
1331 *fLog << "--------------------------------------------------" << endl;
1332
1333 return kTRUE;
1334
1335}
1336
1337
1338
1339Bool_t MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges()
1340{
1341
1342
1343 *fLog << "--------------------------------------------------------------" << endl
1344 << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges function starts "
1345 << endl;
1346
1347
1348
1349 if (fDataONRootFilename.IsNull())
1350 {
1351 *fLog << "MFindSupercutsONOFF::LoopOverThetaRanges; root file for ON data is not defined... "
1352 << "aborting..."
1353 << endl;
1354 return kFALSE;
1355 }
1356
1357
1358 if (fDataOFFRootFilename.IsNull())
1359 {
1360 *fLog << "MFindSupercutsONOFF:: LoopOverThetaRanges; root file for OFF data is not defined... aborting"
1361 << endl;
1362 return kFALSE;
1363 }
1364
1365
1366 if (fHadronnessName.IsNull())
1367 {
1368 *fLog << "MFindSupercutsONOFF::LoopOverThetaRanges; hadronness name for ON data is not defined... aborting"
1369 << endl;
1370 return kFALSE;
1371 }
1372
1373
1374 if (fHadronnessNameOFF.IsNull())
1375 {
1376 *fLog << "MFindSupercutsONOFF::LoopOverThetaRanges; hadronness name for OFF data is not defined... aborting"
1377 << endl;
1378 return kFALSE;
1379 }
1380
1381
1382 if (fCosThetaRangeVector.GetSize() < 2)
1383
1384 {
1385 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1386 << "Dimension of vector fCosThetaRangeVector "
1387 << "is smaller than 2"
1388 << "Theta ranges are not well defined... aborting ..." << endl;
1389
1390 return kFALSE;
1391 }
1392
1393
1394 // Check if pointers to root file names (opt SC parameters and matrices)
1395 // are empty
1396
1397
1398 if (!fTrainMatrixONFilenameVector ||
1399 !fTestMatrixONFilenameVector ||
1400 !fTrainMatrixOFFFilenameVector ||
1401 !fTestMatrixOFFFilenameVector ||
1402 !fOptSCParamFilenameVector)
1403 {
1404 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1405 << "Pointers to root files names for Opt SC parameters and data "
1406 << "matrices are not well defined. Aborting... " << endl;
1407
1408 return kFALSE;
1409 }
1410
1411
1412 if (fAlphaDistributionsRootFilename.IsNull())
1413 {
1414 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1415 << "Root filename where to store alpha distributions "
1416 << "(fAlphaDistributionsRootFilename) is not defined." << endl
1417 << "Loop over theta angles can not go on..." << endl;
1418
1419 return kTRUE;
1420 }
1421
1422
1423
1424
1425 if (!fPsFilename)
1426 {
1427 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1428 << "Object from Class TPostScript is not defined. "
1429 << "Plots will be stored in a a file named 'PsFilesTmp.ps' located "
1430 << "in current directory" << endl;
1431
1432 // fPsFilename = new TPostScript("PsFilesTmp.ps" ,111);
1433
1434
1435 }
1436
1437
1438
1439
1440
1441
1442 // All drinks are in house... the party starts !!
1443
1444 Int_t ThetaLoopIterations = fCosThetaRangeVector.GetSize() - 1;
1445
1446 *fLog << "looping over the following Cos Theta Angle Ranges : "<< endl;
1447 for (Int_t i = 0; i < ThetaLoopIterations; i++)
1448 {
1449 *fLog << i+1 << " -> "
1450 << fCosThetaRangeVector[i] << "-" << fCosThetaRangeVector[i+1]
1451 << endl;
1452 }
1453
1454 *fLog << "--------------------------------------------------------------------"
1455 << endl;
1456
1457
1458
1459 // Let's create histograms where SigmaLiMa, Nex and Normfactors will be stored
1460
1461 if (!SetNexSigmaLiMaNormFactorNEvtsTrainTestHist())
1462 {
1463 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1464 << "Histograms where to store Significances, Nex and Normalization "
1465 << "factors, NEvts for all theta bins could not be created... Aborting."
1466 << endl;
1467
1468 return kFALSE;
1469 }
1470
1471
1472 // Let's create histogram where successfully optimized theta bins will be stored
1473
1474 if (!SetSuccessfulThetaBinsHist())
1475 {
1476 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1477 << "Histogram where to store successfully optimized theta bins "
1478 << " could not be created... Aborting."
1479 << endl;
1480
1481 return kFALSE;
1482 }
1483
1484
1485
1486
1487
1488 Double_t NormFactor = 0.0;
1489 Double_t Nex = 0.0;
1490 Double_t SigmaLiMa = 0.0;
1491 Int_t NEvtsON = 0;
1492 Int_t NEvtsOFF = 0;
1493
1494 Bool_t MinimizationSuccessful = kTRUE;
1495 Double_t MinimizationValue = 1.0; // Value quantifying the minimization, and stored in
1496 // histogram fSuccessfulThetaBinsHist. For the time being it is 0.0 (minimization failed).
1497 // and 1.0 (minimization succeeded). Yet in future it may take other values quantifying
1498 // other aspects...
1499
1500
1501 for (Int_t i = 0; i < ThetaLoopIterations; i++)
1502 {
1503
1504 // Get the Theta range that will be used in this iteration
1505 // Remember that argument of the function is iteration number,
1506 // and it starts at 1...
1507
1508 if (!SetThetaRange(i+1))
1509 {
1510 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; "
1511 << endl
1512 << "Values for fThetaMin and fThetaMax could NOT "
1513 << "be computed for theta bin "
1514 << (i+1) << "; LoopOverThetaRanges is aborting."
1515 << endl;
1516
1517 return kFALSE;
1518 }
1519
1520
1521
1522 // Object of class MFindSupercutsONOFF will be created and used
1523
1524
1525 MFindSupercutsONOFF FindSupercuts("MFindSupercutsONOFF",
1526 "Optimizer for the supercuts");
1527
1528
1529 // Set Theta range that will be used
1530 FindSupercuts.SetThetaRange(fThetaMin, fThetaMax);
1531
1532
1533 // Alphamax where signal is expected is set
1534 FindSupercuts.SetAlphaSig(fAlphaSig);
1535
1536 // Bkg alpha region is set
1537 FindSupercuts.SetAlphaBkgMin(fAlphaBkgMin);
1538 FindSupercuts.SetAlphaBkgMax(fAlphaBkgMax);
1539
1540 // alpha bkg and signal region set in object FindSupercuts
1541 // are re-checked in order to be sure that make sense
1542
1543 FindSupercuts.CheckAlphaSigBkg();
1544
1545
1546 // Degree of polynomials used to fit ON and OFF alpha
1547 // distributions is set
1548
1549 FindSupercuts.SetDegreeON(fDegreeON);
1550 FindSupercuts.SetDegreeOFF(fDegreeOFF);
1551
1552
1553
1554
1555 // bining for alpha plots is set
1556
1557 FindSupercuts.SetAlphaPlotBinining(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
1558
1559
1560
1561
1562 // Bool variable NormFactorFromAlphaBkg is set in FindSupercuts
1563 // in order to decide the way in which mormalization factors
1564 // for TRAIN and TEST samples are computed
1565
1566
1567 FindSupercuts.SetNormFactorFromAlphaBkg(fNormFactorFromAlphaBkg);
1568
1569
1570 // Define size range of events used to fill the data matrices
1571
1572 FindSupercuts.SetSizeRange(fSizeCutLow, fSizeCutUp);
1573
1574 FindSupercuts.SetFilters(fLeakMax, fDistMax, fDistMin);
1575
1576
1577
1578 // Boolean variable that controls wether the optimization of the
1579 // parameters (MMinuitInterface::CallMinuit(..) in function FindParams(..))
1580 // takes place or not. kTRUE will skip such optimization.
1581 // This variable is useful to test the optmized parameters (previously found
1582 // and stored in root file) on the TRAIN sample.
1583
1584 FindSupercuts.SetSkipOptimization(fSkipOptimization);
1585
1586 // Boolean variable that allows the user to write the initial parameters
1587 // into the root file that will be used to store the optimum cuts.
1588 // If fUseInitialSCParams = kTRUE , parameters are written.
1589 // In this way, the initial SC parameters can be applied on the data (train/test)
1590
1591 FindSupercuts.SetUseInitialSCParams(fUseInitialSCParams);
1592
1593
1594 // Set boolean variable that controls wether the theta variable
1595 // is used or not in the computation of the dynamica cuts
1596
1597 FindSupercuts.SetVariableNotUseTheta(fNotUseTheta);
1598
1599 // Set boolean variable that controls wether cuts are dynamical
1600 // or static.
1601
1602 FindSupercuts.SetVariableUseStaticCuts(fUseStaticCuts);
1603
1604
1605 FindSupercuts.SetUseFittedQuantities(fUseFittedQuantities);
1606
1607
1608 // fGammaEfficiency is set in object FindSupercuts
1609
1610 FindSupercuts.SetGammaEfficiency (fGammaEfficiency);
1611
1612
1613 // Filename with SuperCuts parameters is set
1614
1615 FindSupercuts.SetFilenameParam(fOptSCParamFilenameVector[i]);
1616
1617 // Hadroness containers are set
1618
1619 FindSupercuts.SetHadronnessName(fHadronnessName);
1620 FindSupercuts.SetHadronnessNameOFF(fHadronnessNameOFF);
1621
1622 // Rootfilename where to store alpha distribution histograms is set
1623
1624 FindSupercuts.SetAlphaDistributionsRootFilename(fAlphaDistributionsRootFilename);
1625
1626
1627 // Set the names of all trees
1628 // used to store the supercuts applied to ON/OFF Train/Test samples
1629
1630 // In future, names might be controlled by user, but for the time
1631 // being I'll make things simple, and names are set to some defaults
1632 // taking into account the theta range used
1633
1634 // Containers will be stored in rootfile specified by variable
1635 // fAlphaDistributionsRootFilename
1636
1637 FindSupercuts.SetSupercutsAppliedTreeNames();
1638
1639
1640
1641
1642
1643 // Object of the class TPostScritp is set...
1644 // BUT THAT'S NOT WORKING PROPERLY !!!!!
1645
1646 // FindSupercuts.SetPostScriptFile (fPsFilename);
1647
1648 // ********************************************************
1649 // Due to the failure of the use of object TPostScript
1650 // to make a Ps document with all plots, I decided to use the
1651 // standard way (SaveAs(filename.ps)) to store plots related
1652 // to alpha distributions
1653 // for ON and OFF and BEFORE and AFTER cuts (VERY IMPORTANT
1654 // TO COMPUTE EFFICIENCIES IN CUTS)
1655
1656 // Psfilename is set inside function LoopOverThetaRanges()
1657 // and given to the object MFindSupercutsONOFF created
1658 // within this loop.
1659
1660 // This will have to be removed as soon as the TPostScript
1661 // solutions works...
1662 // ********************************************************
1663
1664 TString PsFilenameString = (fPathForFiles);
1665 PsFilenameString += ("SC");
1666 PsFilenameString += (fThetaRangeStringVector[i]);
1667 // Name to be completed inside object MFindSupercutsONOFF
1668
1669 FindSupercuts.SetPsFilenameString(PsFilenameString);
1670
1671
1672 if (fReadMatricesFromFile)
1673 {
1674 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1675 << "Reading Data matrices from root files... " << endl;
1676
1677 /*
1678 FindSupercuts.ReadMatrix(fTrainMatrixONFilenameVector[i],
1679 fTestMatrixONFilenameVector[i]);
1680 FindSupercuts.ReadMatrixOFF(fTrainMatrixOFFFilenameVector[i],
1681 fTestMatrixOFFFilenameVector[i]);
1682 */
1683
1684 FindSupercuts.ReadMatrixTrain(fTrainMatrixONFilenameVector[i],
1685 fTrainMatrixOFFFilenameVector[i]);
1686
1687 if (fTestParameters)
1688 {
1689 FindSupercuts.ReadMatrixTest(fTestMatrixONFilenameVector[i],
1690 fTestMatrixOFFFilenameVector[i]);
1691 }
1692
1693 }
1694
1695 else
1696 {
1697 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1698 << "Reading Data from files " << fDataONRootFilename
1699 << " and " << fDataOFFRootFilename << "(for ON and OFF data "
1700 << "respectively "
1701 << " and producing data matrices that will be stored "
1702 << "in root files. " << endl;
1703
1704 // ON data
1705
1706 FindSupercuts.DefineTrainTestMatrixThetaRange(
1707 fDataONRootFilename,
1708 fWhichFractionTrain,
1709 fWhichFractionTest,
1710 fThetaMin, fThetaMax,
1711 fTrainMatrixONFilenameVector[i],
1712 fTestMatrixONFilenameVector[i]);
1713
1714 // OFF data
1715
1716 FindSupercuts.DefineTrainTestMatrixOFFThetaRange(
1717 fDataOFFRootFilename,
1718 fWhichFractionTrainOFF,
1719 fWhichFractionTestOFF,
1720 fThetaMin, fThetaMax,
1721 fTrainMatrixOFFFilenameVector[i],
1722 fTestMatrixOFFFilenameVector[i]);
1723
1724 }
1725
1726
1727 MinimizationSuccessful = kTRUE;
1728
1729 if (fOptimizeParameters)
1730 {
1731 MinimizationSuccessful =
1732 FindSupercuts.FindParams("",fInitSCPar, fInitSCParSteps);
1733
1734 if (MinimizationSuccessful)
1735 {
1736 // Minimization was successful
1737 // NormFactor, Nex, SigmaLiMa and NEvts Train histograms are updated
1738
1739 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; "
1740 << endl
1741 << "Minimization was successful" << endl
1742 << "Updating values for NormFactor, Nex, SigmaLiMa and NEvts "
1743 << "Train histograms for " << fThetaRangeStringVector[i] << endl;
1744
1745
1746 NormFactor = FindSupercuts.GetNormFactorTrain();
1747 SigmaLiMa = FindSupercuts.GetSigmaLiMaTrain();
1748 Nex = FindSupercuts.GetNexTrain();
1749
1750 NEvtsON = FindSupercuts.GetMatrixTrain()->GetM().GetNrows();
1751 NEvtsOFF = FindSupercuts.GetMatrixTrainOFF()->GetM().GetNrows();
1752
1753 // MinimizationValue is updated...
1754
1755 MinimizationValue = 1.0;
1756
1757 }
1758 else
1759 {
1760 // Minimization was NOT successful
1761 // NormFactor, Nex, SigmaLiMa and NEvts Train histograms are updated
1762
1763 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; "
1764 << endl
1765 << "Minimization was NOT successful" << endl
1766 << "Updating values for NormFactor, and NEvts "
1767 << "Train histograms for " << fThetaRangeStringVector[i] << endl
1768 << "Nex and SigmaLiMa are set to ZERO" << endl;
1769
1770
1771 NormFactor = FindSupercuts.GetNormFactorTrain();
1772
1773 NEvtsON = FindSupercuts.GetMatrixTrain()->GetM().GetNrows();
1774 NEvtsOFF = FindSupercuts.GetMatrixTrainOFF()->GetM().GetNrows();
1775
1776 SigmaLiMa = 0.0;
1777 Nex = 0.0;
1778
1779 // MinimizationValue is updated...
1780
1781 MinimizationValue = 0.0;
1782
1783 }
1784
1785
1786
1787 cout << "---- Train Histograms will be updated with the following numbers: ------ "
1788 << endl
1789 << "Mean CosTheta for this Thetabin = "
1790 << fCosThetaBinCenterVector[i] << endl
1791 << "Minimization status (1.0 is successful, 0.0 is UNsuccessful) = "
1792 << MinimizationValue << endl
1793 << "Number of ON Events in Train sample used = "
1794 << NEvtsON << endl
1795 << "Number of OFF Events in Train sample used = "
1796 << NEvtsOFF << endl
1797 << "Normalization factor (Non/Noff) for Train sample = "
1798 << NormFactor << endl
1799 << "Excess events and SigmaLiMa: "
1800 << Nex << "; " << SigmaLiMa << endl << endl;
1801
1802
1803
1804 fNormFactorTrainHist -> Fill(fCosThetaBinCenterVector[i],
1805 NormFactor);
1806
1807 fNexTrainHist -> Fill(fCosThetaBinCenterVector[i],
1808 Nex);
1809
1810 fSigmaLiMaTrainHist -> Fill(fCosThetaBinCenterVector[i],
1811 SigmaLiMa);
1812
1813
1814 fNEvtsInTrainMatrixONHist -> Fill(fCosThetaBinCenterVector[i],
1815 NEvtsON);
1816
1817 fNEvtsInTrainMatrixOFFHist -> Fill(fCosThetaBinCenterVector[i],
1818 NEvtsOFF);
1819
1820
1821 fSuccessfulThetaBinsHist -> Fill(fCosThetaBinCenterVector[i],
1822 MinimizationValue);
1823
1824 cout << "Histograms (TRAIN) updated successfully " << endl;
1825
1826
1827
1828
1829 }
1830
1831
1832 if (fTestParameters && MinimizationSuccessful)
1833 {
1834
1835 // Cuts are applied to test sample
1836
1837 FindSupercuts.TestParamsOnTestSample();
1838
1839 // NormFactor, Nex and SigmaLiMa Test histograms are updated
1840
1841 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1842 << "Updating values for NormFactor, Nex, SigmaLiMa and NEvts "
1843 << "Test histograms for " << fThetaRangeStringVector[i] << endl;
1844
1845
1846 NormFactor = FindSupercuts.GetNormFactorTest();
1847 SigmaLiMa = FindSupercuts.GetSigmaLiMaTest();
1848 Nex = FindSupercuts.GetNexTest();
1849
1850 NEvtsON = FindSupercuts.GetMatrixTest()->GetM().GetNrows();
1851 NEvtsOFF = FindSupercuts.GetMatrixTestOFF()->GetM().GetNrows();
1852
1853
1854 cout << "---- Test Histograms will be updated with the following numbers: ------ "
1855 << endl
1856 << "Mean CosTheta for this Thetabin = "
1857 << fCosThetaBinCenterVector[i] << endl
1858 << "Number of ON Events in Test sample used = "
1859 << NEvtsON << endl
1860 << "Number of OFF Events in Test sample used = "
1861 << NEvtsOFF << endl
1862 << "Normalization factor (Non/Noff) for Test sample = "
1863 << NormFactor << endl
1864 << "Excess events and SigmaLiMa: "
1865 << Nex << "; " << SigmaLiMa << endl << endl;
1866
1867
1868
1869
1870 fNormFactorTestHist -> Fill(fCosThetaBinCenterVector[i],
1871 NormFactor);
1872
1873 fNexTestHist -> Fill(fCosThetaBinCenterVector[i],
1874 Nex);
1875
1876 fSigmaLiMaTestHist -> Fill(fCosThetaBinCenterVector[i],
1877 SigmaLiMa);
1878
1879
1880 fNEvtsInTestMatrixONHist -> Fill(fCosThetaBinCenterVector[i],
1881 NEvtsON);
1882
1883 fNEvtsInTestMatrixOFFHist -> Fill(fCosThetaBinCenterVector[i],
1884 NEvtsOFF);
1885
1886
1887 cout << "Histograms (TEST) updated successfully " << endl;
1888
1889
1890
1891
1892 }
1893
1894 }
1895
1896
1897 // PsFile is closed
1898
1899
1900 //*fLog << "Closing PsFile ";
1901 //*fLog << fPsFilename -> GetName() << endl;
1902
1903 //fPsFilename -> Close();
1904
1905
1906
1907 if (fOptimizeParameters || fTestParameters)
1908 {
1909 // Histograms cotanining SigmaLiMa, Nex and NormFactors and NEvts are written
1910 // in the same root file where alpha distributions where stored
1911 // i.e. fAlphaDistributionsRootFilename
1912
1913 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1914 << "Executing WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile()"
1915 << endl;
1916
1917 WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile();
1918
1919
1920 if (fOptimizeParameters)
1921 { // Histogram containing minimization status for all theta bins defined
1922 // in vector fCosThetaVector is stored in file fAlphaDistributionsRootFilename
1923
1924 *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
1925 << "Executing WriteSuccessfulThetaBinsHistToFile()"
1926 << endl;
1927
1928 WriteSuccessfulThetaBinsHistToFile();
1929
1930 }
1931
1932 }
1933
1934
1935 return kTRUE;
1936}
1937
1938
1939
1940
1941
1942 // Function that loops over the alpha distributions (ON-OFF)
1943 // stored in root file defined by fAlphaDistributionsRootFilename
1944 // and computes the significance and Nex (using MHFindSignificanceONOFF::FindSigma)
1945 // for several cuts in alpha (0-fAlphaSig; in bins defined for alpha distributions
1946 // by user).
1947
1948 // It initializes the histograms (memeber data), fills them and store them
1949 // in root file defined by fAlphaDistributionsRootFilename. A single histogram
1950 // for each theta bin.
1951
1952
1953
1954Bool_t MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()
1955{
1956
1957
1958 *fLog << endl << endl
1959 << "****************************************************************************" << endl
1960 << "Computing number of excess events and significance vs " << endl
1961 << " the cut in alpha parameter (fAlphaSig) using function " << endl
1962 << " MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig() " << endl
1963 << "****************************************************************************" << endl
1964 << endl;
1965
1966
1967 // Check that all elements needed by the function are available
1968
1969
1970 if (fAlphaDistributionsRootFilename.IsNull())
1971 {
1972 *fLog << err
1973 << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig; "
1974 << endl
1975 << "fAlphaDistributionsRootFilename is not defined yet, hence "
1976 << "histograms containing alpha distributions for the different theta "
1977 << "ranges can not be retrieved. " << endl;
1978
1979 return kFALSE;
1980 }
1981
1982
1983
1984 if (fCosThetaRangeVector.GetSize() < 2)
1985 {
1986 *fLog << err
1987 << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig; "
1988 << endl
1989 << "Size of fCosThetaRangeVector is smaller than 2... "
1990 << "fCosThetaRangeVector is not defined properly, "
1991 << "and thus, number of theta loops and names "
1992 << " of the histograms containing alpha distributions "
1993 << "for the different theta "
1994 << "ranges can not defined properly." << endl;
1995
1996 return kFALSE;
1997 }
1998
1999
2000
2001 // It seems that i have all I need for the time being. Let's start the party...
2002
2003
2004 Int_t tmpint = 0; // variable that will be used to store temporally integers
2005 Double_t tmpdouble = 0.0; // variable that will be used to store temporally doubles
2006
2007 tmpint = fCosThetaRangeVector.GetSize() - 1;
2008 const Int_t NThetaBins = tmpint;
2009
2010
2011 // Definition of the bins for Nex and Significance vs Alpha histograms
2012 // The bins in Alpha will start with ZERO and will end
2013 // with 5 bins further after fAlphaSig.
2014
2015 Double_t AlphaBinWidth = double((fAlphaBinUp-fAlphaBinLow)/fNAlphaBins);
2016 Double_t FirstBinLow = 0.0;
2017 tmpint = int((fAlphaSig/AlphaBinWidth) + 5.0);
2018
2019 while ((FirstBinLow + AlphaBinWidth*tmpint) > fAlphaBkgMin)
2020 {
2021 tmpint--;
2022 }
2023
2024 *fLog << "Information about the binning of the histograms Nex and Significance vs alphasig "
2025 << endl
2026 << "FirstBinLow, AlphaBinWidth, NBins : "
2027 << FirstBinLow << ", " << AlphaBinWidth << ", " << tmpint-1 << endl;
2028
2029
2030
2031 const Int_t NexSignVSAlphaXBinsVectorDimension = tmpint;
2032
2033
2034 Double_t* XBinsVector = new Double_t [NexSignVSAlphaXBinsVectorDimension];
2035
2036 for (Int_t i = 0; i < NexSignVSAlphaXBinsVectorDimension; i++)
2037 {
2038 XBinsVector[i] = FirstBinLow + i*AlphaBinWidth;
2039 // cout << XBinsVector[i] << endl;
2040 }
2041
2042
2043
2044 // Axis labels for Nex and Significance vs Alpha histograms
2045
2046 TString x_axis_label = ("Cut in alpha [\\circ]");
2047 TString y_axis_label_nex = ("Excess Events");
2048 TString y_axis_label_significance = ("Significance (LiMa formula 17)");
2049
2050
2051 // Arguments for function MHFindSignificanceONOFF::FindSigmaONOFF()
2052 // are defined in here. The same arguments (obviously) will be used
2053 // in all theta bins and TRAIN and TEST samples.
2054
2055 // variable alphasig will be defined and varied withing the theta loop
2056
2057 Double_t alphasig = 0.0;
2058
2059 // Minimum value of alpha for bkg region in ON data
2060 const Double_t alphabkgmin = fAlphaBkgMin;
2061
2062 // Maximum value of alpha for bkg region in ON data
2063 const Double_t alphabkgmax = fAlphaBkgMax;
2064
2065
2066 const Bool_t drawpoly = kFALSE;
2067 const Bool_t fitgauss = kFALSE;
2068 const Bool_t print = kFALSE;
2069
2070 const Bool_t saveplots = kFALSE; // No plots are saved
2071
2072 const Bool_t RebinHistograms = kTRUE;
2073 const Bool_t ReduceDegree = kFALSE;
2074
2075
2076
2077 Double_t Nex = 0.0;
2078 Double_t Significance = 0.0;
2079
2080
2081 // Names of histograms containing Nex and Significance vs alphasig
2082
2083 TString NexVSAlphaName = ("NexVSAlpha"); // Names to be completed in theta loop
2084 TString SignificanceVSAlphaName = ("SignificanceVSAlpha");// Names to be completed in theta loop
2085
2086
2087
2088 // Names of histograms to be retrieved from root file
2089 // fAlphaDistributionsRootFilename
2090
2091 TString NormFactorTrainHistoName = ("NormFactorTrainHist");
2092 TString NormFactorTestHistoName = ("NormFactorTestHist");
2093
2094 TString SuccessfulThetaBinsHistoName = ("SuccessfulThetaBinsHist");
2095
2096
2097
2098 // Histogram containing information about the success of the optimization
2099 // for the several theta bins. The content of those bins where optimization was not
2100 // successful, are set to ZERO.
2101
2102
2103 TH1F* SuccessfulThetaBinsHisto;
2104
2105
2106
2107 // Histograms to store the normalization factors
2108
2109 TH1F* NormFactorTrainHisto;
2110 TH1F* NormFactorTestHisto;
2111
2112
2113
2114 // Histograms that will be used to store the alpha distributions
2115 // retrieved from fAlphaDistributionsRootFilename, and that will be use as arguments in
2116 // function MHFindSignificanceONOFF::FindSigma
2117
2118 // The same histograms (pointer to histograms actually) will be used for TRAIN and TEST
2119 // and for the different ThetaNBinsUsed.
2120
2121 TH1F* AlphaONAfterCuts;
2122 TH1F* AlphaOFFAfterCuts;
2123
2124
2125
2126
2127 // Vector containing the histogram index of the successfully
2128 // optimized theta bins.
2129 // It is filled with the contents of the histogram
2130 // fSuccessfulThetaBinsHist retrieved from the root
2131 // file fAlphaDistributionsRootFilename
2132 // Remember that histogram index START AT 1 !!!
2133 TArrayI SuccessfulThetaBinsVector;
2134
2135
2136
2137
2138
2139
2140
2141 // Open root file for retrieving information
2142
2143 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl
2144 << "Opening root file " << fAlphaDistributionsRootFilename << endl;
2145
2146 TFile rootfile (fAlphaDistributionsRootFilename, "READ");
2147
2148 // Retrieving SuccessfulThetaBinsHisto from root file.
2149
2150 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl
2151 << "Retrieving SuccessfulThetaBinsHisto from root file and filling vector "
2152 << " SuccessfulThetaBinsVector" << endl;
2153
2154
2155 cout << "Getting histogram with name " << SuccessfulThetaBinsHistoName << endl;
2156
2157 SuccessfulThetaBinsHisto = (TH1F*) rootfile.Get(SuccessfulThetaBinsHistoName);
2158
2159
2160
2161 // Check that bins in this histo correspond with the
2162 // CosTheta intervals defined by vector fCosThetaRangeVector
2163
2164
2165 if (SuccessfulThetaBinsHisto -> GetNbinsX() != NThetaBins)
2166 {
2167 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()"
2168 << endl
2169 << "Number of theta bins defined by vector fCosThetaRangeVector ("
2170 << NThetaBins << ") does not match with the bins in histogram "
2171 << SuccessfulThetaBinsHistoName << "("
2172 << SuccessfulThetaBinsHisto -> GetNbinsX() << ")" << endl
2173 << "Function execution will be ABORTED." << endl;
2174 return kFALSE;
2175 }
2176
2177
2178 // Filling vector SuccessfulThetaBinsVector
2179
2180 tmpint = 0;
2181 for (Int_t i = 0; i < NThetaBins; i++)
2182 {
2183 if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5)
2184 {
2185 tmpint++;
2186 }
2187 }
2188
2189 SuccessfulThetaBinsVector.Set(tmpint);
2190
2191 tmpint = 0;
2192 for (Int_t i = 0; i < NThetaBins; i++)
2193 {
2194 if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5)
2195 {
2196 if(tmpint < SuccessfulThetaBinsVector.GetSize())
2197 {
2198 SuccessfulThetaBinsVector[tmpint] = i+1;
2199 tmpint++;
2200 }
2201 else
2202 {
2203 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()"
2204 << endl
2205 << "Problem when filling vector 'SuccessfulThetaBinsVector'"
2206 << endl
2207 << "Function execution will be ABORTED." << endl;
2208 return kFALSE;
2209 }
2210
2211
2212 }
2213 }
2214
2215 // variable defining the theta bins where optimization was successful, and
2216 // hence, alpha distributions stored in root file are meaningful (up to some extent)
2217
2218 tmpint = SuccessfulThetaBinsVector.GetSize();
2219 const Int_t SuccessfulNThetaBins = tmpint;
2220
2221 // TEMP
2222 // cout << "const Int_t SuccessfulNThetaBins = " << SuccessfulNThetaBins << endl;
2223 // ENDTEMP
2224
2225 // Vectors of pointers to histograms are defined to store the histograms
2226 // containing Nex and significance for TRAIN and TEST sample for all the
2227 // successful theta bins
2228
2229 // Histograms will be initialized in the loop over theta bins
2230
2231 TH1F* NexVSAlphaSigTrainVector[SuccessfulNThetaBins];
2232 TH1F* NexVSAlphaSigTestVector[SuccessfulNThetaBins];
2233
2234
2235 TH1F* SignificanceVSAlphaSigTrainVector[SuccessfulNThetaBins];
2236 TH1F* SignificanceVSAlphaSigTestVector[SuccessfulNThetaBins];
2237
2238
2239 // ***************************************************
2240 // HISTOGRAMS FOR TRAIN SAMPLE ARE COMPUTED
2241 // ***************************************************
2242
2243
2244
2245 if (fOptimizeParameters)
2246 { // Computation for TRAIN sample
2247
2248
2249 *fLog << endl << endl
2250 << "****************************************************************************" << endl
2251 << "Computing Nex and significance vs cut in alpha (alphasig) for TRAIN sample" << endl
2252 << "****************************************************************************" << endl
2253 << endl << endl;
2254
2255
2256 // Retrieving NormFactorTrainHisto from root file.
2257
2258 cout << "Getting histogram with name " << NormFactorTrainHistoName << endl;
2259
2260 NormFactorTrainHisto = (TH1F*) rootfile.Get(NormFactorTrainHistoName);
2261
2262
2263 // Check that bins in this histo correspond with the
2264 // CosTheta intervals defined by vector fCosThetaRangeVector
2265
2266 ////////////////// START CHECK //////////////////////////////////////
2267
2268 if (NormFactorTrainHisto -> GetNbinsX() != NThetaBins)
2269 {
2270 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()"
2271 << endl
2272 << "Number of theta bins defined by vector fCosThetaRangeVector ("
2273 << NThetaBins << ") does not match with the bins in histogram "
2274 << NormFactorTrainHistoName << "("
2275 << NormFactorTrainHisto -> GetNbinsX() << ")" << endl
2276 << "Function execution will be ABORTED." << endl;
2277 return kFALSE;
2278 }
2279
2280 // histo test
2281 /*
2282 cout << "NORMFACTORTRAINHIST Test: BinCenter and Value" << endl;
2283 for (Int_t k = 0; k < NThetaBins; k++)
2284 {
2285
2286 cout << NormFactorTrainHisto -> GetBinCenter(k+1)
2287 << "; " << NormFactorTrainHisto -> GetBinContent(k+1)<< endl;
2288 }
2289 */
2290
2291 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
2292 {
2293 tmpdouble = NormFactorTrainHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]);
2294 if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble ||
2295 fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble)
2296 {
2297 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()"
2298 << endl
2299 << "Bins defined by vector fCosThetaRangeVector "
2300 << "do not match with the bins of histogram "
2301 << NormFactorTrainHistoName << endl
2302 << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1]
2303 << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl
2304 << "NormFactorHist Bin: " << tmpdouble << endl
2305 << "Function execution will be ABORTED." << endl;
2306
2307 return kFALSE;
2308 }
2309
2310 }
2311
2312 ////////////////// END CHECK //////////////////////////////////////
2313
2314
2315
2316
2317 // Check that the theta range string (needed to compute names for alpha histograms)
2318 // is well defined
2319
2320 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
2321 {
2322 if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
2323 {
2324 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()"
2325 << endl
2326 << "Component " << SuccessfulThetaBinsVector[i]-1
2327 << " of fThetaRangeStringVector is EMPTY"
2328 << endl
2329 << "Function execution will be ABORTED." << endl;
2330
2331 return kFALSE;
2332
2333
2334 }
2335 }
2336
2337
2338 // *************************************************************
2339 // *************************************************************
2340
2341 // Normfactors and successfultheta bins are ok.
2342 // I can now compute the Nex and Sigma for train sample (ON-OFF)
2343
2344 // *************************************************************
2345 // *************************************************************
2346
2347
2348
2349
2350 // Loop over the several successful theta bins doing the following
2351
2352 // 0) Initializing histograms Nex and Significance (with names!!)
2353 // 1) Retrieving alpha histograms (ON-OFF), and normalization factor
2354 // 2) Computing Nex and Signficance vs alpha and filling histograms
2355
2356
2357 *fLog << endl
2358 << "Starting loop over the successfully optimized theta ranges " << endl
2359 << endl;
2360
2361 Double_t NormalizationFactor = 0.0;
2362
2363 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
2364 {
2365
2366 *fLog << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] << endl << endl;
2367
2368 // 0) Initializing histograms Nex and Significance (with names!!)
2369
2370 TString NexVSAlphaNameTmp = (NexVSAlphaName);
2371 TString SignificanceVSAlphaNameTmp = (SignificanceVSAlphaName);
2372
2373 NexVSAlphaNameTmp += "Train";
2374 NexVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
2375
2376 SignificanceVSAlphaNameTmp += "Train";
2377 SignificanceVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
2378
2379
2380 // TEMP
2381 //cout << NexVSAlphaNameTmp << endl;
2382 //cout << SignificanceVSAlphaNameTmp << endl;
2383 // ENDTEMP
2384
2385
2386
2387
2388 NexVSAlphaSigTrainVector[i] = new TH1F (NexVSAlphaNameTmp.Data(),
2389 NexVSAlphaNameTmp.Data(),
2390 NexSignVSAlphaXBinsVectorDimension - 1,
2391 XBinsVector);
2392
2393 // It is important to unlink the the created histograms from
2394 // the current directory (defined now by "rootfile"),
2395 // so that we can do whatever we want to
2396 // with them regardless of such directory is opened or closed
2397 // The reference of the histogram is removed from the current
2398 // directory by means of the member function "SetDirectory(0)"
2399
2400
2401 NexVSAlphaSigTrainVector[i] -> SetDirectory(0);
2402
2403
2404 SignificanceVSAlphaSigTrainVector[i] = new TH1F (SignificanceVSAlphaNameTmp.Data(),
2405 SignificanceVSAlphaNameTmp.Data(),
2406 NexSignVSAlphaXBinsVectorDimension - 1,
2407 XBinsVector);
2408
2409 SignificanceVSAlphaSigTrainVector[i] -> SetDirectory(0);
2410
2411
2412
2413 // axis are defined
2414
2415 NexVSAlphaSigTrainVector[i] -> SetXTitle(x_axis_label.Data());
2416 NexVSAlphaSigTrainVector[i] -> SetYTitle(y_axis_label_nex.Data());
2417
2418 SignificanceVSAlphaSigTrainVector[i] ->SetXTitle(x_axis_label.Data());
2419 SignificanceVSAlphaSigTrainVector[i] ->SetYTitle(y_axis_label_significance.Data());
2420
2421
2422 // 1) Retrieving alpha histograms (ON-OFF) and normalization factor
2423
2424 NormalizationFactor =
2425 NormFactorTrainHisto -> GetBinContent(SuccessfulThetaBinsVector[i]);
2426
2427
2428 TString TmpAlphaHistoName = ("TrainAlphaAfterCuts");
2429 TmpAlphaHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
2430
2431 *fLog << "Getting Histogram with name " << endl
2432 << TmpAlphaHistoName << endl;
2433
2434 AlphaONAfterCuts = (TH1F*) rootfile.Get(TmpAlphaHistoName);
2435
2436
2437 TString TmpAlphaOFFHistoName = ("TrainAlphaOFFAfterCuts");
2438 TmpAlphaOFFHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
2439
2440 *fLog << "Getting Histogram with name " << endl
2441 << TmpAlphaOFFHistoName << endl;
2442
2443 AlphaOFFAfterCuts = (TH1F*) rootfile.Get(TmpAlphaOFFHistoName);
2444
2445
2446
2447
2448
2449 // 2) Computing Nex and Signficance vs alpha and filling histograms
2450 // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF
2451
2452
2453 *fLog << endl << endl
2454 << "Starting loop to compute Nex and Significance vs alphasig " << endl
2455 << "for TRAIN sample and "
2456 << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]
2457 << " :" << endl;
2458
2459
2460 for (Int_t j = 1; j < NexSignVSAlphaXBinsVectorDimension; j++)
2461 {
2462
2463
2464 alphasig = XBinsVector[j];
2465
2466 MHFindSignificanceONOFF findsig;
2467 findsig.SetRebin(RebinHistograms);
2468 findsig.SetReduceDegree(ReduceDegree);
2469 findsig.SetUseFittedQuantities(fUseFittedQuantities);
2470
2471 // Dummy name for psfile
2472 TString psfilename = ("DummyPs");
2473
2474 findsig.FindSigmaONOFF(AlphaONAfterCuts, AlphaOFFAfterCuts,
2475 NormalizationFactor,
2476 alphabkgmin, alphabkgmax,
2477 fDegreeON, fDegreeOFF,
2478 alphasig,
2479 drawpoly, fitgauss, print,
2480 saveplots, psfilename);
2481
2482
2483 Significance = double(findsig.GetSignificance());
2484
2485 // so far there is not a variable that contains
2486 // NexONOFF or NexONOFFFitted (in class MHFindSignificanceONOFF)
2487 // according to the value of the variable fUseFittedQuantities.
2488 // As soon as I have some time I will implement it (with the
2489 // corresponging GET member function) and will change
2490 // MFindSignificanceONOFF accordingly to make use of it.
2491
2492 // For the time being, I will survive with an "if statement".
2493
2494
2495 if(fUseFittedQuantities)
2496 {
2497 Nex = double(findsig.GetNexONOFFFitted());
2498 }
2499 else
2500 {
2501 Nex = double(findsig.GetNexONOFF());
2502 }
2503
2504
2505
2506 *fLog << endl << "Cut in alpha, Nex and significance : "
2507 << alphasig << ", " << Nex << ", " << Significance << endl << endl;
2508
2509
2510 // updating histograms
2511
2512 tmpdouble = alphasig - AlphaBinWidth/2.0;
2513
2514 NexVSAlphaSigTrainVector[i] -> Fill(tmpdouble, Nex);
2515 SignificanceVSAlphaSigTrainVector[i] -> Fill(tmpdouble, Significance);
2516
2517
2518
2519
2520
2521
2522
2523 }
2524
2525
2526 // Dynamic memory allocated in this loop is released
2527
2528
2529 *fLog << "Memory allocated by temporal alpha histograms "
2530 << "will be released" << endl;
2531
2532 delete AlphaONAfterCuts;
2533 delete AlphaOFFAfterCuts;
2534
2535 }
2536
2537
2538 }
2539
2540
2541
2542 // ***************************************************
2543 // HISTOGRAMS FOR TEST SAMPLE ARE COMPUTED
2544 // ***************************************************
2545
2546
2547
2548
2549 if(fTestParameters)
2550 {// computation for TEST sample
2551
2552 *fLog << endl << endl
2553 << "**************************************************************************" << endl
2554 << "Computing Nex and significance vs cut in alpha (alphasig) for TEST sample" << endl
2555 << "**************************************************************************" << endl
2556 << endl << endl;
2557
2558 // Retrieving NormFactorTestHisto from root file.
2559
2560 *fLog << "Getting histogram with name " << NormFactorTestHistoName << endl;
2561
2562 NormFactorTestHisto = (TH1F*) rootfile.Get(NormFactorTestHistoName);
2563
2564
2565 // Check that bins in this histo correspond with the
2566 // CosTheta intervals defined by vector fCosThetaRangeVector
2567
2568 ////////////////// START CHECK //////////////////////////////////////
2569
2570 if (NormFactorTestHisto -> GetNbinsX() != NThetaBins)
2571 {
2572 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()"
2573 << endl
2574 << "Number of theta bins defined by vector fCosThetaRangeVector ("
2575 << NThetaBins << ") does not match with the bins in histogram "
2576 << NormFactorTestHistoName << "("
2577 << NormFactorTestHisto -> GetNbinsX() << ")" << endl
2578 << "Function execution will be ABORTED." << endl;
2579 return kFALSE;
2580 }
2581
2582 // histo test
2583 /*
2584 cout << "NORMFACTORTRAINHIST Test: BinCenter and Value" << endl;
2585 for (Int_t k = 0; k < NThetaBins; k++)
2586 {
2587
2588 cout << NormFactorTestHisto -> GetBinCenter(k+1)
2589 << "; " << NormFactorTestHisto -> GetBinContent(k+1)<< endl;
2590 }
2591 */
2592
2593 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
2594 {
2595 tmpdouble = NormFactorTestHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]);
2596 if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble ||
2597 fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble)
2598 {
2599 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()"
2600 << endl
2601 << "Bins defined by vector fCosThetaRangeVector "
2602 << "do not match with the bins of histogram "
2603 << NormFactorTestHistoName << endl
2604 << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1]
2605 << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl
2606 << "NormFactorHist Bin: " << tmpdouble << endl
2607 << "Function execution will be ABORTED." << endl;
2608
2609 return kFALSE;
2610 }
2611
2612 }
2613
2614 ////////////////// END CHECK //////////////////////////////////////
2615
2616
2617
2618
2619 // Check that the theta range string (needed to compute names for alpha histograms)
2620 // is well defined
2621
2622 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
2623 {
2624 if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
2625 {
2626 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()"
2627 << endl
2628 << "Component " << SuccessfulThetaBinsVector[i]-1
2629 << " of fThetaRangeStringVector is EMPTY"
2630 << endl
2631 << "Function execution will be ABORTED." << endl;
2632
2633 return kFALSE;
2634
2635
2636 }
2637 }
2638
2639
2640 // *************************************************************
2641 // *************************************************************
2642
2643 // Normfactors and successfultheta bins are ok.
2644 // I can now compute the Nex and Sigma for test sample (ON-OFF)
2645
2646 // *************************************************************
2647 // *************************************************************
2648
2649 // Loop over the several successful theta bins doing the following
2650
2651 // 0) Initializing histograms Nex and Significance (with names!!)
2652 // 1) Retrieving alpha histograms (ON-OFF), and normalization factor
2653 // 2) Computing Nex and Signficance vs alpha and filling histograms
2654
2655
2656 *fLog << endl
2657 << "Starting loop over the successfully optimized theta ranges " << endl
2658 << endl;
2659
2660
2661 Double_t NormalizationFactor = 0.0;
2662
2663 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
2664 {
2665
2666 *fLog << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] << endl << endl;
2667
2668 // 0) Initializing histograms Nex and Significance (with names!!)
2669
2670 TString NexVSAlphaNameTmp = (NexVSAlphaName);
2671 TString SignificanceVSAlphaNameTmp = (SignificanceVSAlphaName);
2672
2673 NexVSAlphaNameTmp += "Test";
2674 NexVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
2675
2676 SignificanceVSAlphaNameTmp += "Test";
2677 SignificanceVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
2678
2679
2680 // TEMP
2681 // cout << NexVSAlphaNameTmp << endl;
2682 // cout << SignificanceVSAlphaNameTmp << endl;
2683 // ENDTEMP
2684
2685
2686
2687
2688 NexVSAlphaSigTestVector[i] = new TH1F (NexVSAlphaNameTmp.Data(),
2689 NexVSAlphaNameTmp.Data(),
2690 NexSignVSAlphaXBinsVectorDimension - 1,
2691 XBinsVector);
2692
2693 // It is important to unlink the the created histograms from
2694 // the current directory (defined now by "rootfile"),
2695 // so that we can do whatever we want to
2696 // with them regardless of such directory is opened or closed
2697 // The reference of the histogram is removed from the current
2698 // directory by means of the member function "SetDirectory(0)"
2699
2700
2701 NexVSAlphaSigTestVector[i] -> SetDirectory(0);
2702
2703
2704 SignificanceVSAlphaSigTestVector[i] = new TH1F (SignificanceVSAlphaNameTmp.Data(),
2705 SignificanceVSAlphaNameTmp.Data(),
2706 NexSignVSAlphaXBinsVectorDimension - 1,
2707 XBinsVector);
2708
2709 SignificanceVSAlphaSigTestVector[i] -> SetDirectory(0);
2710
2711
2712
2713 // axis are defined
2714
2715 NexVSAlphaSigTestVector[i] -> SetXTitle(x_axis_label.Data());
2716 NexVSAlphaSigTestVector[i] -> SetYTitle(y_axis_label_nex.Data());
2717
2718 SignificanceVSAlphaSigTestVector[i] ->SetXTitle(x_axis_label.Data());
2719 SignificanceVSAlphaSigTestVector[i] ->SetYTitle(y_axis_label_significance.Data());
2720
2721
2722 // 1) Retrieving alpha histograms (ON-OFF) and normalization factor
2723
2724 NormalizationFactor =
2725 NormFactorTestHisto -> GetBinContent(SuccessfulThetaBinsVector[i]);
2726
2727
2728 TString TmpAlphaHistoName = ("TestAlphaAfterCuts");
2729 TmpAlphaHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
2730
2731 *fLog << "Getting Histogram with name " << endl
2732 << TmpAlphaHistoName << endl;
2733
2734 AlphaONAfterCuts = (TH1F*) rootfile.Get(TmpAlphaHistoName);
2735
2736
2737 TString TmpAlphaOFFHistoName = ("TestAlphaOFFAfterCuts");
2738 TmpAlphaOFFHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
2739
2740 *fLog << "Getting Histogram with name " << endl
2741 << TmpAlphaOFFHistoName << endl;
2742
2743 AlphaOFFAfterCuts = (TH1F*) rootfile.Get(TmpAlphaOFFHistoName);
2744
2745
2746 // 2) Computing Nex and Signficance vs alpha and filling histograms
2747 // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF
2748
2749
2750 *fLog << endl << endl
2751 << "Starting loop to compute Nex and Significance vs alphasig " << endl
2752 << "for TRAIN sample and "
2753 << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]
2754 << " :" << endl;
2755
2756
2757 for (Int_t j = 1; j < NexSignVSAlphaXBinsVectorDimension; j++)
2758 {
2759
2760
2761 alphasig = XBinsVector[j];
2762
2763 MHFindSignificanceONOFF findsig;
2764 findsig.SetRebin(RebinHistograms);
2765 findsig.SetReduceDegree(ReduceDegree);
2766 findsig.SetUseFittedQuantities(fUseFittedQuantities);
2767
2768 // Dummy name for psfile
2769 TString psfilename = ("DummyPs");
2770
2771 findsig.FindSigmaONOFF(AlphaONAfterCuts, AlphaOFFAfterCuts,
2772 NormalizationFactor,
2773 alphabkgmin, alphabkgmax,
2774 fDegreeON, fDegreeOFF,
2775 alphasig,
2776 drawpoly, fitgauss, print,
2777 saveplots, psfilename);
2778
2779
2780
2781 Significance = double(findsig.GetSignificance());
2782
2783 // so far there is not a variable that contains
2784 // NexONOFF or NexONOFFFitted (in class MHFindSignificanceONOFF)
2785 // according to the value of the variable fUseFittedQuantities.
2786 // As soon as I have some time I will implement it (with the
2787 // corresponging GET member function) and will change
2788 // MFindSignificanceONOFF accordingly to make use of it.
2789
2790 // For the time being, I will survive with an "if statement".
2791
2792
2793 if(fUseFittedQuantities)
2794 {
2795 Nex = double(findsig.GetNexONOFFFitted());
2796 }
2797 else
2798 {
2799 Nex = double(findsig.GetNexONOFF());
2800 }
2801
2802
2803 *fLog << endl << "Cut in alpha, Nex and significance : "
2804 << alphasig << ", " << Nex << ", " << Significance << endl << endl;
2805
2806
2807 // updating histograms
2808
2809 tmpdouble = alphasig - AlphaBinWidth/2.0;
2810
2811 NexVSAlphaSigTestVector[i] -> Fill(tmpdouble, Nex);
2812 SignificanceVSAlphaSigTestVector[i] -> Fill(tmpdouble, Significance);
2813
2814
2815
2816
2817
2818
2819
2820 }
2821
2822
2823 // Dynamic memory allocated in this loop is released
2824
2825
2826 *fLog << "Memory allocated by temporal alpha histograms "
2827 << "will be released" << endl;
2828
2829 delete AlphaONAfterCuts;
2830 delete AlphaOFFAfterCuts;
2831
2832 }
2833
2834
2835
2836
2837
2838
2839
2840
2841 }
2842
2843 rootfile.Close();
2844
2845
2846
2847 // *************************************************************
2848 // Histograms Nex and Significance vs alpha are written into
2849 // root file fAlphaDistributionsRootFilename
2850
2851
2852
2853
2854 *fLog <<"Histograms containing Nex and significance vs alphasig "
2855 << " are saved into root file " << endl
2856 << fAlphaDistributionsRootFilename << endl;
2857
2858
2859 TFile rootfiletowrite (fAlphaDistributionsRootFilename, "UPDATE",
2860 "Alpha Distributions for several Theta bins");
2861
2862
2863 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
2864 {
2865
2866 if (fOptimizeParameters)
2867 {// Histograms for train sample are written to root file
2868
2869 if(NexVSAlphaSigTrainVector[i])
2870 {
2871 if(NexVSAlphaSigTrainVector[i]-> GetEntries() > 0.5)
2872 NexVSAlphaSigTrainVector[i]-> Write();
2873
2874 }
2875
2876 if(SignificanceVSAlphaSigTrainVector[i])
2877 {
2878 if(SignificanceVSAlphaSigTrainVector[i]-> GetEntries() > 0.5)
2879 SignificanceVSAlphaSigTrainVector[i]-> Write();
2880
2881 }
2882 }
2883
2884
2885
2886 if (fTestParameters)
2887 {// Histograms for test sample are written to root file
2888
2889 if(NexVSAlphaSigTestVector[i])
2890 {
2891 if(NexVSAlphaSigTestVector[i]-> GetEntries() > 0.5)
2892 NexVSAlphaSigTestVector[i]-> Write();
2893
2894 }
2895
2896 if(SignificanceVSAlphaSigTestVector[i])
2897 {
2898 if(SignificanceVSAlphaSigTestVector[i]-> GetEntries() > 0.5)
2899 SignificanceVSAlphaSigTestVector[i]-> Write();
2900
2901 }
2902 }
2903
2904
2905
2906 }
2907
2908
2909 rootfiletowrite.Close();
2910
2911
2912
2913
2914 // Dynamic memory allocated is released
2915
2916 delete XBinsVector;
2917
2918 // function finished successfully... hopefully...
2919 return kTRUE;
2920
2921}
2922
2923
2924
2925// Function that gets the histograms with the alpha distributions
2926// (for all the theta bins specified by fCosThetaRangeVector)
2927// stored in fAlphaDistributionsRootFilename, and combine them
2928// (correcting OFF histograms with the normalization factors stored
2929// in NormFactorTrain or NormFactorTest) to get one single
2930// Alpha distribution for ON and another one for OFF.
2931// Then these histograms are given as arguments to
2932// the function MHFindSignificanceONOFF::FindSigmaONOFF,
2933// (Object of this class is created) to compute the
2934// Overall Excess events and significance, that will be
2935// stored in variables fOverallNexTrain/Test and fOverallSigmaLiMaTrain/Test
2936
2937
2938
2939
2940Bool_t MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance(Bool_t CombineTrainData,
2941 Bool_t CombineTestData)
2942{
2943
2944
2945 *fLog << endl << endl
2946 << "****************************************************************************" << endl
2947 << "Combining number of excess events and significance for the " << endl
2948 << "selected bins in theta angle using the function " << endl
2949 << " MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance() " << endl
2950 << "****************************************************************************" << endl
2951 << endl;
2952
2953
2954
2955
2956
2957
2958 // First of all let's check that we have the proper the "drinks" for the party
2959
2960 if (fAlphaDistributionsRootFilename.IsNull())
2961 {
2962 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl
2963 << "Root file containing all alpha distributions "
2964 <<" (fAlphaDistributionsRootFilename) is not well defined." << endl
2965 << "Execution of the function is ABORTED" << endl;
2966 return kFALSE;
2967 }
2968
2969
2970 if (fCosThetaRangeVector.GetSize() < 2)
2971
2972 {
2973 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl
2974 << "Dimension of vector fCosThetaRangeVector is smaller than 2."
2975 << "So, vector fCosThetaRangeVector is NOT properly defined" << endl
2976 << "Function execution will be ABORTED" << endl;
2977
2978 return kFALSE;
2979 }
2980
2981
2982 TArrayI SuccessfulThetaBinsVector;
2983 // Vector containing the histogram index of the successfully
2984 // optimized theta bins.
2985 // It is filled with the contents of the histogram
2986 // fSuccessfulThetaBinsHist retrieved from the root
2987 // file fAlphaDistributionsRootFilename
2988 // Remember that histogram index START AT 1 !!!
2989
2990
2991
2992 // Histograms that will contain the alpha distributions
2993 // for the combined theta bins
2994
2995
2996 TH1F* CombinedAlphaTrainON = NULL;
2997 TString CombinedAlphaTrainONName = ("TrainSampleAlphaONCombinedThetaBins");
2998 TH1F* CombinedAlphaTrainOFF = NULL;
2999 TString CombinedAlphaTrainOFFName = ("TrainSampleAlphaOFFCombinedThetaBins");
3000
3001 TH1F* CombinedAlphaTestON = NULL;
3002 TString CombinedAlphaTestONName = ("TestSampleAlphaONCombinedThetaBins");
3003
3004 TH1F* CombinedAlphaTestOFF = NULL;
3005 TString CombinedAlphaTestOFFName = ("TestSampleAlphaONCombinedThetaBins");
3006
3007
3008
3009
3010
3011 Int_t NAlphaBins = 0; // Amount of bins of alpha histograms is expected to
3012 // be equal for all alpha histograms
3013
3014 // Vectors storing the computed error for each of the NAlphaBins
3015 // Only 2 vectors are needed for OFF Train and OFF Test sample, whose
3016 // alpha distributions are corrected with normalization factors.
3017
3018 // The ON sample is just added (no normalization factors are applied)
3019 // and therefore, the error of the bin is the assumed by default
3020 // (i.e Sqrt(BinContent)
3021
3022 TArrayD ErrorInCombinedAlphaTrainOFF;
3023 TArrayD ErrorInCombinedAlphaTestOFF;
3024
3025 TArrayD CombinedAlphaTrainOFFWithoutNormalization;
3026 TArrayD CombinedAlphaTestOFFWithoutNormalization;
3027
3028
3029
3030 const Double_t SmallQuantity = 0.01; // This quantity will be used as the
3031 // maximum difference between quantities (double) that are expected to be equal,
3032 // as the center of the bins in the several alpha distributions
3033
3034 TH1F* NormFactorTrainHisto;
3035 TString NormFactorTrainHistoName = ("NormFactorTrainHist");
3036 TH1F* NormFactorTestHisto;
3037 TString NormFactorTestHistoName = ("NormFactorTestHist");
3038 TString SuccessfulThetaBinsHistoName = ("SuccessfulThetaBinsHist");
3039
3040
3041 TH1F* SuccessfulThetaBinsHisto;
3042
3043 TH1F* TmpTH1FHisto;
3044
3045
3046 // Retrieve number of theta bins that have to be combined
3047 // from TArrayD fCosThetaRangeVector
3048
3049 Int_t tmpdim = fCosThetaRangeVector.GetSize() - 1;
3050 const Int_t NThetaBins = tmpdim;
3051
3052 // Variables to store the mean normalization factor of the
3053 // combined OFF sample.
3054 // MeanNormFactor = Sum (NormFactor_i*Noff_i)/Sum (Noff_i)
3055
3056
3057 Double_t MeanNormFactorTrain = 0.0;
3058 Double_t MeanNormFactorTest = 0.0;
3059
3060 Double_t TotalNumberOFFTrain = 0.0;
3061 Double_t TotalNumberOFFTest = 0.0;
3062
3063
3064
3065 Double_t EventCounter = 0.0;// variable used to count events used to fill histograms
3066 // it is double (and not Int) because the contents of an histogram bin might not
3067 // be an integer value... due to normalization factors !!
3068 Int_t tmpint = 0;
3069 Double_t tmpdouble = 0.0; // tmp variable double
3070 Double_t tmpdouble2 = 0.0; // tmp variable double
3071 Double_t tmperror = 0.0; // quantity used to store
3072 // temporaly the bin error in OFF histogram
3073 Double_t tmpnormfactor = 0.0; // variable used to store normalization factor
3074 // of the theta bin i (used for TRAIN and TEST sample)
3075
3076
3077
3078
3079 // Open root file
3080
3081
3082
3083
3084 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl
3085 << "Opening root file " << fAlphaDistributionsRootFilename << endl;
3086
3087
3088
3089 TFile rootfile (fAlphaDistributionsRootFilename, "READ");
3090
3091 // Retrieving SuccessfulThetaBinsHisto from root file.
3092
3093
3094 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl
3095 << "Retrieving SuccessfulThetaBinsHisto from root file and filling vector "
3096 << " SuccessfulThetaBinsVector" << endl;
3097
3098
3099 cout << "Getting histogram " << SuccessfulThetaBinsHistoName << endl;
3100
3101 SuccessfulThetaBinsHisto = (TH1F*) rootfile.Get(SuccessfulThetaBinsHistoName);
3102
3103
3104
3105 // Check that bins in this histo correspond with the
3106 // CosTheta intervals defined by vector fCosThetaRangeVector
3107
3108
3109 if (SuccessfulThetaBinsHisto -> GetNbinsX() != NThetaBins)
3110 {
3111 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3112 << endl
3113 << "Number of theta bins defined by vector fCosThetaRangeVector ("
3114 << NThetaBins << ") does not match with the bins in histogram "
3115 << SuccessfulThetaBinsHistoName << "("
3116 << SuccessfulThetaBinsHisto -> GetNbinsX() << ")" << endl
3117 << "Function execution will be ABORTED." << endl;
3118 return kFALSE;
3119 }
3120
3121
3122 // Filling vector SuccessfulThetaBinsVector
3123
3124 tmpint = 0;
3125 for (Int_t i = 0; i < NThetaBins; i++)
3126 {
3127 if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5)
3128 {
3129 tmpint++;
3130 }
3131 }
3132
3133 SuccessfulThetaBinsVector.Set(tmpint);
3134
3135 tmpint = 0;
3136 for (Int_t i = 0; i < NThetaBins; i++)
3137 {
3138 if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5)
3139 {
3140 if(tmpint < SuccessfulThetaBinsVector.GetSize())
3141 {
3142 SuccessfulThetaBinsVector[tmpint] = i+1;
3143 tmpint++;
3144 }
3145 else
3146 {
3147 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3148 << endl
3149 << "Problem when filling vector 'SuccessfulThetaBinsVector'"
3150 << endl
3151 << "Function execution will be ABORTED." << endl;
3152 return kFALSE;
3153 }
3154
3155
3156 }
3157 }
3158
3159
3160
3161
3162
3163 // ********* HISTOGRAMS FOR TRAIN SAMPLE WILL BE FILLED *********///
3164
3165 if (CombineTrainData)
3166 {
3167
3168
3169
3170 *fLog << endl
3171 << "****************************************************************************" << endl
3172 << "Combining data for the TRAIN sample " << endl
3173 << "****************************************************************************" << endl
3174 << endl;
3175
3176 // Retrieving NormFactorTrainHisto from root file.
3177
3178 NormFactorTrainHisto = (TH1F*) rootfile.Get(NormFactorTrainHistoName);
3179 // NormFactorTrainHisto-> DrawCopy();
3180
3181 // Check that bins in this histo correspond with the
3182 // CosTheta intervals defined by vector fCosThetaRangeVector
3183
3184
3185 if (NormFactorTrainHisto -> GetNbinsX() != NThetaBins)
3186 {
3187 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3188 << endl
3189 << "Number of theta bins defined by vector fCosThetaRangeVector ("
3190 << NThetaBins << ") does not match with the bins in histogram "
3191 << NormFactorTrainHistoName << "("
3192 << NormFactorTrainHisto -> GetNbinsX() << ")" << endl
3193 << "Function execution will be ABORTED." << endl;
3194 return kFALSE;
3195 }
3196
3197 // histo test
3198 /*
3199 cout << "NORMFACTORTRAINHIST TEST: BinCenter and Value" << endl;
3200 for (Int_t k = 0; k < NThetaBins; k++)
3201 {
3202
3203 cout << NormFactorTrainHisto -> GetBinCenter(k+1)
3204 << "; " << NormFactorTrainHisto -> GetBinContent(k+1)<< endl;
3205 }
3206 */
3207
3208 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
3209 {
3210 tmpdouble = NormFactorTrainHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]);
3211 if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble ||
3212 fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble)
3213 {
3214 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3215 << endl
3216 << "Bins defined by vector fCosThetaRangeVector "
3217 << "do not match with the bins of histogram "
3218 << NormFactorTrainHistoName << endl
3219 << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1]
3220 << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl
3221 << "NormFactorHist Bin: " << tmpdouble << endl
3222 << "Function execution will be ABORTED." << endl;
3223
3224 return kFALSE;
3225 }
3226
3227 }
3228
3229
3230 // Let's summ up all ON alpha distributions
3231
3232 EventCounter = 0;
3233
3234 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
3235 {
3236 if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
3237 {
3238 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3239 << endl
3240 << "Component " << SuccessfulThetaBinsVector[i]-1
3241 << " of fThetaRangeStringVector is EMPTY"
3242 << endl
3243 << "Function execution will be ABORTED." << endl;
3244 return kFALSE;
3245
3246
3247 }
3248
3249
3250
3251 TString TmpHistoName = ("TrainAlphaAfterCuts");
3252 TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
3253
3254 *fLog << "Getting Histogram with name " << endl
3255 << TmpHistoName << endl;
3256
3257 if (i == 0)
3258 {
3259 CombinedAlphaTrainON = (TH1F*) rootfile.Get(TmpHistoName);
3260 NAlphaBins = CombinedAlphaTrainON -> GetNbinsX();
3261
3262 // Name of histogram is set
3263 CombinedAlphaTrainON -> SetName(CombinedAlphaTrainONName);
3264
3265 // update EventCounter
3266 EventCounter += CombinedAlphaTrainON -> GetEntries();
3267
3268 // tmp
3269 cout << "NAlphaBins in histo CombinedAlphaTrainON: "
3270 << NAlphaBins << endl;
3271 // endtmp
3272 }
3273 else
3274 {
3275 TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName);
3276
3277 // update EventCounter
3278 EventCounter += TmpTH1FHisto -> GetEntries();
3279
3280 for (Int_t j = 1; j <= NAlphaBins; j++)
3281 {
3282 // At some time alpha bins might not be the same for
3283 // the different alpha distributions. That's why
3284 // I will check it before summing the alpha histo contents
3285 tmpdouble = CombinedAlphaTrainON->GetBinCenter(j) -
3286 TmpTH1FHisto->GetBinCenter(j);
3287 tmpdouble = TMath::Abs(tmpdouble);
3288
3289
3290 if (tmpdouble > SmallQuantity)
3291 {
3292 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3293 << endl
3294 << "Bins among the several alpha ON histograms "
3295 << "for TRAIN sample do not match"
3296 << endl
3297 << "Function execution will be ABORTED." << endl;
3298 return kFALSE;
3299 }
3300
3301 tmpdouble = CombinedAlphaTrainON->GetBinContent(j) +
3302 TmpTH1FHisto->GetBinContent(j);
3303
3304
3305 CombinedAlphaTrainON-> SetBinContent(j,tmpdouble);
3306
3307 }
3308
3309 // Dynamic memory allocated for TmpTH1FHisto is released
3310
3311
3312 //tmp
3313 cout << "Memory allocated by temporal histogram TmpTH1FHisto "
3314 << "will be released" << endl;
3315 //endtmp
3316 delete TmpTH1FHisto;
3317
3318 }
3319
3320
3321 }
3322
3323 // Entries of histogram CombinedAlphaTrainON is set to
3324 // the events counted in EventCounter
3325
3326 tmpint = int (EventCounter);
3327 CombinedAlphaTrainON->SetEntries(tmpint);
3328
3329 // Let's summ up all OFF alpha distributions
3330
3331 EventCounter = 0;
3332
3333 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
3334 {
3335 if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
3336 {
3337 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3338 << endl
3339 << "Component " << SuccessfulThetaBinsVector[i]-1
3340 << " of fThetaRangeStringVector is EMPTY"
3341 << endl
3342 << "Function execution will be ABORTED." << endl;
3343 return kFALSE;
3344
3345
3346 }
3347
3348 TString TmpHistoName = ("TrainAlphaOFFAfterCuts");
3349 TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
3350
3351 // Normalization constant for this theta bin is stored now in tmpnormfactor
3352 tmpnormfactor = NormFactorTrainHisto -> GetBinContent(SuccessfulThetaBinsVector[i]);
3353
3354 // temp
3355 cout << "Normalization factor for bin "
3356 << SuccessfulThetaBinsVector[i] << " is "
3357 << tmpnormfactor << endl;
3358
3359 // endtemp
3360
3361 if (i == 0)
3362 {
3363 CombinedAlphaTrainOFF = (TH1F*) rootfile.Get(TmpHistoName);
3364 NAlphaBins = CombinedAlphaTrainOFF -> GetNbinsX();
3365
3366
3367 // Name of histogram is set
3368 CombinedAlphaTrainOFF -> SetName(CombinedAlphaTrainOFFName);
3369
3370 // Event counter is updated taking into account the
3371 // normalization factor
3372
3373 EventCounter += tmpnormfactor * CombinedAlphaTrainOFF -> GetEntries();
3374
3375 // Contribution to total number of OFF events and mean normfactor
3376 // of this theta bin is computed
3377
3378 TotalNumberOFFTrain = CombinedAlphaTrainOFF -> GetEntries();
3379 MeanNormFactorTrain = tmpnormfactor *
3380 CombinedAlphaTrainOFF -> GetEntries();
3381
3382
3383
3384 // ****
3385 // Dimension of ErrorInCombinedAlphaTrainOFF is set
3386 ErrorInCombinedAlphaTrainOFF.Set(NAlphaBins);
3387
3388 CombinedAlphaTrainOFFWithoutNormalization.Set(NAlphaBins);
3389
3390 // Histogram now is normalized (to the respective ON data)
3391 // by using the normalization constants stored in
3392 // histogram NormFactorTrainHisto
3393
3394 // Error of the normalized bins is computed and stored in
3395 // vector ErrorInCombinedAlphaTrainOFF
3396
3397 for (Int_t j = 1; j <= NAlphaBins; j++)
3398 {
3399
3400 // Number of events (without normalization) are
3401 // stored in vector CombinedAlphaTrainOFFWithoutNormalization
3402
3403 CombinedAlphaTrainOFFWithoutNormalization[j-1] =
3404 CombinedAlphaTrainOFF -> GetBinContent(j);
3405
3406 // Bin content is set
3407 tmpdouble2 = tmpnormfactor *
3408 CombinedAlphaTrainOFF -> GetBinContent(j);
3409
3410 CombinedAlphaTrainOFF -> SetBinContent(j,tmpdouble2);
3411
3412 // (Bin Error)^2 is computed and stored in
3413 // ErrorInCombinedAlphaTrainOFF
3414 // Bin error = Sqrt(BinContent) X Normalization Factor
3415
3416 tmperror = tmpdouble2*tmpnormfactor;
3417 ErrorInCombinedAlphaTrainOFF[j-1] = tmperror;
3418
3419 }
3420
3421 }
3422 else
3423 {
3424 TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName);
3425
3426
3427 // Event counter is updated taking into account the
3428 // normalization factor
3429
3430 EventCounter += tmpnormfactor * TmpTH1FHisto -> GetEntries();
3431
3432 // Contribution to total number of OFF events and mean normfactor
3433 // of this theta bin is computed
3434
3435 TotalNumberOFFTrain += TmpTH1FHisto -> GetEntries();
3436 MeanNormFactorTrain += tmpnormfactor *
3437 TmpTH1FHisto -> GetEntries();
3438
3439
3440
3441 for (Int_t j = 1; j <= NAlphaBins; j++)
3442 {
3443 // At some time alpha bins might not be the same for
3444 // the different alpha distributions. That's why
3445 // I will check it before summing the alpha histo contents
3446 tmpdouble2 = CombinedAlphaTrainOFF->GetBinCenter(j) -
3447 TmpTH1FHisto->GetBinCenter(j);
3448 tmpdouble2 = TMath::Abs(tmpdouble2);
3449
3450 if (tmpdouble2 > SmallQuantity)
3451 {
3452 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3453 << endl
3454 << "Bins among the several alpha OFF histograms "
3455 << "for TRAIN sample do not match"
3456 << endl
3457 << "Function execution will be ABORTED." << endl;
3458 return kFALSE;
3459 }
3460
3461
3462 // Number of events (without normalization) are
3463 // added in vector CombinedAlphaTrainOFFWithoutNormalization
3464
3465 CombinedAlphaTrainOFFWithoutNormalization[j-1] +=
3466 TmpTH1FHisto -> GetBinContent(j);
3467
3468 // Histogram bin contents must be normalized (to the respective ON data)
3469 // by using the normalization constants stored in
3470 // histogram NormFactorTrainHisto before being added
3471 // to CombinedAlphaTrainOFF
3472
3473
3474 tmpdouble2 = CombinedAlphaTrainOFF->GetBinContent(j) +
3475 (TmpTH1FHisto->GetBinContent(j) * tmpnormfactor);
3476
3477
3478 CombinedAlphaTrainOFF-> SetBinContent(j,tmpdouble2);
3479
3480
3481 // (Bin Error)^2 is computed and added to
3482 // ErrorInCombinedAlphaTrainOFF
3483 // Bin error = Sqrt(BinContent) X Normalization Factor
3484
3485 tmperror = TmpTH1FHisto->GetBinContent(j) * tmpnormfactor;
3486 ErrorInCombinedAlphaTrainOFF[j-1] += tmperror;
3487
3488 }
3489
3490 // Dynamic memory allocated for TmpTH1FHisto is released
3491
3492 delete TmpTH1FHisto;
3493
3494 }
3495
3496
3497
3498
3499
3500 }
3501
3502
3503
3504
3505 // Mean Normalization factor is computed for the combined OFF histogram.
3506 // This factor will be used when computing the significance via
3507 // MHFindSignificance::FindSigmaONOFF().
3508
3509 // The bin contents (and errors) of teh combined OFF TRAIN histo will be
3510 // normalized (divided) by teh mean norm constant, so that
3511 // the histogram for OFF data given as argument in FindSigmaONOFF
3512 // is equivalent to a non-normalized histogram.
3513
3514 MeanNormFactorTrain = MeanNormFactorTrain/TotalNumberOFFTrain;
3515
3516
3517
3518 // Sqrt(Contents) is performed in ErrorInCombinedAlphaTrainOFF,
3519 // in order to compute the real error in the bin content of histogram
3520 // CombinedAlphaTrainOFF.
3521
3522 // Then error and contents (normalized with mean normfctor) are
3523 // set to histogram bin
3524
3525 for (Int_t j = 1; j <= NAlphaBins; j++)
3526 {
3527 tmpdouble2 =
3528 (CombinedAlphaTrainOFF -> GetBinContent(j))/MeanNormFactorTrain;
3529
3530 CombinedAlphaTrainOFF -> SetBinContent(j,tmpdouble2);
3531
3532 ErrorInCombinedAlphaTrainOFF[j-1] =
3533 TMath::Sqrt(ErrorInCombinedAlphaTrainOFF[j-1])/MeanNormFactorTrain;
3534 // Proper error treatment when adding quantities to histogram bin
3535
3536 CombinedAlphaTrainOFF -> SetBinError(j,ErrorInCombinedAlphaTrainOFF[j-1]);
3537
3538 // ****** TMP **************
3539 /*
3540 // Error computation assuming a given (no normalized) bin content
3541
3542 tmpdouble = CombinedAlphaTrainOFF -> GetBinContent(j);
3543
3544 if (CombinedAlphaTrainOFFWithoutNormalization[j-1] < SmallQuantity)
3545 {tmpnormfactor = 0;}
3546 else
3547 {tmpnormfactor = tmpdouble/CombinedAlphaTrainOFFWithoutNormalization[j-1];}
3548
3549 tmperror = TMath::Sqrt(CombinedAlphaTrainOFFWithoutNormalization[j-1]) *
3550 tmpnormfactor;
3551
3552 CombinedAlphaTrainOFF -> SetBinError(j,tmperror);
3553
3554 */
3555 // ****** END TMP **************
3556
3557 // tmp
3558 /*
3559 cout << CombinedAlphaTrainOFF -> GetBinContent(j) << " +/- "
3560 << CombinedAlphaTrainOFF -> GetBinError(j)
3561 << " Number Events without normalization: "
3562 << CombinedAlphaTrainOFFWithoutNormalization[j-1] << endl;
3563 */
3564 // endtmp
3565
3566 }
3567
3568
3569 // Entries of histogram CombinedAlphaTrainON is set to
3570 // the events counted during histogram filling
3571
3572 EventCounter = EventCounter/MeanNormFactorTrain;
3573 tmpint = int (EventCounter);
3574 CombinedAlphaTrainOFF->SetEntries(tmpint);
3575
3576 // Check that the number of events in the combined normalized histogram is
3577 // that of the initial sample
3578
3579 /*
3580 cout << "SILLY INFO FOR TRAIN SAMPLE: Total number of events Before nomralization and "
3581 << " re-re normalized" << endl
3582 << TotalNumberOFFTrain << " - " << EventCounter << endl;
3583
3584 */
3585 }
3586
3587 // ********* COMBINED HISTOGRAMS FOR TRAIN SAMPLE ALREADY FILLED ********* ///
3588
3589
3590
3591
3592 // ********* HISTOGRAMS FOR TEST SAMPLE WILL BE FILLED ********* ///
3593
3594
3595 if (CombineTestData)
3596 {
3597
3598 *fLog << endl
3599 << "****************************************************************************" << endl
3600 << "Combining data for the TEST sample " << endl
3601 << "****************************************************************************" << endl
3602 << endl;
3603
3604
3605 // Retrieving NormFactorTestHisto from root file.
3606
3607 NormFactorTestHisto = (TH1F*) rootfile.Get(NormFactorTestHistoName);
3608 // NormFactorTestHisto-> DrawCopy();
3609
3610 // Check that bins in this histo correspond with the
3611 // CosTheta intervals defined by vector fCosThetaRangeVector
3612
3613
3614 if (NormFactorTestHisto -> GetNbinsX() != NThetaBins)
3615 {
3616 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3617 << endl
3618 << "Number of theta bins defined by vector fCosThetaRangeVector ("
3619 << NThetaBins << ") does not match with the bins in histogram "
3620 << NormFactorTestHistoName << "("
3621 << NormFactorTestHisto -> GetNbinsX() << ")" << endl
3622 << "Function execution will be ABORTED." << endl;
3623 return kFALSE;
3624 }
3625
3626 // histo test
3627 /*
3628 cout << "NORMFACTORTRAINHIST TEST: BinCenter and Value" << endl;
3629 for (Int_t k = 0; k < NThetaBins; k++)
3630 {
3631
3632 cout << NormFactorTestHisto -> GetBinCenter(k+1)
3633 << "; " << NormFactorTestHisto -> GetBinContent(k+1)<< endl;
3634 }
3635 */
3636
3637
3638
3639 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
3640 {
3641 tmpdouble = NormFactorTestHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]);
3642 if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble ||
3643 fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble)
3644 {
3645 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3646 << endl
3647 << "Bins defined by vector fCosThetaRangeVector "
3648 << "do not match with the bins of histogram "
3649 << NormFactorTestHistoName << endl
3650 << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1]
3651 << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl
3652 << "NormFactorHist Bin: " << tmpdouble << endl
3653 << "Function execution will be ABORTED." << endl;
3654
3655 return kFALSE;
3656 }
3657
3658 }
3659
3660
3661 // Let's summ up all ON alpha distributions
3662
3663 // Event counter (counts events used in the histogram filling) is
3664 // initialized to zero
3665
3666 EventCounter = 0;
3667
3668 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
3669 {
3670 if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
3671 {
3672 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3673 << endl
3674 << "Component " << SuccessfulThetaBinsVector[i]-1
3675 << " of fThetaRangeStringVector is EMPTY"
3676 << endl
3677 << "Function execution will be ABORTED." << endl;
3678 return kFALSE;
3679
3680
3681 }
3682
3683
3684
3685 TString TmpHistoName = ("TestAlphaAfterCuts");
3686 TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
3687
3688 *fLog << "Getting Histogram with name " << endl
3689 << TmpHistoName << endl;
3690
3691 if (i == 0)
3692 {
3693 CombinedAlphaTestON = (TH1F*) rootfile.Get(TmpHistoName);
3694 NAlphaBins = CombinedAlphaTestON -> GetNbinsX();
3695
3696 // Name of histogram is set
3697 CombinedAlphaTestON -> SetName(CombinedAlphaTestONName);
3698
3699 // EventCounter is updated
3700
3701 EventCounter += CombinedAlphaTestON -> GetEntries();
3702
3703 // tmp
3704 cout << "NAlphaBins in histo CombinedAlphaTestON: "
3705 << NAlphaBins << endl;
3706 // endtmp
3707 }
3708 else
3709 {
3710 TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName);
3711
3712 // EventCounter is updated
3713
3714 EventCounter += TmpTH1FHisto -> GetEntries();
3715
3716 for (Int_t j = 1; j <= NAlphaBins; j++)
3717 {
3718 // At some time alpha bins might not be the same for
3719 // the different alpha distributions. That's why
3720 // I will check it before summing the alpha histo contents
3721 tmpdouble = CombinedAlphaTestON->GetBinCenter(j) -
3722 TmpTH1FHisto->GetBinCenter(j);
3723 tmpdouble = TMath::Abs(tmpdouble);
3724
3725
3726 if (tmpdouble > SmallQuantity)
3727 {
3728 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3729 << endl
3730 << "Bins among the several alpha ON histograms "
3731 << "for TEST sample do not match"
3732 << endl
3733 << "Function execution will be ABORTED." << endl;
3734 return kFALSE;
3735 }
3736
3737 tmpdouble = CombinedAlphaTestON->GetBinContent(j) +
3738 TmpTH1FHisto->GetBinContent(j);
3739
3740
3741 CombinedAlphaTestON-> SetBinContent(j,tmpdouble);
3742
3743 }
3744
3745 // Dynamic memory allocated for TmpTH1FHisto is released
3746
3747
3748 //tmp
3749 cout << "Memory allocated by temporal histogram TmpTH1FHisto "
3750 << "will be released" << endl;
3751 //endtmp
3752 delete TmpTH1FHisto;
3753
3754 }
3755
3756
3757
3758
3759 }
3760
3761 // Entries of histogram CombinedAlphaTestON is set to
3762 // the events counted during histogram filling
3763 tmpint = int (EventCounter);
3764 CombinedAlphaTestON->SetEntries(tmpint);
3765
3766
3767 // Let's summ up all OFF alpha distributions
3768
3769 EventCounter = 0;
3770
3771 for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
3772 {
3773 if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
3774 {
3775 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3776 << endl
3777 << "Component " << SuccessfulThetaBinsVector[i]-1
3778 << " of fThetaRangeStringVector is EMPTY"
3779 << endl
3780 << "Function execution will be ABORTED." << endl;
3781 return kFALSE;
3782
3783
3784 }
3785
3786 TString TmpHistoName = ("TestAlphaOFFAfterCuts");
3787 TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];
3788
3789 // Normalization constant for this theta bin is stored now in tmpnormfactor
3790 tmpnormfactor = NormFactorTestHisto -> GetBinContent(SuccessfulThetaBinsVector[i]);
3791
3792 // temp
3793 cout << "Normalization factor for bin "
3794 << SuccessfulThetaBinsVector[i] << " is "
3795 << tmpnormfactor << endl;
3796
3797 // endtemp
3798
3799 if (i == 0)
3800 {
3801 CombinedAlphaTestOFF = (TH1F*) rootfile.Get(TmpHistoName);
3802 NAlphaBins = CombinedAlphaTestOFF -> GetNbinsX();
3803
3804
3805 // Name of histogram is set
3806 CombinedAlphaTestOFF -> SetName(CombinedAlphaTestOFFName);
3807
3808 // EventCounter is updated taking into account
3809 // normalization factor
3810
3811 EventCounter += tmpnormfactor*CombinedAlphaTestOFF -> GetEntries();
3812
3813
3814 // Contribution to total number of OFF events and mean normfactor
3815 // of this theta bin is computed
3816
3817 TotalNumberOFFTest = CombinedAlphaTestOFF -> GetEntries();
3818 MeanNormFactorTest = tmpnormfactor *
3819 CombinedAlphaTestOFF -> GetEntries();
3820
3821
3822 // Dimension of ErrorInCombinedAlphaTrainOFF is set
3823 ErrorInCombinedAlphaTestOFF.Set(NAlphaBins);
3824
3825 CombinedAlphaTestOFFWithoutNormalization.Set(NAlphaBins);
3826
3827
3828 // Histogram now is normalized (to the respective ON data)
3829 // by using the normalization constants stored in
3830 // histogram NormFactorTestHisto
3831
3832
3833 for (Int_t j = 1; j <= NAlphaBins; j++)
3834 {
3835
3836 // Vector containing number of events without
3837 // normalization is filled
3838
3839 CombinedAlphaTestOFFWithoutNormalization[j-1] =
3840 CombinedAlphaTestOFF -> GetBinContent(j);
3841
3842
3843
3844 tmpdouble2 = tmpnormfactor *
3845 CombinedAlphaTestOFF -> GetBinContent(j);
3846
3847 CombinedAlphaTestOFF -> SetBinContent(j,tmpdouble2);
3848
3849 // (Bin Error)^2 is computed and stored in
3850 // ErrorInCombinedAlphaTestOFF
3851 // Bin error = Sqrt(BinContent) X Normalization Factor
3852
3853 tmperror = tmpdouble2*tmpnormfactor;
3854 ErrorInCombinedAlphaTestOFF[j-1] = tmperror;
3855
3856
3857 }
3858
3859 }
3860 else
3861 {
3862 TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName);
3863
3864
3865 // EventCounter is updated taking into account
3866 // normalization factor
3867
3868 EventCounter += tmpnormfactor* TmpTH1FHisto-> GetEntries();
3869
3870 // Contribution to total number of OFF events and mean normfactor
3871 // of this theta bin is computed
3872
3873
3874 TotalNumberOFFTest += TmpTH1FHisto -> GetEntries();
3875 MeanNormFactorTest += tmpnormfactor *
3876 TmpTH1FHisto -> GetEntries();
3877
3878
3879
3880 for (Int_t j = 1; j <= NAlphaBins; j++)
3881 {
3882 // At some time alpha bins might not be the same for
3883 // the different alpha distributions. That's why
3884 // I will check it before summing the alpha histo contents
3885 tmpdouble2 = CombinedAlphaTestOFF->GetBinCenter(j) -
3886 TmpTH1FHisto->GetBinCenter(j);
3887 tmpdouble2 = TMath::Abs(tmpdouble2);
3888
3889 if (tmpdouble2 > SmallQuantity)
3890 {
3891 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()"
3892 << endl
3893 << "Bins among the several alpha OFF histograms "
3894 << "for TRAIN sample do not match"
3895 << endl
3896 << "Function execution will be ABORTED." << endl;
3897 return kFALSE;
3898 }
3899
3900
3901 // Vector containing number of events without
3902 // normalization is updated
3903
3904 // Vector containing number of events without
3905 // normalization is filled
3906
3907 CombinedAlphaTestOFFWithoutNormalization[j-1] +=
3908 TmpTH1FHisto -> GetBinContent(j);
3909
3910
3911
3912
3913
3914
3915 // Histogram bin contents must be normalized (to the respective ON data)
3916 // by using the normalization constants stored in
3917 // histogram NormFactorTestHisto before being added
3918 // to CombinedAlphaTestOFF
3919
3920
3921 tmpdouble2 = CombinedAlphaTestOFF->GetBinContent(j) +
3922 (TmpTH1FHisto->GetBinContent(j) * tmpnormfactor);
3923
3924
3925 CombinedAlphaTestOFF-> SetBinContent(j,tmpdouble2);
3926
3927
3928 // (Bin Error)^2 is computed and added to
3929 // ErrorInCombinedAlphaTestOFF
3930 // Bin error = Sqrt(BinContent) X Normalization Factor
3931
3932 tmperror = TmpTH1FHisto->GetBinContent(j) * tmpnormfactor;
3933 ErrorInCombinedAlphaTestOFF[j-1] += tmperror;
3934
3935
3936 }
3937
3938 // Dynamic memory allocated for TmpTH1FHisto is released
3939
3940 delete TmpTH1FHisto;
3941
3942 }
3943
3944
3945
3946
3947
3948 }
3949
3950
3951
3952 // Mean Normalization factor is computed for the combined OFF histogram.
3953 // This factor will be used when computing the significance via
3954 // MHFindSignificance::FindSigmaONOFF().
3955
3956 // The bin contents (and errors) of teh combined OFF TEST histo will be
3957 // normalized (divided) by teh mean norm constant, so that
3958 // the histogram for OFF data given as argument in FindSigmaONOFF
3959 // is equivalent to a non-normalized histogram.
3960
3961 MeanNormFactorTest = MeanNormFactorTest/TotalNumberOFFTest;
3962
3963
3964
3965 // Sqrt(Contents) is performed in ErrorInCombinedAlphaTestOFF,
3966 // in order to compute the real error in the bin content of histogram
3967 // CombinedAlphaTestOFF.
3968
3969 // Then error and contents (normalized with mean normfctor) are
3970 // set to histogram bin
3971
3972
3973
3974
3975
3976 for (Int_t j = 1; j <= NAlphaBins; j++)
3977 {
3978
3979 tmpdouble2 =
3980 (CombinedAlphaTestOFF -> GetBinContent(j))/MeanNormFactorTest;
3981
3982 CombinedAlphaTestOFF -> SetBinContent(j,tmpdouble2);
3983
3984 ErrorInCombinedAlphaTestOFF[j-1] =
3985 TMath::Sqrt(ErrorInCombinedAlphaTestOFF[j-1])/MeanNormFactorTest;
3986
3987 // Proper error treatment for the bin contents
3988 CombinedAlphaTestOFF -> SetBinError(j,ErrorInCombinedAlphaTestOFF[j-1]);
3989
3990
3991 // ****** TMP **************
3992 /*
3993 // Error computation assuming a given (no normalized) bin content
3994
3995 tmpdouble = CombinedAlphaTestOFF -> GetBinContent(j);
3996
3997 if (CombinedAlphaTestOFFWithoutNormalization[j-1] < SmallQuantity)
3998 {tmpnormfactor = 0;}
3999 else
4000 {tmpnormfactor = tmpdouble/CombinedAlphaTestOFFWithoutNormalization[j-1];}
4001
4002 tmperror = (TMath::Sqrt(CombinedAlphaTestOFFWithoutNormalization[j-1])) * tmpnormfactor;
4003
4004 CombinedAlphaTestOFF -> SetBinError(j,tmperror);
4005
4006 */
4007
4008 // ****** ENDTMP **************
4009
4010 // tmp
4011 // cout << CombinedAlphaTestOFF -> GetBinContent(j) << " +/- "
4012 // << CombinedAlphaTestOFF -> GetBinError(j) << endl;
4013
4014 // endtmp
4015
4016 }
4017
4018
4019
4020 // Entries of histogram CombinedAlphaTestOFF is set to
4021 // the events counted during histogram filling
4022
4023 EventCounter = EventCounter/MeanNormFactorTest;
4024 tmpint = int(EventCounter);
4025 CombinedAlphaTestOFF->SetEntries(tmpint);
4026
4027
4028 // Check that the number of events in the combined normalized histogram is
4029 // that of the initial sample
4030
4031 /*
4032 cout << "SILLY INFO FOR TEST SAMPLE: Total number of events Before nomralization and "
4033 << " re-re normalized" << endl
4034 << TotalNumberOFFTest << " - " << EventCounter << endl;
4035
4036 */
4037
4038
4039 /*
4040 // Sumw2() test
4041
4042 cout << "******** TEST ********* " << endl;
4043 cout << "Sumw2() is executed: " << endl;
4044 CombinedAlphaTestOFF -> Sumw2();
4045 for (Int_t j = 1; j <= NAlphaBins; j++)
4046 {
4047 // tmp
4048 cout << CombinedAlphaTestOFF -> GetBinContent(j) << " +/- "
4049 << CombinedAlphaTestOFF -> GetBinError(j) << endl;
4050
4051 // endtmp
4052
4053 }
4054
4055 */
4056
4057
4058
4059 }
4060
4061
4062
4063 // ********* COMBINED HISTOGRAMS FOR TEST SAMPLE ALREADY FILLED ********* ///
4064
4065
4066 /*
4067 CombinedAlphaTrainON-> DrawCopy();
4068 gPad -> SaveAs("tmpONTrain.ps");
4069
4070 CombinedAlphaTrainOFF-> DrawCopy();
4071 gPad -> SaveAs("tmpOFFTrain.ps");
4072
4073
4074 CombinedAlphaTestON-> DrawCopy();
4075 gPad -> SaveAs("tmpONTest.ps");
4076
4077 CombinedAlphaTestOFF-> DrawCopy();
4078 gPad -> SaveAs("tmpOFFTest.ps");
4079
4080 */
4081
4082
4083 // **** COMPUTATION OF SIGNIFICANCES FOR COMBINED ALPHA HISTOGRAMS ****** ///
4084
4085 if (SuccessfulThetaBinsVector.GetSize() >= 1)
4086 {
4087
4088
4089 *fLog << endl
4090 << "**************************************************************************************" << endl
4091 << "Computation of excess events and singnificance for the combined alpha histograms " << endl
4092 << "**************************************************************************************" << endl
4093 << endl;
4094
4095 // There is, at least, ONE theta bin for which optimization of
4096 // supercuts was SUCCESSFUL, and thus, it is worth to
4097 // compute significance of signal
4098
4099
4100 // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF
4101
4102
4103
4104 // Maximum value of alpha below which signal is expected
4105 const Double_t alphasig = fAlphaSig;
4106
4107 // Minimum value of alpha for bkg region in ON data
4108 const Double_t alphabkgmin = fAlphaBkgMin;
4109
4110 // Maximum value of alpha for bkg region in ON data
4111 const Double_t alphabkgmax = fAlphaBkgMax;
4112
4113
4114 const Bool_t drawpoly = kTRUE;
4115 const Bool_t fitgauss = kTRUE;
4116 const Bool_t print = kTRUE;
4117
4118 const Bool_t saveplots = kTRUE; // Save final alpha plot for ON and OFF
4119 // (with Nex and significance) in Psfile
4120
4121
4122
4123
4124 MHFindSignificanceONOFF findsig;
4125 findsig.SetRebin(kTRUE);
4126 findsig.SetReduceDegree(kFALSE);
4127
4128
4129
4130 if (CombineTrainData)
4131 {
4132
4133 // Name of psfile where Alpha plot will be stored
4134
4135 TString psfilename = (fPathForFiles);
4136 psfilename += ("TRAINSampleCombinedThetaBins");
4137
4138 findsig.FindSigmaONOFF(CombinedAlphaTrainON, CombinedAlphaTrainOFF,
4139 MeanNormFactorTrain,
4140 alphabkgmin, alphabkgmax,
4141 fDegreeON, fDegreeOFF,
4142 alphasig,
4143 drawpoly, fitgauss, print,
4144 saveplots, psfilename);
4145
4146
4147 fOverallSigmaLiMaTrain = double(findsig.GetSignificance());
4148 fOverallNexTrain = double(findsig.GetNexONOFF());
4149
4150 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl
4151 << "After performing the combined analysis for TRAIN sample it was found "
4152 << "the following Nex and Significance: " << fOverallNexTrain
4153 << ", " << fOverallSigmaLiMaTrain << endl;
4154
4155
4156 }
4157
4158
4159 if (CombineTestData)
4160 {
4161
4162 // Name of psfile where Alpha plot will be stored
4163 TString psfilename = (fPathForFiles);
4164 psfilename += ("TESTSampleCombinedThetaBins");
4165
4166 findsig.FindSigmaONOFF(CombinedAlphaTestON, CombinedAlphaTestOFF,
4167 MeanNormFactorTest,
4168 alphabkgmin, alphabkgmax,
4169 fDegreeON, fDegreeOFF,
4170 alphasig,
4171 drawpoly, fitgauss, print,
4172 saveplots, psfilename);
4173
4174 fOverallSigmaLiMaTest = double(findsig.GetSignificance());
4175 fOverallNexTest = double(findsig.GetNexONOFF());
4176
4177 *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl
4178 << "After performing the combined analysis for TEST sample it was found "
4179 << "the following Nex and Significance: " << fOverallNexTest
4180 << ", " << fOverallSigmaLiMaTest << endl;
4181
4182 }
4183
4184
4185
4186 *fLog << endl
4187 << "**************************************************************************************" << endl
4188 << "Alpha distributions (Train/Test and ON/OFF) for the combined " << endl
4189 << " theta bins are stored in the root file " << endl
4190 << fAlphaDistributionsRootFilename << endl
4191 << "**************************************************************************************" << endl
4192 << endl;
4193
4194
4195
4196 TFile rootfiletowrite(fAlphaDistributionsRootFilename, "UPDATE");
4197
4198
4199
4200 // Histograms that will contain the Norm factors for the
4201 // combined TRAIN/TEST theta bins are created
4202
4203
4204 TString CombinedNormFactorTrainName = ("TrainSampleNormFactorCombinedThetaBins");
4205 TH1F CombinedNormFactorTrain (CombinedNormFactorTrainName,
4206 CombinedNormFactorTrainName,
4207 1, 0.0, 1.0);
4208
4209 TString CombinedNormFactorTestName = ("TestSampleNormFactorCombinedThetaBins");
4210 TH1F CombinedNormFactorTest (CombinedNormFactorTestName,
4211 CombinedNormFactorTestName,
4212 1, 0.0, 1.0);
4213
4214
4215 if (CombineTrainData)
4216 {
4217
4218 // Histogram that will contain the Norm factors for the
4219 // combined TRAIN theta bins are filled
4220
4221
4222 CombinedNormFactorTrain.Fill(0.5, MeanNormFactorTrain);
4223
4224
4225 // Train histos are written to root file
4226
4227 CombinedAlphaTrainON -> Write();
4228 CombinedAlphaTrainOFF -> Write();
4229 CombinedNormFactorTrain.Write();
4230
4231
4232
4233 }
4234
4235
4236 if (CombineTestData)
4237 {
4238
4239 // Histogram that will contain the Norm factors for the
4240 // combined TEST theta bins are filled
4241
4242
4243 CombinedNormFactorTest.Fill(0.5, MeanNormFactorTest);
4244
4245 // Test histos are written to root file
4246
4247 CombinedAlphaTestON -> Write();
4248 CombinedAlphaTestOFF -> Write();
4249 CombinedNormFactorTest.Write();
4250
4251
4252 }
4253
4254
4255 // Root file is closed
4256
4257 rootfiletowrite.Close();
4258
4259
4260 }
4261
4262
4263 return kTRUE;
4264}
4265
4266
4267
4268
4269
4270
4271
4272
Note: See TracBrowser for help on using the repository browser.