/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): David Paneque 12/2003 ! ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MFindSupercutsONOFFThetaLoop // // // // Class for optimizing the parameters of the supercuts // // Using ON and OFF data runing over data which are subdivided in // several Theta ranges. // // // // // ///////////////////////////////////////////////////////////////////////////// #include "MFindSupercutsONOFFThetaLoop.h" #include // fabs #include #include #include #include #include #include #include #include #include #include #include #include "MFindSupercutsONOFF.h" #include "MBinning.h" #include "MContinue.h" #include "MSupercuts.h" #include "MSupercutsCalcONOFF.h" #include "MDataElement.h" #include "MDataMember.h" #include #include "MEvtLoop.h" #include "MFCT1SelFinal.h" #include "MF.h" #include "MFEventSelector.h" #include "MFEventSelector2.h" #include "MFillH.h" //#include "MGeomCamCT1Daniel.h" //#include "MGeomCamCT1.h" #include "MGeomCamMagic.h" // only for magic #include "MFRandomSplit.h" #include "MH3.h" #include "MHCT1Supercuts.h" #include "MHFindSignificance.h" // To be removed at some point... #include "MHFindSignificanceONOFF.h" #include "MHMatrix.h" #include "MHOnSubtraction.h" #include "MDataValue.h" // #include "MDataString.h" #include "MLog.h" #include "MLogManip.h" #include "MMatrixLoop.h" #include "MMinuitInterface.h" #include "MParList.h" #include "MProgressBar.h" #include "MReadMarsFile.h" #include "MReadTree.h" #include "MTaskList.h" ClassImp(MFindSupercutsONOFF); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. // MFindSupercutsONOFFThetaLoop::MFindSupercutsONOFFThetaLoop(const char *name, const char *title) { fName = name ? name : "MFindSupercutsONOFF"; fTitle = title ? title : "Optimizer of the supercuts"; //--------------------------- // ********************************* // ** NOT TRUE ANY LONGER... ** // // Cuts (low and up) in variable ThetaOrig.fVal // The default is not cut, i.e. all values (0-1) are taken // fThetaOrig.fVal is measured in Radians; thus 1 = 57 degrees. // fThetaMin = 0.0; // fThetaMax = 1.0; // ********************************* // For MAGIC new (July 2004) Theta is measured in degrees fThetaMin = 0.0; fThetaMax = 90.0; fAlphaSig = 20; // By default, signal is expected in alpha<20 for CT1 // By default, bkg region is set to alpha range 30-90 for CT1 fAlphaBkgMin = 30; fAlphaBkgMax = 90; fActualCosThetaBinCenter = 0; fNormFactorTrainHist = NULL; //fNormFactorTrainHist -> SetDirectory(0); fNormFactorTestHist = NULL; //fNormFactorTestHist -> SetDirectory(0); fSigmaLiMaTrainHist = NULL; //fSigmaLiMaTrainHist -> SetDirectory(0); fSigmaLiMaTestHist = NULL; //fSigmaLiMaTestHist -> SetDirectory(0); fNexTrainHist = NULL; //fSigmaLiMaTestHist -> SetDirectory(0); fNexTestHist = NULL; //fNexTestHist -> SetDirectory(0); fSuccessfulThetaBinsHist = NULL; /* fNexVSAlphaSigTrain = NULL; fNexVSAlphaSigTest = NULL; fSignificanceVSAlphaSigTrain = NULL; fSignificanceVSAlphaSigTest = NULL; */ fNEvtsInTrainMatrixONHist = NULL; fNEvtsInTestMatrixONHist = NULL; fNEvtsInTrainMatrixOFFHist = NULL; fNEvtsInTestMatrixOFFHist = NULL; fTrainMatrixONFilenameVector = NULL; fTestMatrixONFilenameVector= NULL; fTrainMatrixOFFFilenameVector= NULL; fTestMatrixOFFFilenameVector= NULL; fOptSCParamFilenameVector = NULL; fThetaRangeStringVector = NULL; fPsFilename = NULL; fTuneNormFactor = kTRUE; // Norm factors will be corrected using the total amount of OFF events before cuts and the estimated excess events // fNormFactorTrain = fNormFactorTrain - Ngammas/EventsInTrainMatrixOFF fGammaEfficiency = 0.5; // Fraction of gammas that remain after cuts // Quantity that will have to be determined with MC, yet for the // time being I set it to 0.5 // boolean variables are set to kFALSE as default fReadMatricesFromFile = kFALSE; fOptimizeParameters = kFALSE; fTestParameters = kFALSE; // Fraction of Train/Test ON/OFF samples is set to 0.5 as default fWhichFractionTrain = 0.5; // number <= 1; specifying fraction of ON Train events fWhichFractionTest = 0.5; // number <= 1; specifying fraction of ON Test events fWhichFractionTrainOFF = 0.5; // number <= 1; specifying fraction of OFF Train events fWhichFractionTestOFF = 0.5; // number <= 1; specifying fraction of OFF Test events // use quantities computed from the fits // The variable allows the user to NOT use these quantities when there is not // enough statistics and fit not always is possible fUseFittedQuantities = kTRUE; // Boolean variable that controls wether the optimization of the // parameters (MMinuitInterface::CallMinuit(..) in function FindParams(..)) // takes place or not. kTRUE will skip such optimization. // This variable is useful to test the optmized parameters (previously found // and stored in root file) on the TRAIN sample. fSkipOptimization = kFALSE; // Boolean variable that allows the user to write the initial parameters // into the root file that will be used to store the optimum cuts. // If fUseInitialSCParams = kTRUE , parameters are written. // In this way, the initial SC parameters can be applied on the data (train/test) // The initial parameters are ONLY written to the root file if // there is NO SC params optimization, i.e., if variable // fSkipOptimization = kTRUE; // The default value is obviously kFALSE. fUseInitialSCParams = kTRUE; // fInitSCPar and fInitSCParSteps initialized to zero fInitSCPar.Set(0); fInitSCParSteps.Set(0); } // -------------------------------------------------------------------------- // // Default destructor. // MFindSupercutsONOFFThetaLoop::~MFindSupercutsONOFFThetaLoop() { *fLog << "Destructor of MFindSupercutsONOFFThetaLoop is called" << endl; fPsFilename = NULL; delete fNormFactorTrainHist; delete fNormFactorTestHist; delete fSigmaLiMaTrainHist; delete fSigmaLiMaTestHist; delete fNexTrainHist; delete fNexTestHist; delete fSuccessfulThetaBinsHist; delete fNEvtsInTrainMatrixONHist; delete fNEvtsInTestMatrixONHist; delete fNEvtsInTrainMatrixOFFHist; delete fNEvtsInTestMatrixOFFHist; delete [] fTrainMatrixONFilenameVector; delete [] fTestMatrixONFilenameVector; delete [] fTrainMatrixOFFFilenameVector; delete [] fTestMatrixOFFFilenameVector; delete [] fOptSCParamFilenameVector; delete [] fThetaRangeStringVector; *fLog << "Destructor of MFindSupercutsONOFFThetaLoop finished " << "successfully." << endl; } // -------------------------------------------------------------------------- // // Function that sets the name of the PostScript file where alpha distributions // for the different Theta bins will be stored. // It also initializes void MFindSupercutsONOFFThetaLoop::SetPostScriptFile (TPostScript* PsFile) { fPsFilename = PsFile; *fLog << "MFindSupercutsONOFF : Results (alpha distributions with excess and significances) will be stored in PostScript file " << fPsFilename -> GetName() << endl; } Bool_t MFindSupercutsONOFFThetaLoop::SetGammaEfficiency (Double_t gammaeff) { if (gammaeff <= 0.0 || gammaeff > 1) { *fLog << " ****** ERROR *********" << endl << "MFindSupercutsONOFFThetaLoop :: SetGammaEfficiency; " << endl << "The requested Gamma efficiency " << gammaeff << " can not occur" << endl << "Member data fGammaEfficiency can not be set properly." << endl << "Default value for fGammaEfficiency, which is " << fGammaEfficiency << " will remain." << endl; return kFALSE; } fGammaEfficiency = gammaeff; return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::ReadSCParamsFromAsciiFile(const char* filename, Int_t Nparams) { const Int_t NColumns = 2; Double_t tmp_vector[NColumns]; Int_t counter = 0; Double_t tmpdouble = 0.0; // File is read and data stored in data_vector ifstream infile (filename, ios::in); if ( !infile ) { cout << "Input file could not be opened" << endl; return kFALSE; // exit (1); } counter = 0; while (infile >> tmp_vector [0]) { for (Int_t i = 1; i < NColumns; i++) { infile >> tmp_vector[i]; } counter++; } infile.close(); // Length of TArray vectors is set to the length of the // vectors defined in ascii file. // If such length does not match with the specified // by variable Nparams, the function complains and // aborts. if (counter != Nparams) { *fLog << "MFindSupercutsONOFFThetaLoop::ReadSCParamsFromAsciiFile" << endl << "Length of vectors in ascii file " << filename << " is " << counter << " , which does not match with the expected length " << Nparams << endl << " fInitSCPar and fInitSCParSteps are NOT retrieved from ascii file" << endl; return kFALSE; } fInitSCPar.Set(counter); fInitSCParSteps.Set(counter); // Vectors are filled ifstream infile2 (filename, ios::in); if ( !infile2 ) { cout << "Input file could not be opened" << endl; return kFALSE; // exit (1); } for (Int_t i = 0; i < fInitSCPar.GetSize(); i++) { infile2 >> tmpdouble; fInitSCPar[i] = tmpdouble; infile2 >> tmpdouble; fInitSCParSteps[i] = tmpdouble; } // Print the vectors read *fLog << "***** INITIAL SC PARAMETERS AND STEPS TAKEN FROM ASCII FILE ********* " << endl; for (Int_t i = 0; i < fInitSCPar.GetSize(); i++) { cout << fInitSCPar[i] << setw(20) << fInitSCParSteps[i] << endl; } // // File close infile2.close(); return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetInitSCPar (TArrayD &d) { *fLog << "MFindSupercutsONOFFThetaLoop; Initial SC parameters are " << "set using function SetInitSCPar (TArrayD &d)" << endl; fInitSCPar = d; return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetInitSCParSteps (TArrayD &d) { *fLog << "MFindSupercutsONOFFThetaLoop; Initial SC parameters STEPS are " << "set using function SetInitSCParSteps (TArrayD &d)" << endl; fInitSCParSteps = d; return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaSig(Double_t alphasig) { // check that alpha is within the limits 0-90 if (alphasig <= 0 || alphasig > 90) { *fLog << "MFindSupercutsONOFFThetaLoop ::SetAlphaSig; " << "value " << alphasig << " is not within the the " << "logical limits of alpha; 0-90" << endl; return kFALSE; } fAlphaSig = alphasig; return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaBkgMin(Double_t alpha) { // check that alpha is within the limits 0-90 if (alpha <= 0 || alpha >= 90) { *fLog << "MFindSupercutsONOFFThetaLoop::SetAlphaBkgMin; " << "value " << alpha << " is not within the the " << "logical limits of alpha; 0-90" << endl; return kFALSE; } fAlphaBkgMin = alpha; return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaBkgMax(Double_t alpha) { // check that alpha is within the limits 0-90 if (alpha <= 0 || alpha > 90.001) { *fLog << "MFindSupercutsONOFFThetaLoop ::SetAlphaBkgMax; " << "value " << alpha << " is not within the the " << "logical limits of alpha; 0-90" << endl; return kFALSE; } fAlphaBkgMax = alpha; return kTRUE; } // Function that checks that the values of the member data // fAlphaSig, fAlphaBkgMin and fAlphaBkgMax make sense // (ie, fAlphaSig < fAlphaBkgMin < fAlphaBkgMax) Bool_t MFindSupercutsONOFFThetaLoop::CheckAlphaSigBkg() { if (fAlphaSig > fAlphaBkgMin) { *fLog << "MFindSupercutsONOFFThetaLoop ::CheckAlphaSigBkg(); " << endl << "fAlphaSig > fAlphaBkgMin, which should not occur..." << endl << "fAlphaSig = " << fAlphaSig << ", fAlphaBkgMin = " << fAlphaBkgMin << endl; return kFALSE; } if (fAlphaBkgMax < fAlphaBkgMin) { *fLog << "MFindSupercutsONOFFThetaLoop::CheckAlphaSigBkg(); " << endl << "fAlphaBkgMin > fAlphaBkgMax, which should not occur..." << endl << "fAlphaBkgMin = " << fAlphaBkgMin << ", fAlphaBkgMax = " << fAlphaBkgMax << endl; return kFALSE; } return kTRUE; } // Function that sets the value of fCosThetaRangeVector and // fCosThetaBinCenterVector Bool_t MFindSupercutsONOFFThetaLoop::SetCosThetaRangeVector( const TArrayD &d) { if (d.GetSize() < 2) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetCosThetaRangeVector: " << endl << "Size of d is smaller than 2... " << "fCosThetaRangeVector can not be defined properly" << endl; return kFALSE; } fCosThetaRangeVector.Set(d.GetSize()); fCosThetaRangeVector = d; // fCosThetaBinCenterVector is defined Int_t dim = fCosThetaRangeVector.GetSize() - 1; fCosThetaBinCenterVector.Set(dim); for (Int_t i = 0; i < dim; i++) { fCosThetaBinCenterVector[i] = (fCosThetaRangeVector[i+1]+fCosThetaRangeVector[i])/2; } return kTRUE; } // Function that sets the range of Theta (fThetaMin, fThetaMax) // where data will be used in the iteration // thetabin. This theta range is used when geting data from // file and filling matrices. If matrices are already built // and just read from root files, this theta range is useless. // **** IMPORTANT ********* // thetabin is an integer value that goes from 1 to // fCosThetaBinCenterVector.GetSize() // ************************ Bool_t MFindSupercutsONOFFThetaLoop::SetThetaRange (Int_t thetabin) { // Check that iteration thetabin is smaller than // components of fCosThetaBinCenterVector if (fCosThetaRangeVector.GetSize() < 2 || fCosThetaBinCenterVector.GetSize() < 1) { *fLog << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl << "Dimension of vectors fCosThetaRangeVector or/and " << "fCosThetaBinCenterVector are smaller than 2 and 1 " << "respectively. Range for theta can not be defined." << endl; return kFALSE; } if (thetabin < 1 || thetabin > fCosThetaBinCenterVector.GetSize()) { *fLog << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl << "thetabin value is " << thetabin << " Is is outside the possible values given by " << "vector fCosThetaRangeVector. " << "Range for theta can not be defined." << endl; return kFALSE; } fThetaMin = fCosThetaRangeVector[thetabin-1]; fThetaMax = fCosThetaRangeVector[thetabin]; // Now costheta values have to be converted to Radians // Function TMath::ACos(x) gives angle in Radians fThetaMin = TMath::ACos(fThetaMin); fThetaMax = TMath::ACos(fThetaMax); // For new MAGIC (July 2004), angles are given in degrees. // So I convert angles into degrees by multiplying by a the // conversion rad/degree fThetaMin = fThetaMin * TMath::RadToDeg(); fThetaMax = fThetaMax * TMath::RadToDeg(); // If fThetaMin > fThetaMax means that fCosThetaRangeVector // is defined with increasing values of CosTheta, which // means decreasing values of Theta // In such case, values of fThetaMin and fThetaMax are swapped if (fThetaMin > fThetaMax) { Double_t tmp = fThetaMin; fThetaMin = fThetaMax; fThetaMax = tmp; } // Info for user *fLog << "-------------------------------------------------"<< endl << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl << "Theta range for this iteration is set to " << fThetaMin << "-" << fThetaMax << " degrees." << endl << "-------------------------------------------------"<< endl; return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetNormFactorTrainHist() { if (fCosThetaRangeVector.GetSize() <= 1) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNormFactorTrainHist: " << endl << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin" << endl; return kFALSE; } Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1; Double_t* BinsVector = fCosThetaRangeVector.GetArray(); fNormFactorTrainHist = new TH1F ("NormFactorTrainHist", "NormFactorTrainHist", HistoNbins,BinsVector ); return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetNormFactorTestHist() { if (fCosThetaRangeVector.GetSize() <= 1) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNormFactorTestHist: " << endl << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin" << endl; return kFALSE; } Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1; Double_t* BinsVector = fCosThetaRangeVector.GetArray(); fNormFactorTestHist = new TH1F ("NormFactorTestHist", "NormFactorTestHist", HistoNbins, BinsVector); return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTrainHist() { if (fCosThetaRangeVector.GetSize() <= 1) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTrainHist: " << endl << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin" << endl; return kFALSE; } // At some point, definition of bins should be able to // contain variable bins... Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1; Double_t* BinsVector = fCosThetaRangeVector.GetArray(); fSigmaLiMaTrainHist = new TH1F ("SigmaLiMaTrainHist", "SigmaLiMaTrainHist", HistoNbins, BinsVector); return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTestHist() { if (fCosThetaRangeVector.GetSize() <= 1) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTestHist: " << endl << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin" << endl; return kFALSE; } Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1; Double_t* BinsVector = fCosThetaRangeVector.GetArray(); fSigmaLiMaTestHist = new TH1F ("SigmaLiMaTestHist", "SigmaLiMaTestHist", HistoNbins, BinsVector); return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetNexTrainHist() { if (fCosThetaRangeVector.GetSize() <= 1) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNexTrainHist: " << endl << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin" << endl; return kFALSE; } Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1; Double_t* BinsVector = fCosThetaRangeVector.GetArray(); fNexTrainHist = new TH1F ("NexTrainHist", "NexTrainHist", HistoNbins, BinsVector); return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetNexTestHist() { if (fCosThetaRangeVector.GetSize() <= 1) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNexTestHist: " << endl << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin" << endl; return kFALSE; } Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1; Double_t* BinsVector = fCosThetaRangeVector.GetArray(); fNexTestHist = new TH1F ("NexTestHist", "NexTestHist", HistoNbins, BinsVector); return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetSuccessfulThetaBinsHist() { if (fCosThetaRangeVector.GetSize() <= 1) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSuccessfulThetaBinsHist(): " << endl << "Size of fCosThetaRangeVector is <=1 ... " << "fCosThetaRangeVector must exist and have >= 2 components (defining a single bin)" << endl; return kFALSE; } Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1; Double_t* BinsVector = fCosThetaRangeVector.GetArray(); fSuccessfulThetaBinsHist = new TH1F ("SuccessfulThetaBinsHist", "SuccessfulThetaBinsHist", HistoNbins, BinsVector); return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::SetNexSigmaLiMaNormFactorNEvtsTrainTestHist() { if (fCosThetaRangeVector.GetSize() <= 1) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNexSigmaLiMaNormFactorTrainTestHist: " << endl << "Size of fCosThetaRangeVector is <=1 ... " << "fCosThetaRangeVector must exist and have, at least, " << "2 components defining a single bin" << endl; return kFALSE; } Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1; Double_t* BinsVector = fCosThetaRangeVector.GetArray(); char* x_axis_title = {"Cos(Theta)"}; fNormFactorTrainHist = new TH1F ("NormFactorTrainHist", "NormFactorTrainHist", HistoNbins, BinsVector); fNormFactorTrainHist -> SetXTitle(x_axis_title); fNormFactorTrainHist -> SetYTitle("Normalization Factor (Non/Noff)"); fNormFactorTestHist = new TH1F ("NormFactorTestHist", "NormFactorTestHist", HistoNbins, BinsVector); fNormFactorTestHist -> SetXTitle(x_axis_title); fNormFactorTestHist -> SetYTitle("Normalization Factor (Non/Noff)"); fSigmaLiMaTrainHist = new TH1F ("SigmaLiMaTrainHist", "SigmaLiMaTrainHist", HistoNbins, BinsVector); fSigmaLiMaTrainHist -> SetXTitle(x_axis_title); fSigmaLiMaTrainHist -> SetYTitle("Significance"); fSigmaLiMaTestHist = new TH1F ("SigmaLiMaTestHist", "SigmaLiMaTestHist", HistoNbins, BinsVector); fSigmaLiMaTestHist -> SetXTitle(x_axis_title); fSigmaLiMaTestHist -> SetYTitle("Significance"); fNexTrainHist = new TH1F ("NexTrainHist", "NexTrainHist", HistoNbins, BinsVector); fNexTrainHist -> SetXTitle(x_axis_title); fNexTrainHist -> SetYTitle("Excess Events"); fNexTestHist = new TH1F ("NexTestHist", "NexTestHist", HistoNbins, BinsVector); fNexTestHist -> SetXTitle(x_axis_title); fNexTestHist -> SetYTitle("Excess Events"); fNEvtsInTrainMatrixONHist = new TH1F("NEvtsInTrainMatrixONHist", "NEvtsInTrainMatrixONHist", HistoNbins, BinsVector); fNEvtsInTrainMatrixONHist -> SetXTitle(x_axis_title); fNEvtsInTrainMatrixONHist -> SetYTitle("Number of ON Events"); fNEvtsInTestMatrixONHist = new TH1F("NEvtsInTestMatrixONHist", "NEvtsInTestMatrixONHist", HistoNbins, BinsVector); fNEvtsInTestMatrixONHist -> SetXTitle(x_axis_title); fNEvtsInTestMatrixONHist -> SetYTitle("Number of ON Events"); fNEvtsInTrainMatrixOFFHist = new TH1F("NEvtsInTrainMatrixOFFHist", "NEvtsInTrainMatrixOFFHist", HistoNbins, BinsVector); fNEvtsInTrainMatrixOFFHist -> SetXTitle(x_axis_title); fNEvtsInTrainMatrixOFFHist -> SetYTitle("Number of OFF Events"); fNEvtsInTestMatrixOFFHist = new TH1F("NEvtsInTestMatrixOFFHist", "NEvtsInTestMatrixOFFHist", HistoNbins, BinsVector); fNEvtsInTestMatrixOFFHist -> SetXTitle(x_axis_title); fNEvtsInTestMatrixOFFHist -> SetYTitle("Number of OFF Events"); return kTRUE; } void MFindSupercutsONOFFThetaLoop::WriteSuccessfulThetaBinsHistToFile() { *fLog << "MFindSupercutsONOFFThetaLoop : Writing histogram " << " 'SuccessfulThetaBinsHist'" << " (if contents != 0) into root file " << fAlphaDistributionsRootFilename << endl; TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE", "Alpha Distributions for several Theta bins"); if (fSuccessfulThetaBinsHist) { if (fSuccessfulThetaBinsHist->GetEntries() > 0.5) fSuccessfulThetaBinsHist -> Write(); } rootfile.Close(); } void MFindSupercutsONOFFThetaLoop::WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile() { *fLog << "MFindSupercutsONOFFThetaLoop : Writing histograms 'NexSigmaLiMaNormFactorNEvtsTrainTest'" << " (if they have contents != 0) into root file " << fAlphaDistributionsRootFilename << endl; TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE", "Alpha Distributions for several Theta bins"); if (fNormFactorTrainHist) { if (fNormFactorTrainHist->GetEntries() > 0.5) fNormFactorTrainHist -> Write(); } if (fNormFactorTestHist) { if (fNormFactorTestHist->GetEntries() > 0.5) fNormFactorTestHist -> Write(); } if (fSigmaLiMaTrainHist) { if (fSigmaLiMaTrainHist->GetEntries() > 0.5) fSigmaLiMaTrainHist -> Write(); } if (fSigmaLiMaTestHist) { if (fSigmaLiMaTestHist->GetEntries() > 0.5) fSigmaLiMaTestHist -> Write(); } if (fNexTrainHist) { if (fNexTrainHist->GetEntries() > 0.5) fNexTrainHist -> Write(); } if (fNexTestHist) { if (fNexTestHist->GetEntries() > 0.5) fNexTestHist -> Write(); } if (fNEvtsInTrainMatrixONHist) { if (fNEvtsInTrainMatrixONHist->GetEntries() > 0.5) fNEvtsInTrainMatrixONHist -> Write(); } if (fNEvtsInTestMatrixONHist) { if (fNEvtsInTestMatrixONHist->GetEntries() > 0.5) fNEvtsInTestMatrixONHist -> Write(); } if (fNEvtsInTrainMatrixOFFHist) { if (fNEvtsInTrainMatrixOFFHist->GetEntries() > 0.5) fNEvtsInTrainMatrixOFFHist -> Write(); } if (fNEvtsInTestMatrixOFFHist) { if (fNEvtsInTestMatrixOFFHist->GetEntries() > 0.5) fNEvtsInTestMatrixOFFHist -> Write(); } rootfile.Close(); } void MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents( Double_t fontrain, Double_t fontest, Double_t fofftrain, Double_t fofftest) { // Check that fractions set by user are within the expected limits if (fontrain < 0.00 || fontrain > 1.00 || fontest < 0.00 || fontest > 1.00 || fofftrain < 0.00 || fofftrain > 1.00 || fofftest < 0.00 || fontest > 1.00) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents: " << endl << "Fractions of ON/OFF data for TRAIN/TEST sample which were specified by user " << endl << "are outside the natural limits 0.0-1.0. " << endl << "Set reasonable fraction numbers for TRAIN/TEST please. The program will be aborted." << endl; exit(1); } if((fontrain+fontest) > 1.00 || (fofftrain+fofftest) >1.00) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents: " << endl << "The sum of the fractions of TRAIN/TEST for ON or OFF data which were specified by user " << endl << "are larger than 1.0. " << endl << "Set reasonable fraction numbers for TRAIN/TEST please. The program will be aborted." << endl; exit(1); } fWhichFractionTrain = fontrain; fWhichFractionTest = fontest; fWhichFractionTrainOFF = fofftrain; fWhichFractionTestOFF = fofftest; } void MFindSupercutsONOFFThetaLoop::SetReadMatricesFromFile (Bool_t b) { *fLog << "---------------------------------------------------" << endl << "MFindSupercutsONOFFThetaLoop: " << endl; if (b) { *fLog << "Matrices are read from root files." << endl;} else { *fLog << "Matrices are produced and stored into root files." << endl; } fReadMatricesFromFile = b; } Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaDistributionsRootFilename( const TString &name) { if (fPathForFiles.IsNull()) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetAlphaDistributionsRootFilename; " << endl << "fPathForFiles is not defined yet, hence " << "name for rootfile containing alpha distributions " << "can not be defined properly" << endl; return kFALSE; } fAlphaDistributionsRootFilename = (fPathForFiles); fAlphaDistributionsRootFilename += (name); *fLog << "MFindSupercutsONOFFThetaLoop: fAlphaDistributionsRootFilename = " << fAlphaDistributionsRootFilename << endl; return kTRUE; } // Function to set Names of root files containing matrices and // optimizedparameters and also fThetaRangeStringVector vector Bool_t MFindSupercutsONOFFThetaLoop::SetALLNames() { // Names for files only possible if CosThetaRangeVector // is defined properly if (fCosThetaRangeVector.GetSize() < 2) { *fLog << err << "MFindSupercutsONOFFThetaLoop::SetALLNames; " << endl << "Size of fCosThetaRangeVector is smaller than 2... " << "fCosThetaRangeVector is not defined properly, " << "and thus, names for rootfiles containing data matrices " << "and optimized parameters for the different theta " << "ranges can not defined properly." << endl; return kFALSE; } if (fPathForFiles.IsNull()) { *fLog << err << "MFindSupercutsONOFFThetaLoop::Setnames; " << endl << "fPathForFiles is not defined yet, hence " << "names for rootfiles containing data matrices " << "and optimized parameters for the different theta " << "ranges can not defined properly." << endl; return kFALSE; } // Name of rootfiles for Supercuts parameters follow the // follwing structure // Path + "OptSCParametersONOFF" + "ThetaRange" // + int {(fThetaMin_fThetaMax)} + ".root" // Name of rootfiles for Matrices follow the follwing structure // Path + "ON(OFF)Data" + "Train(Test)" + "Matrix" + Fraction // + "ThetaRange" + int {(fThetaMin_fThetaMax)} + ".root" Int_t tmpdim = fCosThetaRangeVector.GetSize() - 1; const Int_t VectorDimension = tmpdim; fThetaRangeStringVector = new TString[VectorDimension]; fOptSCParamFilenameVector = new TString[VectorDimension]; fTrainMatrixONFilenameVector = new TString[VectorDimension]; fTestMatrixONFilenameVector = new TString[VectorDimension]; fTrainMatrixOFFFilenameVector = new TString[VectorDimension]; fTestMatrixOFFFilenameVector = new TString[VectorDimension]; // Int_t ThetaMinTimes1000 = 0; // Int_t ThetaMaxTimes1000 = 0; // For MAGIC new (July 2004) Theta is given in deg Int_t ThetaMinInt = 0; Int_t ThetaMaxInt = 0; Int_t FractionONTrainTimes100 = 0; Int_t FractionONTestTimes100 = 0; Int_t FractionOFFTrainTimes100 = 0; Int_t FractionOFFTestTimes100 = 0; char ThetaRangeString[100]; for (Int_t i = 0; i < VectorDimension; i++) { if (!SetThetaRange(i+1)) { *fLog << "MFindSupercutsONOFFThetaLoop::SetALLNames; " << endl << "Values for fThetaMin and fThetaMax could NOT " << "be computed for theta bin " << (i+1) << "; SetNames can not go on successfully." << endl; return kFALSE; } //ThetaMinTimes1000 = int(fThetaMin*1000); //ThetaMaxTimes1000 = int(fThetaMax*1000); //sprintf(ThetaRangeString, "ThetaRange%d_%dmRad", // ThetaMinTimes1000, ThetaMaxTimes1000); // For New MAGIC (July 2004) // Theta is given in deg, no need to multiply by 1000 // 0.5 is added so that the rounding to integer is correct ThetaMinInt = int(fThetaMin + 0.5); ThetaMaxInt = int(fThetaMax + 0.5); sprintf(ThetaRangeString, "ThetaRange%d_%d_Degrees", ThetaMinInt, ThetaMaxInt); fThetaRangeStringVector[i] = (ThetaRangeString); // tmp *fLog << "ThetaRangeString = " << fThetaRangeStringVector[i] << endl; // endtmp // 0.5 is added so that the rounding to integer is correct FractionONTrainTimes100 = int(fWhichFractionTrain*100 + 0.5); FractionONTestTimes100 = int(fWhichFractionTest*100 + 0.5); FractionOFFTrainTimes100 = int(fWhichFractionTrainOFF*100 + 0.5); FractionOFFTestTimes100 = int(fWhichFractionTestOFF*100 + 0.5); // File for SC parameters fOptSCParamFilenameVector[i] = (fPathForFiles); fOptSCParamFilenameVector[i] += ("OptSCParametersONOFF"); fOptSCParamFilenameVector[i] += (fThetaRangeStringVector[i]); fOptSCParamFilenameVector[i] += (".root"); // File for Data Matrices fTrainMatrixONFilenameVector[i] = (fPathForFiles); fTrainMatrixONFilenameVector[i] += ("ONDataTrainMatrixFraction"); fTrainMatrixONFilenameVector[i] += (FractionONTrainTimes100); fTrainMatrixONFilenameVector[i] += (fThetaRangeStringVector[i]); fTrainMatrixONFilenameVector[i] += (".root"); fTestMatrixONFilenameVector[i] = (fPathForFiles); fTestMatrixONFilenameVector[i] += ("ONDataTestMatrixFraction"); fTestMatrixONFilenameVector[i] += (FractionONTestTimes100); fTestMatrixONFilenameVector[i] += (fThetaRangeStringVector[i]); fTestMatrixONFilenameVector[i] += (".root"); fTrainMatrixOFFFilenameVector[i] = (fPathForFiles); fTrainMatrixOFFFilenameVector[i] += ("OFFDataTrainMatrixFraction"); fTrainMatrixOFFFilenameVector[i] += (FractionOFFTrainTimes100); fTrainMatrixOFFFilenameVector[i] += (fThetaRangeStringVector[i]); fTrainMatrixOFFFilenameVector[i] += (".root"); fTestMatrixOFFFilenameVector[i] = (fPathForFiles); fTestMatrixOFFFilenameVector[i] += ("OFFDataTestMatrixFraction"); fTestMatrixOFFFilenameVector[i] += (FractionOFFTestTimes100); fTestMatrixOFFFilenameVector[i] += (fThetaRangeStringVector[i]); fTestMatrixOFFFilenameVector[i] += (".root"); } // Info concerning names is printed *fLog << "--------------------------------------------------" << endl << "MFindSupercutsONOFFThetaLoop::Setnames; " << endl << "Names for root files containing Opt SC Param. and Data " << "matrices for the different theta bins are the " << "following ones: " << endl << "MeanCosTheta, OptSCParm, ONDataTrainMatrix, ..." << endl; for (Int_t i = 0; i < VectorDimension; i++) { *fLog << fCosThetaBinCenterVector[i] << ", " << endl << fOptSCParamFilenameVector[i] << ", "<< endl << fTrainMatrixONFilenameVector[i] << ", " << endl << fTestMatrixONFilenameVector[i] << ", " << endl << fTrainMatrixOFFFilenameVector[i] << ", " << endl << fTestMatrixOFFFilenameVector[i] << ", " << endl << endl << endl; } *fLog << "--------------------------------------------------" << endl; return kTRUE; } Bool_t MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges() { *fLog << "--------------------------------------------------------------" << endl << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges function starts " << endl; if (fDataONRootFilename.IsNull()) { *fLog << "MFindSupercutsONOFF::LoopOverThetaRanges; root file for ON data is not defined... " << "aborting..." << endl; return kFALSE; } if (fDataOFFRootFilename.IsNull()) { *fLog << "MFindSupercutsONOFF:: LoopOverThetaRanges; root file for OFF data is not defined... aborting" << endl; return kFALSE; } if (fHadronnessName.IsNull()) { *fLog << "MFindSupercutsONOFF::LoopOverThetaRanges; hadronness name for ON data is not defined... aborting" << endl; return kFALSE; } if (fHadronnessNameOFF.IsNull()) { *fLog << "MFindSupercutsONOFF::LoopOverThetaRanges; hadronness name for OFF data is not defined... aborting" << endl; return kFALSE; } if (fCosThetaRangeVector.GetSize() < 2) { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Dimension of vector fCosThetaRangeVector " << "is smaller than 2" << "Theta ranges are not well defined... aborting ..." << endl; return kFALSE; } // Check if pointers to root file names (opt SC parameters and matrices) // are empty if (!fTrainMatrixONFilenameVector || !fTestMatrixONFilenameVector || !fTrainMatrixOFFFilenameVector || !fTestMatrixOFFFilenameVector || !fOptSCParamFilenameVector) { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Pointers to root files names for Opt SC parameters and data " << "matrices are not well defined. Aborting... " << endl; return kFALSE; } if (fAlphaDistributionsRootFilename.IsNull()) { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Root filename where to store alpha distributions " << "(fAlphaDistributionsRootFilename) is not defined." << endl << "Loop over theta angles can not go on..." << endl; return kTRUE; } if (!fPsFilename) { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Object from Class TPostScript is not defined. " << "Plots will be stored in a a file named 'PsFilesTmp.ps' located " << "in current directory" << endl; // fPsFilename = new TPostScript("PsFilesTmp.ps" ,111); } // All drinks are in house... the party starts !! Int_t ThetaLoopIterations = fCosThetaRangeVector.GetSize() - 1; *fLog << "looping over the following Cos Theta Angle Ranges : "<< endl; for (Int_t i = 0; i < ThetaLoopIterations; i++) { *fLog << i+1 << " -> " << fCosThetaRangeVector[i] << "-" << fCosThetaRangeVector[i+1] << endl; } *fLog << "--------------------------------------------------------------------" << endl; // Let's create histograms where SigmaLiMa, Nex and Normfactors will be stored if (!SetNexSigmaLiMaNormFactorNEvtsTrainTestHist()) { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Histograms where to store Significances, Nex and Normalization " << "factors, NEvts for all theta bins could not be created... Aborting." << endl; return kFALSE; } // Let's create histogram where successfully optimized theta bins will be stored if (!SetSuccessfulThetaBinsHist()) { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Histogram where to store successfully optimized theta bins " << " could not be created... Aborting." << endl; return kFALSE; } Double_t NormFactor = 0.0; Double_t Nex = 0.0; Double_t SigmaLiMa = 0.0; Int_t NEvtsON = 0; Int_t NEvtsOFF = 0; Bool_t MinimizationSuccessful = kTRUE; Double_t MinimizationValue = 1.0; // Value quantifying the minimization, and stored in // histogram fSuccessfulThetaBinsHist. For the time being it is 0.0 (minimization failed). // and 1.0 (minimization succeeded). Yet in future it may take other values quantifying // other aspects... for (Int_t i = 0; i < ThetaLoopIterations; i++) { // Get the Theta range that will be used in this iteration // Remember that argument of the function is iteration number, // and it starts at 1... if (!SetThetaRange(i+1)) { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Values for fThetaMin and fThetaMax could NOT " << "be computed for theta bin " << (i+1) << "; LoopOverThetaRanges is aborting." << endl; return kFALSE; } // Object of class MFindSupercutsONOFF will be created and used MFindSupercutsONOFF FindSupercuts("MFindSupercutsONOFF", "Optimizer for the supercuts"); // Set Theta range that will be used FindSupercuts.SetThetaRange(fThetaMin, fThetaMax); // Alphamax where signal is expected is set FindSupercuts.SetAlphaSig(fAlphaSig); // Bkg alpha region is set FindSupercuts.SetAlphaBkgMin(fAlphaBkgMin); FindSupercuts.SetAlphaBkgMax(fAlphaBkgMax); // alpha bkg and signal region set in object FindSupercuts // are re-checked in order to be sure that make sense FindSupercuts.CheckAlphaSigBkg(); // Degree of polynomials used to fit ON and OFF alpha // distributions is set FindSupercuts.SetDegreeON(fDegreeON); FindSupercuts.SetDegreeOFF(fDegreeOFF); // bining for alpha plots is set FindSupercuts.SetAlphaPlotBinining(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); // Bool variable NormFactorFromAlphaBkg is set in FindSupercuts // in order to decide the way in which mormalization factors // for TRAIN and TEST samples are computed FindSupercuts.SetNormFactorFromAlphaBkg(fNormFactorFromAlphaBkg); // Define size range of events used to fill the data matrices FindSupercuts.SetSizeRange(fSizeCutLow, fSizeCutUp); FindSupercuts.SetFilters(fLeakMax, fDistMax, fDistMin); // Boolean variable that controls wether the optimization of the // parameters (MMinuitInterface::CallMinuit(..) in function FindParams(..)) // takes place or not. kTRUE will skip such optimization. // This variable is useful to test the optmized parameters (previously found // and stored in root file) on the TRAIN sample. FindSupercuts.SetSkipOptimization(fSkipOptimization); // Boolean variable that allows the user to write the initial parameters // into the root file that will be used to store the optimum cuts. // If fUseInitialSCParams = kTRUE , parameters are written. // In this way, the initial SC parameters can be applied on the data (train/test) FindSupercuts.SetUseInitialSCParams(fUseInitialSCParams); // Set boolean variable that controls wether the theta variable // is used or not in the computation of the dynamica cuts FindSupercuts.SetVariableNotUseTheta(fNotUseTheta); // Set boolean variable that controls wether cuts are dynamical // or static. FindSupercuts.SetVariableUseStaticCuts(fUseStaticCuts); FindSupercuts.SetUseFittedQuantities(fUseFittedQuantities); // fGammaEfficiency is set in object FindSupercuts FindSupercuts.SetGammaEfficiency (fGammaEfficiency); // Filename with SuperCuts parameters is set FindSupercuts.SetFilenameParam(fOptSCParamFilenameVector[i]); // Hadroness containers are set FindSupercuts.SetHadronnessName(fHadronnessName); FindSupercuts.SetHadronnessNameOFF(fHadronnessNameOFF); // Rootfilename where to store alpha distribution histograms is set FindSupercuts.SetAlphaDistributionsRootFilename(fAlphaDistributionsRootFilename); // Set the names of all trees // used to store the supercuts applied to ON/OFF Train/Test samples // In future, names might be controlled by user, but for the time // being I'll make things simple, and names are set to some defaults // taking into account the theta range used // Containers will be stored in rootfile specified by variable // fAlphaDistributionsRootFilename FindSupercuts.SetSupercutsAppliedTreeNames(); // Object of the class TPostScritp is set... // BUT THAT'S NOT WORKING PROPERLY !!!!! // FindSupercuts.SetPostScriptFile (fPsFilename); // ******************************************************** // Due to the failure of the use of object TPostScript // to make a Ps document with all plots, I decided to use the // standard way (SaveAs(filename.ps)) to store plots related // to alpha distributions // for ON and OFF and BEFORE and AFTER cuts (VERY IMPORTANT // TO COMPUTE EFFICIENCIES IN CUTS) // Psfilename is set inside function LoopOverThetaRanges() // and given to the object MFindSupercutsONOFF created // within this loop. // This will have to be removed as soon as the TPostScript // solutions works... // ******************************************************** TString PsFilenameString = (fPathForFiles); PsFilenameString += ("SC"); PsFilenameString += (fThetaRangeStringVector[i]); // Name to be completed inside object MFindSupercutsONOFF FindSupercuts.SetPsFilenameString(PsFilenameString); if (fReadMatricesFromFile) { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Reading Data matrices from root files... " << endl; /* FindSupercuts.ReadMatrix(fTrainMatrixONFilenameVector[i], fTestMatrixONFilenameVector[i]); FindSupercuts.ReadMatrixOFF(fTrainMatrixOFFFilenameVector[i], fTestMatrixOFFFilenameVector[i]); */ FindSupercuts.ReadMatrixTrain(fTrainMatrixONFilenameVector[i], fTrainMatrixOFFFilenameVector[i]); if (fTestParameters) { FindSupercuts.ReadMatrixTest(fTestMatrixONFilenameVector[i], fTestMatrixOFFFilenameVector[i]); } } else { *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Reading Data from files " << fDataONRootFilename << " and " << fDataOFFRootFilename << "(for ON and OFF data " << "respectively " << " and producing data matrices that will be stored " << "in root files. " << endl; // ON data FindSupercuts.DefineTrainTestMatrixThetaRange( fDataONRootFilename, fWhichFractionTrain, fWhichFractionTest, fThetaMin, fThetaMax, fTrainMatrixONFilenameVector[i], fTestMatrixONFilenameVector[i]); // OFF data FindSupercuts.DefineTrainTestMatrixOFFThetaRange( fDataOFFRootFilename, fWhichFractionTrainOFF, fWhichFractionTestOFF, fThetaMin, fThetaMax, fTrainMatrixOFFFilenameVector[i], fTestMatrixOFFFilenameVector[i]); } MinimizationSuccessful = kTRUE; if (fOptimizeParameters) { MinimizationSuccessful = FindSupercuts.FindParams("",fInitSCPar, fInitSCParSteps); if (MinimizationSuccessful) { // Minimization was successful // NormFactor, Nex, SigmaLiMa and NEvts Train histograms are updated *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Minimization was successful" << endl << "Updating values for NormFactor, Nex, SigmaLiMa and NEvts " << "Train histograms for " << fThetaRangeStringVector[i] << endl; NormFactor = FindSupercuts.GetNormFactorTrain(); SigmaLiMa = FindSupercuts.GetSigmaLiMaTrain(); Nex = FindSupercuts.GetNexTrain(); NEvtsON = FindSupercuts.GetMatrixTrain()->GetM().GetNrows(); NEvtsOFF = FindSupercuts.GetMatrixTrainOFF()->GetM().GetNrows(); // MinimizationValue is updated... MinimizationValue = 1.0; } else { // Minimization was NOT successful // NormFactor, Nex, SigmaLiMa and NEvts Train histograms are updated *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Minimization was NOT successful" << endl << "Updating values for NormFactor, and NEvts " << "Train histograms for " << fThetaRangeStringVector[i] << endl << "Nex and SigmaLiMa are set to ZERO" << endl; NormFactor = FindSupercuts.GetNormFactorTrain(); NEvtsON = FindSupercuts.GetMatrixTrain()->GetM().GetNrows(); NEvtsOFF = FindSupercuts.GetMatrixTrainOFF()->GetM().GetNrows(); SigmaLiMa = 0.0; Nex = 0.0; // MinimizationValue is updated... MinimizationValue = 0.0; } cout << "---- Train Histograms will be updated with the following numbers: ------ " << endl << "Mean CosTheta for this Thetabin = " << fCosThetaBinCenterVector[i] << endl << "Minimization status (1.0 is successful, 0.0 is UNsuccessful) = " << MinimizationValue << endl << "Number of ON Events in Train sample used = " << NEvtsON << endl << "Number of OFF Events in Train sample used = " << NEvtsOFF << endl << "Normalization factor (Non/Noff) for Train sample = " << NormFactor << endl << "Excess events and SigmaLiMa: " << Nex << "; " << SigmaLiMa << endl << endl; fNormFactorTrainHist -> Fill(fCosThetaBinCenterVector[i], NormFactor); fNexTrainHist -> Fill(fCosThetaBinCenterVector[i], Nex); fSigmaLiMaTrainHist -> Fill(fCosThetaBinCenterVector[i], SigmaLiMa); fNEvtsInTrainMatrixONHist -> Fill(fCosThetaBinCenterVector[i], NEvtsON); fNEvtsInTrainMatrixOFFHist -> Fill(fCosThetaBinCenterVector[i], NEvtsOFF); fSuccessfulThetaBinsHist -> Fill(fCosThetaBinCenterVector[i], MinimizationValue); cout << "Histograms (TRAIN) updated successfully " << endl; } if (fTestParameters && MinimizationSuccessful) { // Cuts are applied to test sample FindSupercuts.TestParamsOnTestSample(); // NormFactor, Nex and SigmaLiMa Test histograms are updated *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Updating values for NormFactor, Nex, SigmaLiMa and NEvts " << "Test histograms for " << fThetaRangeStringVector[i] << endl; NormFactor = FindSupercuts.GetNormFactorTest(); SigmaLiMa = FindSupercuts.GetSigmaLiMaTest(); Nex = FindSupercuts.GetNexTest(); NEvtsON = FindSupercuts.GetMatrixTest()->GetM().GetNrows(); NEvtsOFF = FindSupercuts.GetMatrixTestOFF()->GetM().GetNrows(); cout << "---- Test Histograms will be updated with the following numbers: ------ " << endl << "Mean CosTheta for this Thetabin = " << fCosThetaBinCenterVector[i] << endl << "Number of ON Events in Test sample used = " << NEvtsON << endl << "Number of OFF Events in Test sample used = " << NEvtsOFF << endl << "Normalization factor (Non/Noff) for Test sample = " << NormFactor << endl << "Excess events and SigmaLiMa: " << Nex << "; " << SigmaLiMa << endl << endl; fNormFactorTestHist -> Fill(fCosThetaBinCenterVector[i], NormFactor); fNexTestHist -> Fill(fCosThetaBinCenterVector[i], Nex); fSigmaLiMaTestHist -> Fill(fCosThetaBinCenterVector[i], SigmaLiMa); fNEvtsInTestMatrixONHist -> Fill(fCosThetaBinCenterVector[i], NEvtsON); fNEvtsInTestMatrixOFFHist -> Fill(fCosThetaBinCenterVector[i], NEvtsOFF); cout << "Histograms (TEST) updated successfully " << endl; } } // PsFile is closed //*fLog << "Closing PsFile "; //*fLog << fPsFilename -> GetName() << endl; //fPsFilename -> Close(); if (fOptimizeParameters || fTestParameters) { // Histograms cotanining SigmaLiMa, Nex and NormFactors and NEvts are written // in the same root file where alpha distributions where stored // i.e. fAlphaDistributionsRootFilename *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Executing WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile()" << endl; WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile(); if (fOptimizeParameters) { // Histogram containing minimization status for all theta bins defined // in vector fCosThetaVector is stored in file fAlphaDistributionsRootFilename *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl << "Executing WriteSuccessfulThetaBinsHistToFile()" << endl; WriteSuccessfulThetaBinsHistToFile(); } } return kTRUE; } // Function that loops over the alpha distributions (ON-OFF) // stored in root file defined by fAlphaDistributionsRootFilename // and computes the significance and Nex (using MHFindSignificanceONOFF::FindSigma) // for several cuts in alpha (0-fAlphaSig; in bins defined for alpha distributions // by user). // It initializes the histograms (memeber data), fills them and store them // in root file defined by fAlphaDistributionsRootFilename. A single histogram // for each theta bin. Bool_t MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig() { *fLog << endl << endl << "****************************************************************************" << endl << "Computing number of excess events and significance vs " << endl << " the cut in alpha parameter (fAlphaSig) using function " << endl << " MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig() " << endl << "****************************************************************************" << endl << endl; // Check that all elements needed by the function are available if (fAlphaDistributionsRootFilename.IsNull()) { *fLog << err << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig; " << endl << "fAlphaDistributionsRootFilename is not defined yet, hence " << "histograms containing alpha distributions for the different theta " << "ranges can not be retrieved. " << endl; return kFALSE; } if (fCosThetaRangeVector.GetSize() < 2) { *fLog << err << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig; " << endl << "Size of fCosThetaRangeVector is smaller than 2... " << "fCosThetaRangeVector is not defined properly, " << "and thus, number of theta loops and names " << " of the histograms containing alpha distributions " << "for the different theta " << "ranges can not defined properly." << endl; return kFALSE; } // It seems that i have all I need for the time being. Let's start the party... Int_t tmpint = 0; // variable that will be used to store temporally integers Double_t tmpdouble = 0.0; // variable that will be used to store temporally doubles tmpint = fCosThetaRangeVector.GetSize() - 1; const Int_t NThetaBins = tmpint; // Definition of the bins for Nex and Significance vs Alpha histograms // The bins in Alpha will start with ZERO and will end // with 5 bins further after fAlphaSig. Double_t AlphaBinWidth = double((fAlphaBinUp-fAlphaBinLow)/fNAlphaBins); Double_t FirstBinLow = 0.0; tmpint = int((fAlphaSig/AlphaBinWidth) + 5.0); while ((FirstBinLow + AlphaBinWidth*tmpint) > fAlphaBkgMin) { tmpint--; } *fLog << "Information about the binning of the histograms Nex and Significance vs alphasig " << endl << "FirstBinLow, AlphaBinWidth, NBins : " << FirstBinLow << ", " << AlphaBinWidth << ", " << tmpint-1 << endl; const Int_t NexSignVSAlphaXBinsVectorDimension = tmpint; Double_t* XBinsVector = new Double_t [NexSignVSAlphaXBinsVectorDimension]; for (Int_t i = 0; i < NexSignVSAlphaXBinsVectorDimension; i++) { XBinsVector[i] = FirstBinLow + i*AlphaBinWidth; // cout << XBinsVector[i] << endl; } // Axis labels for Nex and Significance vs Alpha histograms TString x_axis_label = ("Cut in alpha [\\circ]"); TString y_axis_label_nex = ("Excess Events"); TString y_axis_label_significance = ("Significance (LiMa formula 17)"); // Arguments for function MHFindSignificanceONOFF::FindSigmaONOFF() // are defined in here. The same arguments (obviously) will be used // in all theta bins and TRAIN and TEST samples. // variable alphasig will be defined and varied withing the theta loop Double_t alphasig = 0.0; // Minimum value of alpha for bkg region in ON data const Double_t alphabkgmin = fAlphaBkgMin; // Maximum value of alpha for bkg region in ON data const Double_t alphabkgmax = fAlphaBkgMax; const Bool_t drawpoly = kFALSE; const Bool_t fitgauss = kFALSE; const Bool_t print = kFALSE; const Bool_t saveplots = kFALSE; // No plots are saved const Bool_t RebinHistograms = kTRUE; const Bool_t ReduceDegree = kFALSE; Double_t Nex = 0.0; Double_t Significance = 0.0; // Names of histograms containing Nex and Significance vs alphasig TString NexVSAlphaName = ("NexVSAlpha"); // Names to be completed in theta loop TString SignificanceVSAlphaName = ("SignificanceVSAlpha");// Names to be completed in theta loop // Names of histograms to be retrieved from root file // fAlphaDistributionsRootFilename TString NormFactorTrainHistoName = ("NormFactorTrainHist"); TString NormFactorTestHistoName = ("NormFactorTestHist"); TString SuccessfulThetaBinsHistoName = ("SuccessfulThetaBinsHist"); // Histogram containing information about the success of the optimization // for the several theta bins. The content of those bins where optimization was not // successful, are set to ZERO. TH1F* SuccessfulThetaBinsHisto; // Histograms to store the normalization factors TH1F* NormFactorTrainHisto; TH1F* NormFactorTestHisto; // Histograms that will be used to store the alpha distributions // retrieved from fAlphaDistributionsRootFilename, and that will be use as arguments in // function MHFindSignificanceONOFF::FindSigma // The same histograms (pointer to histograms actually) will be used for TRAIN and TEST // and for the different ThetaNBinsUsed. TH1F* AlphaONAfterCuts; TH1F* AlphaOFFAfterCuts; // Vector containing the histogram index of the successfully // optimized theta bins. // It is filled with the contents of the histogram // fSuccessfulThetaBinsHist retrieved from the root // file fAlphaDistributionsRootFilename // Remember that histogram index START AT 1 !!! TArrayI SuccessfulThetaBinsVector; // Open root file for retrieving information *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Opening root file " << fAlphaDistributionsRootFilename << endl; TFile rootfile (fAlphaDistributionsRootFilename, "READ"); // Retrieving SuccessfulThetaBinsHisto from root file. *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Retrieving SuccessfulThetaBinsHisto from root file and filling vector " << " SuccessfulThetaBinsVector" << endl; cout << "Getting histogram with name " << SuccessfulThetaBinsHistoName << endl; SuccessfulThetaBinsHisto = (TH1F*) rootfile.Get(SuccessfulThetaBinsHistoName); // Check that bins in this histo correspond with the // CosTheta intervals defined by vector fCosThetaRangeVector if (SuccessfulThetaBinsHisto -> GetNbinsX() != NThetaBins) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Number of theta bins defined by vector fCosThetaRangeVector (" << NThetaBins << ") does not match with the bins in histogram " << SuccessfulThetaBinsHistoName << "(" << SuccessfulThetaBinsHisto -> GetNbinsX() << ")" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } // Filling vector SuccessfulThetaBinsVector tmpint = 0; for (Int_t i = 0; i < NThetaBins; i++) { if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5) { tmpint++; } } SuccessfulThetaBinsVector.Set(tmpint); tmpint = 0; for (Int_t i = 0; i < NThetaBins; i++) { if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5) { if(tmpint < SuccessfulThetaBinsVector.GetSize()) { SuccessfulThetaBinsVector[tmpint] = i+1; tmpint++; } else { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Problem when filling vector 'SuccessfulThetaBinsVector'" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } } } // variable defining the theta bins where optimization was successful, and // hence, alpha distributions stored in root file are meaningful (up to some extent) tmpint = SuccessfulThetaBinsVector.GetSize(); const Int_t SuccessfulNThetaBins = tmpint; // TEMP // cout << "const Int_t SuccessfulNThetaBins = " << SuccessfulNThetaBins << endl; // ENDTEMP // Vectors of pointers to histograms are defined to store the histograms // containing Nex and significance for TRAIN and TEST sample for all the // successful theta bins // Histograms will be initialized in the loop over theta bins TH1F* NexVSAlphaSigTrainVector[SuccessfulNThetaBins]; TH1F* NexVSAlphaSigTestVector[SuccessfulNThetaBins]; TH1F* SignificanceVSAlphaSigTrainVector[SuccessfulNThetaBins]; TH1F* SignificanceVSAlphaSigTestVector[SuccessfulNThetaBins]; // *************************************************** // HISTOGRAMS FOR TRAIN SAMPLE ARE COMPUTED // *************************************************** if (fOptimizeParameters) { // Computation for TRAIN sample *fLog << endl << endl << "****************************************************************************" << endl << "Computing Nex and significance vs cut in alpha (alphasig) for TRAIN sample" << endl << "****************************************************************************" << endl << endl << endl; // Retrieving NormFactorTrainHisto from root file. cout << "Getting histogram with name " << NormFactorTrainHistoName << endl; NormFactorTrainHisto = (TH1F*) rootfile.Get(NormFactorTrainHistoName); // Check that bins in this histo correspond with the // CosTheta intervals defined by vector fCosThetaRangeVector ////////////////// START CHECK ////////////////////////////////////// if (NormFactorTrainHisto -> GetNbinsX() != NThetaBins) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Number of theta bins defined by vector fCosThetaRangeVector (" << NThetaBins << ") does not match with the bins in histogram " << NormFactorTrainHistoName << "(" << NormFactorTrainHisto -> GetNbinsX() << ")" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } // histo test /* cout << "NORMFACTORTRAINHIST Test: BinCenter and Value" << endl; for (Int_t k = 0; k < NThetaBins; k++) { cout << NormFactorTrainHisto -> GetBinCenter(k+1) << "; " << NormFactorTrainHisto -> GetBinContent(k+1)<< endl; } */ for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { tmpdouble = NormFactorTrainHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]); if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble || fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Bins defined by vector fCosThetaRangeVector " << "do not match with the bins of histogram " << NormFactorTrainHistoName << endl << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl << "NormFactorHist Bin: " << tmpdouble << endl << "Function execution will be ABORTED." << endl; return kFALSE; } } ////////////////// END CHECK ////////////////////////////////////// // Check that the theta range string (needed to compute names for alpha histograms) // is well defined for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull()) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Component " << SuccessfulThetaBinsVector[i]-1 << " of fThetaRangeStringVector is EMPTY" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } } // ************************************************************* // ************************************************************* // Normfactors and successfultheta bins are ok. // I can now compute the Nex and Sigma for train sample (ON-OFF) // ************************************************************* // ************************************************************* // Loop over the several successful theta bins doing the following // 0) Initializing histograms Nex and Significance (with names!!) // 1) Retrieving alpha histograms (ON-OFF), and normalization factor // 2) Computing Nex and Signficance vs alpha and filling histograms *fLog << endl << "Starting loop over the successfully optimized theta ranges " << endl << endl; Double_t NormalizationFactor = 0.0; for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { *fLog << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] << endl << endl; // 0) Initializing histograms Nex and Significance (with names!!) TString NexVSAlphaNameTmp = (NexVSAlphaName); TString SignificanceVSAlphaNameTmp = (SignificanceVSAlphaName); NexVSAlphaNameTmp += "Train"; NexVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; SignificanceVSAlphaNameTmp += "Train"; SignificanceVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; // TEMP //cout << NexVSAlphaNameTmp << endl; //cout << SignificanceVSAlphaNameTmp << endl; // ENDTEMP NexVSAlphaSigTrainVector[i] = new TH1F (NexVSAlphaNameTmp.Data(), NexVSAlphaNameTmp.Data(), NexSignVSAlphaXBinsVectorDimension - 1, XBinsVector); // It is important to unlink the the created histograms from // the current directory (defined now by "rootfile"), // so that we can do whatever we want to // with them regardless of such directory is opened or closed // The reference of the histogram is removed from the current // directory by means of the member function "SetDirectory(0)" NexVSAlphaSigTrainVector[i] -> SetDirectory(0); SignificanceVSAlphaSigTrainVector[i] = new TH1F (SignificanceVSAlphaNameTmp.Data(), SignificanceVSAlphaNameTmp.Data(), NexSignVSAlphaXBinsVectorDimension - 1, XBinsVector); SignificanceVSAlphaSigTrainVector[i] -> SetDirectory(0); // axis are defined NexVSAlphaSigTrainVector[i] -> SetXTitle(x_axis_label.Data()); NexVSAlphaSigTrainVector[i] -> SetYTitle(y_axis_label_nex.Data()); SignificanceVSAlphaSigTrainVector[i] ->SetXTitle(x_axis_label.Data()); SignificanceVSAlphaSigTrainVector[i] ->SetYTitle(y_axis_label_significance.Data()); // 1) Retrieving alpha histograms (ON-OFF) and normalization factor NormalizationFactor = NormFactorTrainHisto -> GetBinContent(SuccessfulThetaBinsVector[i]); TString TmpAlphaHistoName = ("TrainAlphaAfterCuts"); TmpAlphaHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; *fLog << "Getting Histogram with name " << endl << TmpAlphaHistoName << endl; AlphaONAfterCuts = (TH1F*) rootfile.Get(TmpAlphaHistoName); TString TmpAlphaOFFHistoName = ("TrainAlphaOFFAfterCuts"); TmpAlphaOFFHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; *fLog << "Getting Histogram with name " << endl << TmpAlphaOFFHistoName << endl; AlphaOFFAfterCuts = (TH1F*) rootfile.Get(TmpAlphaOFFHistoName); // 2) Computing Nex and Signficance vs alpha and filling histograms // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF *fLog << endl << endl << "Starting loop to compute Nex and Significance vs alphasig " << endl << "for TRAIN sample and " << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] << " :" << endl; for (Int_t j = 1; j < NexSignVSAlphaXBinsVectorDimension; j++) { alphasig = XBinsVector[j]; MHFindSignificanceONOFF findsig; findsig.SetRebin(RebinHistograms); findsig.SetReduceDegree(ReduceDegree); findsig.SetUseFittedQuantities(fUseFittedQuantities); // Dummy name for psfile TString psfilename = ("DummyPs"); findsig.FindSigmaONOFF(AlphaONAfterCuts, AlphaOFFAfterCuts, NormalizationFactor, alphabkgmin, alphabkgmax, fDegreeON, fDegreeOFF, alphasig, drawpoly, fitgauss, print, saveplots, psfilename); Significance = double(findsig.GetSignificance()); // so far there is not a variable that contains // NexONOFF or NexONOFFFitted (in class MHFindSignificanceONOFF) // according to the value of the variable fUseFittedQuantities. // As soon as I have some time I will implement it (with the // corresponging GET member function) and will change // MFindSignificanceONOFF accordingly to make use of it. // For the time being, I will survive with an "if statement". if(fUseFittedQuantities) { Nex = double(findsig.GetNexONOFFFitted()); } else { Nex = double(findsig.GetNexONOFF()); } *fLog << endl << "Cut in alpha, Nex and significance : " << alphasig << ", " << Nex << ", " << Significance << endl << endl; // updating histograms tmpdouble = alphasig - AlphaBinWidth/2.0; NexVSAlphaSigTrainVector[i] -> Fill(tmpdouble, Nex); SignificanceVSAlphaSigTrainVector[i] -> Fill(tmpdouble, Significance); } // Dynamic memory allocated in this loop is released *fLog << "Memory allocated by temporal alpha histograms " << "will be released" << endl; delete AlphaONAfterCuts; delete AlphaOFFAfterCuts; } } // *************************************************** // HISTOGRAMS FOR TEST SAMPLE ARE COMPUTED // *************************************************** if(fTestParameters) {// computation for TEST sample *fLog << endl << endl << "**************************************************************************" << endl << "Computing Nex and significance vs cut in alpha (alphasig) for TEST sample" << endl << "**************************************************************************" << endl << endl << endl; // Retrieving NormFactorTestHisto from root file. *fLog << "Getting histogram with name " << NormFactorTestHistoName << endl; NormFactorTestHisto = (TH1F*) rootfile.Get(NormFactorTestHistoName); // Check that bins in this histo correspond with the // CosTheta intervals defined by vector fCosThetaRangeVector ////////////////// START CHECK ////////////////////////////////////// if (NormFactorTestHisto -> GetNbinsX() != NThetaBins) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Number of theta bins defined by vector fCosThetaRangeVector (" << NThetaBins << ") does not match with the bins in histogram " << NormFactorTestHistoName << "(" << NormFactorTestHisto -> GetNbinsX() << ")" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } // histo test /* cout << "NORMFACTORTRAINHIST Test: BinCenter and Value" << endl; for (Int_t k = 0; k < NThetaBins; k++) { cout << NormFactorTestHisto -> GetBinCenter(k+1) << "; " << NormFactorTestHisto -> GetBinContent(k+1)<< endl; } */ for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { tmpdouble = NormFactorTestHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]); if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble || fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Bins defined by vector fCosThetaRangeVector " << "do not match with the bins of histogram " << NormFactorTestHistoName << endl << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl << "NormFactorHist Bin: " << tmpdouble << endl << "Function execution will be ABORTED." << endl; return kFALSE; } } ////////////////// END CHECK ////////////////////////////////////// // Check that the theta range string (needed to compute names for alpha histograms) // is well defined for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull()) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl << "Component " << SuccessfulThetaBinsVector[i]-1 << " of fThetaRangeStringVector is EMPTY" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } } // ************************************************************* // ************************************************************* // Normfactors and successfultheta bins are ok. // I can now compute the Nex and Sigma for test sample (ON-OFF) // ************************************************************* // ************************************************************* // Loop over the several successful theta bins doing the following // 0) Initializing histograms Nex and Significance (with names!!) // 1) Retrieving alpha histograms (ON-OFF), and normalization factor // 2) Computing Nex and Signficance vs alpha and filling histograms *fLog << endl << "Starting loop over the successfully optimized theta ranges " << endl << endl; Double_t NormalizationFactor = 0.0; for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { *fLog << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] << endl << endl; // 0) Initializing histograms Nex and Significance (with names!!) TString NexVSAlphaNameTmp = (NexVSAlphaName); TString SignificanceVSAlphaNameTmp = (SignificanceVSAlphaName); NexVSAlphaNameTmp += "Test"; NexVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; SignificanceVSAlphaNameTmp += "Test"; SignificanceVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; // TEMP // cout << NexVSAlphaNameTmp << endl; // cout << SignificanceVSAlphaNameTmp << endl; // ENDTEMP NexVSAlphaSigTestVector[i] = new TH1F (NexVSAlphaNameTmp.Data(), NexVSAlphaNameTmp.Data(), NexSignVSAlphaXBinsVectorDimension - 1, XBinsVector); // It is important to unlink the the created histograms from // the current directory (defined now by "rootfile"), // so that we can do whatever we want to // with them regardless of such directory is opened or closed // The reference of the histogram is removed from the current // directory by means of the member function "SetDirectory(0)" NexVSAlphaSigTestVector[i] -> SetDirectory(0); SignificanceVSAlphaSigTestVector[i] = new TH1F (SignificanceVSAlphaNameTmp.Data(), SignificanceVSAlphaNameTmp.Data(), NexSignVSAlphaXBinsVectorDimension - 1, XBinsVector); SignificanceVSAlphaSigTestVector[i] -> SetDirectory(0); // axis are defined NexVSAlphaSigTestVector[i] -> SetXTitle(x_axis_label.Data()); NexVSAlphaSigTestVector[i] -> SetYTitle(y_axis_label_nex.Data()); SignificanceVSAlphaSigTestVector[i] ->SetXTitle(x_axis_label.Data()); SignificanceVSAlphaSigTestVector[i] ->SetYTitle(y_axis_label_significance.Data()); // 1) Retrieving alpha histograms (ON-OFF) and normalization factor NormalizationFactor = NormFactorTestHisto -> GetBinContent(SuccessfulThetaBinsVector[i]); TString TmpAlphaHistoName = ("TestAlphaAfterCuts"); TmpAlphaHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; *fLog << "Getting Histogram with name " << endl << TmpAlphaHistoName << endl; AlphaONAfterCuts = (TH1F*) rootfile.Get(TmpAlphaHistoName); TString TmpAlphaOFFHistoName = ("TestAlphaOFFAfterCuts"); TmpAlphaOFFHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; *fLog << "Getting Histogram with name " << endl << TmpAlphaOFFHistoName << endl; AlphaOFFAfterCuts = (TH1F*) rootfile.Get(TmpAlphaOFFHistoName); // 2) Computing Nex and Signficance vs alpha and filling histograms // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF *fLog << endl << endl << "Starting loop to compute Nex and Significance vs alphasig " << endl << "for TRAIN sample and " << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] << " :" << endl; for (Int_t j = 1; j < NexSignVSAlphaXBinsVectorDimension; j++) { alphasig = XBinsVector[j]; MHFindSignificanceONOFF findsig; findsig.SetRebin(RebinHistograms); findsig.SetReduceDegree(ReduceDegree); findsig.SetUseFittedQuantities(fUseFittedQuantities); // Dummy name for psfile TString psfilename = ("DummyPs"); findsig.FindSigmaONOFF(AlphaONAfterCuts, AlphaOFFAfterCuts, NormalizationFactor, alphabkgmin, alphabkgmax, fDegreeON, fDegreeOFF, alphasig, drawpoly, fitgauss, print, saveplots, psfilename); Significance = double(findsig.GetSignificance()); // so far there is not a variable that contains // NexONOFF or NexONOFFFitted (in class MHFindSignificanceONOFF) // according to the value of the variable fUseFittedQuantities. // As soon as I have some time I will implement it (with the // corresponging GET member function) and will change // MFindSignificanceONOFF accordingly to make use of it. // For the time being, I will survive with an "if statement". if(fUseFittedQuantities) { Nex = double(findsig.GetNexONOFFFitted()); } else { Nex = double(findsig.GetNexONOFF()); } *fLog << endl << "Cut in alpha, Nex and significance : " << alphasig << ", " << Nex << ", " << Significance << endl << endl; // updating histograms tmpdouble = alphasig - AlphaBinWidth/2.0; NexVSAlphaSigTestVector[i] -> Fill(tmpdouble, Nex); SignificanceVSAlphaSigTestVector[i] -> Fill(tmpdouble, Significance); } // Dynamic memory allocated in this loop is released *fLog << "Memory allocated by temporal alpha histograms " << "will be released" << endl; delete AlphaONAfterCuts; delete AlphaOFFAfterCuts; } } rootfile.Close(); // ************************************************************* // Histograms Nex and Significance vs alpha are written into // root file fAlphaDistributionsRootFilename *fLog <<"Histograms containing Nex and significance vs alphasig " << " are saved into root file " << endl << fAlphaDistributionsRootFilename << endl; TFile rootfiletowrite (fAlphaDistributionsRootFilename, "UPDATE", "Alpha Distributions for several Theta bins"); for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { if (fOptimizeParameters) {// Histograms for train sample are written to root file if(NexVSAlphaSigTrainVector[i]) { if(NexVSAlphaSigTrainVector[i]-> GetEntries() > 0.5) NexVSAlphaSigTrainVector[i]-> Write(); } if(SignificanceVSAlphaSigTrainVector[i]) { if(SignificanceVSAlphaSigTrainVector[i]-> GetEntries() > 0.5) SignificanceVSAlphaSigTrainVector[i]-> Write(); } } if (fTestParameters) {// Histograms for test sample are written to root file if(NexVSAlphaSigTestVector[i]) { if(NexVSAlphaSigTestVector[i]-> GetEntries() > 0.5) NexVSAlphaSigTestVector[i]-> Write(); } if(SignificanceVSAlphaSigTestVector[i]) { if(SignificanceVSAlphaSigTestVector[i]-> GetEntries() > 0.5) SignificanceVSAlphaSigTestVector[i]-> Write(); } } } rootfiletowrite.Close(); // Dynamic memory allocated is released delete XBinsVector; // function finished successfully... hopefully... return kTRUE; } // Function that gets the histograms with the alpha distributions // (for all the theta bins specified by fCosThetaRangeVector) // stored in fAlphaDistributionsRootFilename, and combine them // (correcting OFF histograms with the normalization factors stored // in NormFactorTrain or NormFactorTest) to get one single // Alpha distribution for ON and another one for OFF. // Then these histograms are given as arguments to // the function MHFindSignificanceONOFF::FindSigmaONOFF, // (Object of this class is created) to compute the // Overall Excess events and significance, that will be // stored in variables fOverallNexTrain/Test and fOverallSigmaLiMaTrain/Test Bool_t MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance(Bool_t CombineTrainData, Bool_t CombineTestData) { *fLog << endl << endl << "****************************************************************************" << endl << "Combining number of excess events and significance for the " << endl << "selected bins in theta angle using the function " << endl << " MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance() " << endl << "****************************************************************************" << endl << endl; // First of all let's check that we have the proper the "drinks" for the party if (fAlphaDistributionsRootFilename.IsNull()) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Root file containing all alpha distributions " <<" (fAlphaDistributionsRootFilename) is not well defined." << endl << "Execution of the function is ABORTED" << endl; return kFALSE; } if (fCosThetaRangeVector.GetSize() < 2) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Dimension of vector fCosThetaRangeVector is smaller than 2." << "So, vector fCosThetaRangeVector is NOT properly defined" << endl << "Function execution will be ABORTED" << endl; return kFALSE; } TArrayI SuccessfulThetaBinsVector; // Vector containing the histogram index of the successfully // optimized theta bins. // It is filled with the contents of the histogram // fSuccessfulThetaBinsHist retrieved from the root // file fAlphaDistributionsRootFilename // Remember that histogram index START AT 1 !!! // Histograms that will contain the alpha distributions // for the combined theta bins TH1F* CombinedAlphaTrainON = NULL; TString CombinedAlphaTrainONName = ("TrainSampleAlphaONCombinedThetaBins"); TH1F* CombinedAlphaTrainOFF = NULL; TString CombinedAlphaTrainOFFName = ("TrainSampleAlphaOFFCombinedThetaBins"); TH1F* CombinedAlphaTestON = NULL; TString CombinedAlphaTestONName = ("TestSampleAlphaONCombinedThetaBins"); TH1F* CombinedAlphaTestOFF = NULL; TString CombinedAlphaTestOFFName = ("TestSampleAlphaONCombinedThetaBins"); Int_t NAlphaBins = 0; // Amount of bins of alpha histograms is expected to // be equal for all alpha histograms // Vectors storing the computed error for each of the NAlphaBins // Only 2 vectors are needed for OFF Train and OFF Test sample, whose // alpha distributions are corrected with normalization factors. // The ON sample is just added (no normalization factors are applied) // and therefore, the error of the bin is the assumed by default // (i.e Sqrt(BinContent) TArrayD ErrorInCombinedAlphaTrainOFF; TArrayD ErrorInCombinedAlphaTestOFF; TArrayD CombinedAlphaTrainOFFWithoutNormalization; TArrayD CombinedAlphaTestOFFWithoutNormalization; const Double_t SmallQuantity = 0.01; // This quantity will be used as the // maximum difference between quantities (double) that are expected to be equal, // as the center of the bins in the several alpha distributions TH1F* NormFactorTrainHisto; TString NormFactorTrainHistoName = ("NormFactorTrainHist"); TH1F* NormFactorTestHisto; TString NormFactorTestHistoName = ("NormFactorTestHist"); TString SuccessfulThetaBinsHistoName = ("SuccessfulThetaBinsHist"); TH1F* SuccessfulThetaBinsHisto; TH1F* TmpTH1FHisto; // Retrieve number of theta bins that have to be combined // from TArrayD fCosThetaRangeVector Int_t tmpdim = fCosThetaRangeVector.GetSize() - 1; const Int_t NThetaBins = tmpdim; // Variables to store the mean normalization factor of the // combined OFF sample. // MeanNormFactor = Sum (NormFactor_i*Noff_i)/Sum (Noff_i) Double_t MeanNormFactorTrain = 0.0; Double_t MeanNormFactorTest = 0.0; Double_t TotalNumberOFFTrain = 0.0; Double_t TotalNumberOFFTest = 0.0; Double_t EventCounter = 0.0;// variable used to count events used to fill histograms // it is double (and not Int) because the contents of an histogram bin might not // be an integer value... due to normalization factors !! Int_t tmpint = 0; Double_t tmpdouble = 0.0; // tmp variable double Double_t tmpdouble2 = 0.0; // tmp variable double Double_t tmperror = 0.0; // quantity used to store // temporaly the bin error in OFF histogram Double_t tmpnormfactor = 0.0; // variable used to store normalization factor // of the theta bin i (used for TRAIN and TEST sample) // Open root file *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Opening root file " << fAlphaDistributionsRootFilename << endl; TFile rootfile (fAlphaDistributionsRootFilename, "READ"); // Retrieving SuccessfulThetaBinsHisto from root file. *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Retrieving SuccessfulThetaBinsHisto from root file and filling vector " << " SuccessfulThetaBinsVector" << endl; cout << "Getting histogram " << SuccessfulThetaBinsHistoName << endl; SuccessfulThetaBinsHisto = (TH1F*) rootfile.Get(SuccessfulThetaBinsHistoName); // Check that bins in this histo correspond with the // CosTheta intervals defined by vector fCosThetaRangeVector if (SuccessfulThetaBinsHisto -> GetNbinsX() != NThetaBins) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Number of theta bins defined by vector fCosThetaRangeVector (" << NThetaBins << ") does not match with the bins in histogram " << SuccessfulThetaBinsHistoName << "(" << SuccessfulThetaBinsHisto -> GetNbinsX() << ")" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } // Filling vector SuccessfulThetaBinsVector tmpint = 0; for (Int_t i = 0; i < NThetaBins; i++) { if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5) { tmpint++; } } SuccessfulThetaBinsVector.Set(tmpint); tmpint = 0; for (Int_t i = 0; i < NThetaBins; i++) { if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5) { if(tmpint < SuccessfulThetaBinsVector.GetSize()) { SuccessfulThetaBinsVector[tmpint] = i+1; tmpint++; } else { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Problem when filling vector 'SuccessfulThetaBinsVector'" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } } } // ********* HISTOGRAMS FOR TRAIN SAMPLE WILL BE FILLED *********/// if (CombineTrainData) { *fLog << endl << "****************************************************************************" << endl << "Combining data for the TRAIN sample " << endl << "****************************************************************************" << endl << endl; // Retrieving NormFactorTrainHisto from root file. NormFactorTrainHisto = (TH1F*) rootfile.Get(NormFactorTrainHistoName); // NormFactorTrainHisto-> DrawCopy(); // Check that bins in this histo correspond with the // CosTheta intervals defined by vector fCosThetaRangeVector if (NormFactorTrainHisto -> GetNbinsX() != NThetaBins) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Number of theta bins defined by vector fCosThetaRangeVector (" << NThetaBins << ") does not match with the bins in histogram " << NormFactorTrainHistoName << "(" << NormFactorTrainHisto -> GetNbinsX() << ")" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } // histo test /* cout << "NORMFACTORTRAINHIST TEST: BinCenter and Value" << endl; for (Int_t k = 0; k < NThetaBins; k++) { cout << NormFactorTrainHisto -> GetBinCenter(k+1) << "; " << NormFactorTrainHisto -> GetBinContent(k+1)<< endl; } */ for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { tmpdouble = NormFactorTrainHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]); if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble || fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Bins defined by vector fCosThetaRangeVector " << "do not match with the bins of histogram " << NormFactorTrainHistoName << endl << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl << "NormFactorHist Bin: " << tmpdouble << endl << "Function execution will be ABORTED." << endl; return kFALSE; } } // Let's summ up all ON alpha distributions EventCounter = 0; for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull()) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Component " << SuccessfulThetaBinsVector[i]-1 << " of fThetaRangeStringVector is EMPTY" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } TString TmpHistoName = ("TrainAlphaAfterCuts"); TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; *fLog << "Getting Histogram with name " << endl << TmpHistoName << endl; if (i == 0) { CombinedAlphaTrainON = (TH1F*) rootfile.Get(TmpHistoName); NAlphaBins = CombinedAlphaTrainON -> GetNbinsX(); // Name of histogram is set CombinedAlphaTrainON -> SetName(CombinedAlphaTrainONName); // update EventCounter EventCounter += CombinedAlphaTrainON -> GetEntries(); // tmp cout << "NAlphaBins in histo CombinedAlphaTrainON: " << NAlphaBins << endl; // endtmp } else { TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName); // update EventCounter EventCounter += TmpTH1FHisto -> GetEntries(); for (Int_t j = 1; j <= NAlphaBins; j++) { // At some time alpha bins might not be the same for // the different alpha distributions. That's why // I will check it before summing the alpha histo contents tmpdouble = CombinedAlphaTrainON->GetBinCenter(j) - TmpTH1FHisto->GetBinCenter(j); tmpdouble = TMath::Abs(tmpdouble); if (tmpdouble > SmallQuantity) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Bins among the several alpha ON histograms " << "for TRAIN sample do not match" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } tmpdouble = CombinedAlphaTrainON->GetBinContent(j) + TmpTH1FHisto->GetBinContent(j); CombinedAlphaTrainON-> SetBinContent(j,tmpdouble); } // Dynamic memory allocated for TmpTH1FHisto is released //tmp cout << "Memory allocated by temporal histogram TmpTH1FHisto " << "will be released" << endl; //endtmp delete TmpTH1FHisto; } } // Entries of histogram CombinedAlphaTrainON is set to // the events counted in EventCounter tmpint = int (EventCounter); CombinedAlphaTrainON->SetEntries(tmpint); // Let's summ up all OFF alpha distributions EventCounter = 0; for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull()) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Component " << SuccessfulThetaBinsVector[i]-1 << " of fThetaRangeStringVector is EMPTY" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } TString TmpHistoName = ("TrainAlphaOFFAfterCuts"); TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; // Normalization constant for this theta bin is stored now in tmpnormfactor tmpnormfactor = NormFactorTrainHisto -> GetBinContent(SuccessfulThetaBinsVector[i]); // temp cout << "Normalization factor for bin " << SuccessfulThetaBinsVector[i] << " is " << tmpnormfactor << endl; // endtemp if (i == 0) { CombinedAlphaTrainOFF = (TH1F*) rootfile.Get(TmpHistoName); NAlphaBins = CombinedAlphaTrainOFF -> GetNbinsX(); // Name of histogram is set CombinedAlphaTrainOFF -> SetName(CombinedAlphaTrainOFFName); // Event counter is updated taking into account the // normalization factor EventCounter += tmpnormfactor * CombinedAlphaTrainOFF -> GetEntries(); // Contribution to total number of OFF events and mean normfactor // of this theta bin is computed TotalNumberOFFTrain = CombinedAlphaTrainOFF -> GetEntries(); MeanNormFactorTrain = tmpnormfactor * CombinedAlphaTrainOFF -> GetEntries(); // **** // Dimension of ErrorInCombinedAlphaTrainOFF is set ErrorInCombinedAlphaTrainOFF.Set(NAlphaBins); CombinedAlphaTrainOFFWithoutNormalization.Set(NAlphaBins); // Histogram now is normalized (to the respective ON data) // by using the normalization constants stored in // histogram NormFactorTrainHisto // Error of the normalized bins is computed and stored in // vector ErrorInCombinedAlphaTrainOFF for (Int_t j = 1; j <= NAlphaBins; j++) { // Number of events (without normalization) are // stored in vector CombinedAlphaTrainOFFWithoutNormalization CombinedAlphaTrainOFFWithoutNormalization[j-1] = CombinedAlphaTrainOFF -> GetBinContent(j); // Bin content is set tmpdouble2 = tmpnormfactor * CombinedAlphaTrainOFF -> GetBinContent(j); CombinedAlphaTrainOFF -> SetBinContent(j,tmpdouble2); // (Bin Error)^2 is computed and stored in // ErrorInCombinedAlphaTrainOFF // Bin error = Sqrt(BinContent) X Normalization Factor tmperror = tmpdouble2*tmpnormfactor; ErrorInCombinedAlphaTrainOFF[j-1] = tmperror; } } else { TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName); // Event counter is updated taking into account the // normalization factor EventCounter += tmpnormfactor * TmpTH1FHisto -> GetEntries(); // Contribution to total number of OFF events and mean normfactor // of this theta bin is computed TotalNumberOFFTrain += TmpTH1FHisto -> GetEntries(); MeanNormFactorTrain += tmpnormfactor * TmpTH1FHisto -> GetEntries(); for (Int_t j = 1; j <= NAlphaBins; j++) { // At some time alpha bins might not be the same for // the different alpha distributions. That's why // I will check it before summing the alpha histo contents tmpdouble2 = CombinedAlphaTrainOFF->GetBinCenter(j) - TmpTH1FHisto->GetBinCenter(j); tmpdouble2 = TMath::Abs(tmpdouble2); if (tmpdouble2 > SmallQuantity) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Bins among the several alpha OFF histograms " << "for TRAIN sample do not match" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } // Number of events (without normalization) are // added in vector CombinedAlphaTrainOFFWithoutNormalization CombinedAlphaTrainOFFWithoutNormalization[j-1] += TmpTH1FHisto -> GetBinContent(j); // Histogram bin contents must be normalized (to the respective ON data) // by using the normalization constants stored in // histogram NormFactorTrainHisto before being added // to CombinedAlphaTrainOFF tmpdouble2 = CombinedAlphaTrainOFF->GetBinContent(j) + (TmpTH1FHisto->GetBinContent(j) * tmpnormfactor); CombinedAlphaTrainOFF-> SetBinContent(j,tmpdouble2); // (Bin Error)^2 is computed and added to // ErrorInCombinedAlphaTrainOFF // Bin error = Sqrt(BinContent) X Normalization Factor tmperror = TmpTH1FHisto->GetBinContent(j) * tmpnormfactor; ErrorInCombinedAlphaTrainOFF[j-1] += tmperror; } // Dynamic memory allocated for TmpTH1FHisto is released delete TmpTH1FHisto; } } // Mean Normalization factor is computed for the combined OFF histogram. // This factor will be used when computing the significance via // MHFindSignificance::FindSigmaONOFF(). // The bin contents (and errors) of teh combined OFF TRAIN histo will be // normalized (divided) by teh mean norm constant, so that // the histogram for OFF data given as argument in FindSigmaONOFF // is equivalent to a non-normalized histogram. MeanNormFactorTrain = MeanNormFactorTrain/TotalNumberOFFTrain; // Sqrt(Contents) is performed in ErrorInCombinedAlphaTrainOFF, // in order to compute the real error in the bin content of histogram // CombinedAlphaTrainOFF. // Then error and contents (normalized with mean normfctor) are // set to histogram bin for (Int_t j = 1; j <= NAlphaBins; j++) { tmpdouble2 = (CombinedAlphaTrainOFF -> GetBinContent(j))/MeanNormFactorTrain; CombinedAlphaTrainOFF -> SetBinContent(j,tmpdouble2); ErrorInCombinedAlphaTrainOFF[j-1] = TMath::Sqrt(ErrorInCombinedAlphaTrainOFF[j-1])/MeanNormFactorTrain; // Proper error treatment when adding quantities to histogram bin CombinedAlphaTrainOFF -> SetBinError(j,ErrorInCombinedAlphaTrainOFF[j-1]); // ****** TMP ************** /* // Error computation assuming a given (no normalized) bin content tmpdouble = CombinedAlphaTrainOFF -> GetBinContent(j); if (CombinedAlphaTrainOFFWithoutNormalization[j-1] < SmallQuantity) {tmpnormfactor = 0;} else {tmpnormfactor = tmpdouble/CombinedAlphaTrainOFFWithoutNormalization[j-1];} tmperror = TMath::Sqrt(CombinedAlphaTrainOFFWithoutNormalization[j-1]) * tmpnormfactor; CombinedAlphaTrainOFF -> SetBinError(j,tmperror); */ // ****** END TMP ************** // tmp /* cout << CombinedAlphaTrainOFF -> GetBinContent(j) << " +/- " << CombinedAlphaTrainOFF -> GetBinError(j) << " Number Events without normalization: " << CombinedAlphaTrainOFFWithoutNormalization[j-1] << endl; */ // endtmp } // Entries of histogram CombinedAlphaTrainON is set to // the events counted during histogram filling EventCounter = EventCounter/MeanNormFactorTrain; tmpint = int (EventCounter); CombinedAlphaTrainOFF->SetEntries(tmpint); // Check that the number of events in the combined normalized histogram is // that of the initial sample /* cout << "SILLY INFO FOR TRAIN SAMPLE: Total number of events Before nomralization and " << " re-re normalized" << endl << TotalNumberOFFTrain << " - " << EventCounter << endl; */ } // ********* COMBINED HISTOGRAMS FOR TRAIN SAMPLE ALREADY FILLED ********* /// // ********* HISTOGRAMS FOR TEST SAMPLE WILL BE FILLED ********* /// if (CombineTestData) { *fLog << endl << "****************************************************************************" << endl << "Combining data for the TEST sample " << endl << "****************************************************************************" << endl << endl; // Retrieving NormFactorTestHisto from root file. NormFactorTestHisto = (TH1F*) rootfile.Get(NormFactorTestHistoName); // NormFactorTestHisto-> DrawCopy(); // Check that bins in this histo correspond with the // CosTheta intervals defined by vector fCosThetaRangeVector if (NormFactorTestHisto -> GetNbinsX() != NThetaBins) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Number of theta bins defined by vector fCosThetaRangeVector (" << NThetaBins << ") does not match with the bins in histogram " << NormFactorTestHistoName << "(" << NormFactorTestHisto -> GetNbinsX() << ")" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } // histo test /* cout << "NORMFACTORTRAINHIST TEST: BinCenter and Value" << endl; for (Int_t k = 0; k < NThetaBins; k++) { cout << NormFactorTestHisto -> GetBinCenter(k+1) << "; " << NormFactorTestHisto -> GetBinContent(k+1)<< endl; } */ for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { tmpdouble = NormFactorTestHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]); if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble || fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Bins defined by vector fCosThetaRangeVector " << "do not match with the bins of histogram " << NormFactorTestHistoName << endl << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl << "NormFactorHist Bin: " << tmpdouble << endl << "Function execution will be ABORTED." << endl; return kFALSE; } } // Let's summ up all ON alpha distributions // Event counter (counts events used in the histogram filling) is // initialized to zero EventCounter = 0; for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull()) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Component " << SuccessfulThetaBinsVector[i]-1 << " of fThetaRangeStringVector is EMPTY" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } TString TmpHistoName = ("TestAlphaAfterCuts"); TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; *fLog << "Getting Histogram with name " << endl << TmpHistoName << endl; if (i == 0) { CombinedAlphaTestON = (TH1F*) rootfile.Get(TmpHistoName); NAlphaBins = CombinedAlphaTestON -> GetNbinsX(); // Name of histogram is set CombinedAlphaTestON -> SetName(CombinedAlphaTestONName); // EventCounter is updated EventCounter += CombinedAlphaTestON -> GetEntries(); // tmp cout << "NAlphaBins in histo CombinedAlphaTestON: " << NAlphaBins << endl; // endtmp } else { TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName); // EventCounter is updated EventCounter += TmpTH1FHisto -> GetEntries(); for (Int_t j = 1; j <= NAlphaBins; j++) { // At some time alpha bins might not be the same for // the different alpha distributions. That's why // I will check it before summing the alpha histo contents tmpdouble = CombinedAlphaTestON->GetBinCenter(j) - TmpTH1FHisto->GetBinCenter(j); tmpdouble = TMath::Abs(tmpdouble); if (tmpdouble > SmallQuantity) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Bins among the several alpha ON histograms " << "for TEST sample do not match" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } tmpdouble = CombinedAlphaTestON->GetBinContent(j) + TmpTH1FHisto->GetBinContent(j); CombinedAlphaTestON-> SetBinContent(j,tmpdouble); } // Dynamic memory allocated for TmpTH1FHisto is released //tmp cout << "Memory allocated by temporal histogram TmpTH1FHisto " << "will be released" << endl; //endtmp delete TmpTH1FHisto; } } // Entries of histogram CombinedAlphaTestON is set to // the events counted during histogram filling tmpint = int (EventCounter); CombinedAlphaTestON->SetEntries(tmpint); // Let's summ up all OFF alpha distributions EventCounter = 0; for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++) { if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull()) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Component " << SuccessfulThetaBinsVector[i]-1 << " of fThetaRangeStringVector is EMPTY" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } TString TmpHistoName = ("TestAlphaOFFAfterCuts"); TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1]; // Normalization constant for this theta bin is stored now in tmpnormfactor tmpnormfactor = NormFactorTestHisto -> GetBinContent(SuccessfulThetaBinsVector[i]); // temp cout << "Normalization factor for bin " << SuccessfulThetaBinsVector[i] << " is " << tmpnormfactor << endl; // endtemp if (i == 0) { CombinedAlphaTestOFF = (TH1F*) rootfile.Get(TmpHistoName); NAlphaBins = CombinedAlphaTestOFF -> GetNbinsX(); // Name of histogram is set CombinedAlphaTestOFF -> SetName(CombinedAlphaTestOFFName); // EventCounter is updated taking into account // normalization factor EventCounter += tmpnormfactor*CombinedAlphaTestOFF -> GetEntries(); // Contribution to total number of OFF events and mean normfactor // of this theta bin is computed TotalNumberOFFTest = CombinedAlphaTestOFF -> GetEntries(); MeanNormFactorTest = tmpnormfactor * CombinedAlphaTestOFF -> GetEntries(); // Dimension of ErrorInCombinedAlphaTrainOFF is set ErrorInCombinedAlphaTestOFF.Set(NAlphaBins); CombinedAlphaTestOFFWithoutNormalization.Set(NAlphaBins); // Histogram now is normalized (to the respective ON data) // by using the normalization constants stored in // histogram NormFactorTestHisto for (Int_t j = 1; j <= NAlphaBins; j++) { // Vector containing number of events without // normalization is filled CombinedAlphaTestOFFWithoutNormalization[j-1] = CombinedAlphaTestOFF -> GetBinContent(j); tmpdouble2 = tmpnormfactor * CombinedAlphaTestOFF -> GetBinContent(j); CombinedAlphaTestOFF -> SetBinContent(j,tmpdouble2); // (Bin Error)^2 is computed and stored in // ErrorInCombinedAlphaTestOFF // Bin error = Sqrt(BinContent) X Normalization Factor tmperror = tmpdouble2*tmpnormfactor; ErrorInCombinedAlphaTestOFF[j-1] = tmperror; } } else { TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName); // EventCounter is updated taking into account // normalization factor EventCounter += tmpnormfactor* TmpTH1FHisto-> GetEntries(); // Contribution to total number of OFF events and mean normfactor // of this theta bin is computed TotalNumberOFFTest += TmpTH1FHisto -> GetEntries(); MeanNormFactorTest += tmpnormfactor * TmpTH1FHisto -> GetEntries(); for (Int_t j = 1; j <= NAlphaBins; j++) { // At some time alpha bins might not be the same for // the different alpha distributions. That's why // I will check it before summing the alpha histo contents tmpdouble2 = CombinedAlphaTestOFF->GetBinCenter(j) - TmpTH1FHisto->GetBinCenter(j); tmpdouble2 = TMath::Abs(tmpdouble2); if (tmpdouble2 > SmallQuantity) { *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl << "Bins among the several alpha OFF histograms " << "for TRAIN sample do not match" << endl << "Function execution will be ABORTED." << endl; return kFALSE; } // Vector containing number of events without // normalization is updated // Vector containing number of events without // normalization is filled CombinedAlphaTestOFFWithoutNormalization[j-1] += TmpTH1FHisto -> GetBinContent(j); // Histogram bin contents must be normalized (to the respective ON data) // by using the normalization constants stored in // histogram NormFactorTestHisto before being added // to CombinedAlphaTestOFF tmpdouble2 = CombinedAlphaTestOFF->GetBinContent(j) + (TmpTH1FHisto->GetBinContent(j) * tmpnormfactor); CombinedAlphaTestOFF-> SetBinContent(j,tmpdouble2); // (Bin Error)^2 is computed and added to // ErrorInCombinedAlphaTestOFF // Bin error = Sqrt(BinContent) X Normalization Factor tmperror = TmpTH1FHisto->GetBinContent(j) * tmpnormfactor; ErrorInCombinedAlphaTestOFF[j-1] += tmperror; } // Dynamic memory allocated for TmpTH1FHisto is released delete TmpTH1FHisto; } } // Mean Normalization factor is computed for the combined OFF histogram. // This factor will be used when computing the significance via // MHFindSignificance::FindSigmaONOFF(). // The bin contents (and errors) of teh combined OFF TEST histo will be // normalized (divided) by teh mean norm constant, so that // the histogram for OFF data given as argument in FindSigmaONOFF // is equivalent to a non-normalized histogram. MeanNormFactorTest = MeanNormFactorTest/TotalNumberOFFTest; // Sqrt(Contents) is performed in ErrorInCombinedAlphaTestOFF, // in order to compute the real error in the bin content of histogram // CombinedAlphaTestOFF. // Then error and contents (normalized with mean normfctor) are // set to histogram bin for (Int_t j = 1; j <= NAlphaBins; j++) { tmpdouble2 = (CombinedAlphaTestOFF -> GetBinContent(j))/MeanNormFactorTest; CombinedAlphaTestOFF -> SetBinContent(j,tmpdouble2); ErrorInCombinedAlphaTestOFF[j-1] = TMath::Sqrt(ErrorInCombinedAlphaTestOFF[j-1])/MeanNormFactorTest; // Proper error treatment for the bin contents CombinedAlphaTestOFF -> SetBinError(j,ErrorInCombinedAlphaTestOFF[j-1]); // ****** TMP ************** /* // Error computation assuming a given (no normalized) bin content tmpdouble = CombinedAlphaTestOFF -> GetBinContent(j); if (CombinedAlphaTestOFFWithoutNormalization[j-1] < SmallQuantity) {tmpnormfactor = 0;} else {tmpnormfactor = tmpdouble/CombinedAlphaTestOFFWithoutNormalization[j-1];} tmperror = (TMath::Sqrt(CombinedAlphaTestOFFWithoutNormalization[j-1])) * tmpnormfactor; CombinedAlphaTestOFF -> SetBinError(j,tmperror); */ // ****** ENDTMP ************** // tmp // cout << CombinedAlphaTestOFF -> GetBinContent(j) << " +/- " // << CombinedAlphaTestOFF -> GetBinError(j) << endl; // endtmp } // Entries of histogram CombinedAlphaTestOFF is set to // the events counted during histogram filling EventCounter = EventCounter/MeanNormFactorTest; tmpint = int(EventCounter); CombinedAlphaTestOFF->SetEntries(tmpint); // Check that the number of events in the combined normalized histogram is // that of the initial sample /* cout << "SILLY INFO FOR TEST SAMPLE: Total number of events Before nomralization and " << " re-re normalized" << endl << TotalNumberOFFTest << " - " << EventCounter << endl; */ /* // Sumw2() test cout << "******** TEST ********* " << endl; cout << "Sumw2() is executed: " << endl; CombinedAlphaTestOFF -> Sumw2(); for (Int_t j = 1; j <= NAlphaBins; j++) { // tmp cout << CombinedAlphaTestOFF -> GetBinContent(j) << " +/- " << CombinedAlphaTestOFF -> GetBinError(j) << endl; // endtmp } */ } // ********* COMBINED HISTOGRAMS FOR TEST SAMPLE ALREADY FILLED ********* /// /* CombinedAlphaTrainON-> DrawCopy(); gPad -> SaveAs("tmpONTrain.ps"); CombinedAlphaTrainOFF-> DrawCopy(); gPad -> SaveAs("tmpOFFTrain.ps"); CombinedAlphaTestON-> DrawCopy(); gPad -> SaveAs("tmpONTest.ps"); CombinedAlphaTestOFF-> DrawCopy(); gPad -> SaveAs("tmpOFFTest.ps"); */ // **** COMPUTATION OF SIGNIFICANCES FOR COMBINED ALPHA HISTOGRAMS ****** /// if (SuccessfulThetaBinsVector.GetSize() >= 1) { *fLog << endl << "**************************************************************************************" << endl << "Computation of excess events and singnificance for the combined alpha histograms " << endl << "**************************************************************************************" << endl << endl; // There is, at least, ONE theta bin for which optimization of // supercuts was SUCCESSFUL, and thus, it is worth to // compute significance of signal // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF // Maximum value of alpha below which signal is expected const Double_t alphasig = fAlphaSig; // Minimum value of alpha for bkg region in ON data const Double_t alphabkgmin = fAlphaBkgMin; // Maximum value of alpha for bkg region in ON data const Double_t alphabkgmax = fAlphaBkgMax; const Bool_t drawpoly = kTRUE; const Bool_t fitgauss = kTRUE; const Bool_t print = kTRUE; const Bool_t saveplots = kTRUE; // Save final alpha plot for ON and OFF // (with Nex and significance) in Psfile MHFindSignificanceONOFF findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); if (CombineTrainData) { // Name of psfile where Alpha plot will be stored TString psfilename = (fPathForFiles); psfilename += ("TRAINSampleCombinedThetaBins"); findsig.FindSigmaONOFF(CombinedAlphaTrainON, CombinedAlphaTrainOFF, MeanNormFactorTrain, alphabkgmin, alphabkgmax, fDegreeON, fDegreeOFF, alphasig, drawpoly, fitgauss, print, saveplots, psfilename); fOverallSigmaLiMaTrain = double(findsig.GetSignificance()); fOverallNexTrain = double(findsig.GetNexONOFF()); *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl << "After performing the combined analysis for TRAIN sample it was found " << "the following Nex and Significance: " << fOverallNexTrain << ", " << fOverallSigmaLiMaTrain << endl; } if (CombineTestData) { // Name of psfile where Alpha plot will be stored TString psfilename = (fPathForFiles); psfilename += ("TESTSampleCombinedThetaBins"); findsig.FindSigmaONOFF(CombinedAlphaTestON, CombinedAlphaTestOFF, MeanNormFactorTest, alphabkgmin, alphabkgmax, fDegreeON, fDegreeOFF, alphasig, drawpoly, fitgauss, print, saveplots, psfilename); fOverallSigmaLiMaTest = double(findsig.GetSignificance()); fOverallNexTest = double(findsig.GetNexONOFF()); *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl << "After performing the combined analysis for TEST sample it was found " << "the following Nex and Significance: " << fOverallNexTest << ", " << fOverallSigmaLiMaTest << endl; } *fLog << endl << "**************************************************************************************" << endl << "Alpha distributions (Train/Test and ON/OFF) for the combined " << endl << " theta bins are stored in the root file " << endl << fAlphaDistributionsRootFilename << endl << "**************************************************************************************" << endl << endl; TFile rootfiletowrite(fAlphaDistributionsRootFilename, "UPDATE"); // Histograms that will contain the Norm factors for the // combined TRAIN/TEST theta bins are created TString CombinedNormFactorTrainName = ("TrainSampleNormFactorCombinedThetaBins"); TH1F CombinedNormFactorTrain (CombinedNormFactorTrainName, CombinedNormFactorTrainName, 1, 0.0, 1.0); TString CombinedNormFactorTestName = ("TestSampleNormFactorCombinedThetaBins"); TH1F CombinedNormFactorTest (CombinedNormFactorTestName, CombinedNormFactorTestName, 1, 0.0, 1.0); if (CombineTrainData) { // Histogram that will contain the Norm factors for the // combined TRAIN theta bins are filled CombinedNormFactorTrain.Fill(0.5, MeanNormFactorTrain); // Train histos are written to root file CombinedAlphaTrainON -> Write(); CombinedAlphaTrainOFF -> Write(); CombinedNormFactorTrain.Write(); } if (CombineTestData) { // Histogram that will contain the Norm factors for the // combined TEST theta bins are filled CombinedNormFactorTest.Fill(0.5, MeanNormFactorTest); // Test histos are written to root file CombinedAlphaTestON -> Write(); CombinedAlphaTestOFF -> Write(); CombinedNormFactorTest.Write(); } // Root file is closed rootfiletowrite.Close(); } return kTRUE; }