/* ======================================================================== *\ ! ! * ! * 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): Wolfgang Wittek 9/2003 ! David Paneque 11/2003 ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MFindSupercutsONOFF // // // // Class for optimizing the parameters of the supercuts // // Using ON and OFF data // // // // // ///////////////////////////////////////////////////////////////////////////// #include "MFindSupercutsONOFF.h" #include // fabs #include #include #include #include #include #include #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" #include "MFRandomSplit.h" #include "MH3.h" #include "MHCT1Supercuts.h" #include "MHFindSignificance.h" // To be removed at some point... #include "MHFindSignificanceONOFF.h" #include "MTSupercutsApplied.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; // Function that computes the normalization factor using COUNTED events in alpha histogram // for ON data and alpha histogram for OFF data; // in alpha region defined by AlphaBkgMin and AlphaBkgMax // It is defined outside the class MFindSupercutsONOFF so that it can be used // in function fcnSupercuts static Double_t ComputeNormFactorFromAlphaBkg(TH1 *histON, TH1 *histOFF, Double_t AlphaBkgMin, Double_t AlphaBkgMax) { Double_t NormFactor = 0.0; Double_t ONEvents = 0.0; Double_t OFFEvents = 0.0; const Double_t SmallQuantity = 0.01; Double_t xlo = 0.0; Double_t xup = 0.0; Double_t width = 0.0; Int_t BinCounterOFF = 0; Int_t BinCounterON = 0; // I make a copy of the histograms so that nothing happens to the // histograms used in argument TH1* HistON = (TH1*) histON->Clone(); TH1* HistOFF = (TH1*) histOFF->Clone(); if ( !HistON ) { gLog << "ComputeNormFactorFromAlphaBkg; " << endl << "Clone of ON histogram could not be generated" << endl; return 0.0; } if ( !HistOFF ) { gLog << "ComputeNormFactorFromAlphaBkg; " << endl << " Clone of OFF histogram could not be generated" << endl; return 0.0; } // Calculate the number of OFF events in the Bkg region // defined by AlphaBkgMin and AlphaBkgMax // ___________________________________________________ Int_t nbinsOFF = HistOFF -> GetNbinsX(); Double_t binwidthOFF = HistOFF -> GetBinWidth(1); for (Int_t i=1; i<=nbinsOFF; i++) { xlo = HistOFF->GetBinLowEdge(i); xup = HistOFF->GetBinLowEdge(i+1); // bin must be completely contained in the bkg region if ( xlo >= (AlphaBkgMin-SmallQuantity) && xup <= (AlphaBkgMax+SmallQuantity) ) { width = fabs(xup-xlo); if (fabs(width-binwidthOFF) > SmallQuantity) { gLog << "ComputeNormFactorFromAlphaBkg; " << endl << " HistOFF has variable binning, which is not allowed" << endl; return 0.0; } BinCounterOFF++; OFFEvents += HistOFF->GetBinContent(i); } } // Calculate the number of ON events in the Bkg region // defined by AlphaBkgMin and AlphaBkgMax // ___________________________________________________ Int_t nbinsON = HistON -> GetNbinsX(); Double_t binwidthON = HistON -> GetBinWidth(1); for (Int_t i=1; i<=nbinsON; i++) { xlo = HistON->GetBinLowEdge(i); xup = HistON->GetBinLowEdge(i+1); // bin must be completely contained in the bkg region if ( xlo >= (AlphaBkgMin-SmallQuantity) && xup <= (AlphaBkgMax+SmallQuantity) ) { width = fabs(xup-xlo); if (fabs(width-binwidthON) > SmallQuantity) { gLog << "ComputeNormFactorFromAlphaBkg; " << endl << " HistON has variable binning, which is not allowed" << endl; return 0.0; } BinCounterON++; ONEvents += HistON->GetBinContent(i); } } // NormFactor is computed if (ONEvents < SmallQuantity || OFFEvents < SmallQuantity) { gLog << "ComputeNormFactorFromAlphaBkg; " << endl << "ONEvents or OFFEvents computed in bkg region are < " << SmallQuantity << endl; return 0.0; } NormFactor = ONEvents/OFFEvents; Double_t error = 1/ONEvents + 1/OFFEvents; error = TMath::Sqrt(error); error = error * NormFactor; // tmp info gLog << "ComputeNormFactorFromAlphaBkg;" << endl << "ON Events in bkg region = " << ONEvents << " (" << BinCounterON << " bins)" << " , OFF Events in bkg region = " << OFFEvents << " (" << BinCounterOFF << " bins)" << endl <<"NormFactor computed from bkg region = " << NormFactor << " +/- " << error << endl; // end temp return NormFactor; } //------------------------------------------------------------------------ // // fcnSupercuts // // - calculates the quantity to be minimized (using TMinuit) // // - the quantity to be minimized is (-1)*significance of the gamma signal // in the alpha distribution (after cuts) // // - the parameters to be varied in the minimization are the cut parameters // (par) // static void fcnSupercuts(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) { //cout << "entry fcnSupercuts" << endl; //------------------------------------------------------------- // save pointer to the MINUIT object for optimizing the supercuts // because it will be overwritten // when fitting the alpha distribution in MHFindSignificance TMinuit *savePointer = gMinuit; //------------------------------------------------------------- MEvtLoop* evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit(); // Event loops for ON and OFF data are recovered. MEvtLoop* ONDataevtloopfcn = &evtloopfcn[0]; MEvtLoop* OFFDataevtloopfcn = &evtloopfcn[1]; // Parameter list from event loops for ON and OFF data are recovered MParList *ONDataplistfcn = (MParList*) ONDataevtloopfcn->GetParList(); MParList *OFFDataplistfcn = (MParList*) OFFDataevtloopfcn->GetParList(); // MTaskList *ONDataTasklistfcn = (MTaskList*) ONDataevtloopfcn->GetTaskList(); // MTaskList *OFFDataTasklistfcn = (MTaskList*) OFFDataevtloopfcn->GetTaskList(); // Container for supercuts is retrieved from ONDataplistfcn. // NOTE: The same supercuts parameter container is used in OFFDataplistfcn MSupercuts *super = (MSupercuts*) ONDataplistfcn->FindObject("MSupercuts"); if (!super) { gLog << "fcnSupercuts : MSupercuts object '" << "MSupercuts" << "' not found... aborting" << endl; return; } // Normalization factor for train sample (Train ON before cuts/ train OFF before cuts) // is retrieved from the ONDataplistfcn Double_t NormalizationFactorTrain = 0.0; MDataValue* NormFactorContainer = (MDataValue*) ONDataplistfcn->FindObject("NormFactorTrain"); NormalizationFactorTrain = NormFactorContainer -> GetValue(); gLog << "fcnSupercuts : Normalization factor retrieved from ONDataplistfcn: " << endl << "NormalizationFactorTrain = " << NormalizationFactorTrain << endl; Double_t AlphaSig = 0.0; MDataValue* AlphaSigContainer = (MDataValue*) ONDataplistfcn->FindObject("AlphaSigValue"); AlphaSig = AlphaSigContainer -> GetValue(); Double_t AlphaBkgMin = 0.0; MDataValue* AlphaBkgMinContainer = (MDataValue*) ONDataplistfcn->FindObject("AlphaBkgMinValue"); AlphaBkgMin = AlphaBkgMinContainer -> GetValue(); Double_t AlphaBkgMax = 0.0; MDataValue* AlphaBkgMaxContainer = (MDataValue*) ONDataplistfcn->FindObject("AlphaBkgMaxValue"); AlphaBkgMax = AlphaBkgMaxContainer -> GetValue(); gLog << "fcnSupercuts : AlphaSig and AlphaBkgMin-AlphaBkgMax retrieved from ONDataplistfcn: " << endl << "AlphaSig = " << AlphaSig << "; AlphaBkgMin-AlphaBkgMax = " << AlphaBkgMin << "-" << AlphaBkgMax << endl; // Variable fUseFittedQuantities is retrieved from ONDataplistfcn Bool_t UseFittedQuantities; MDataValue* UseFittedQuantitiesContainer = (MDataValue*) ONDataplistfcn->FindObject("UseFittedQuantitiesValue"); UseFittedQuantities = bool(UseFittedQuantitiesContainer -> GetValue()); if (UseFittedQuantities) { gLog << "fcnSupercuts : UseFittedQuantities variable set to kTRUE" << endl; } else { gLog << "fcnSupercuts : UseFittedQuantities variable set to kFALSE" << endl; } Bool_t UseNormFactorFromAlphaBkg; MDataValue* UseNormFactorFromAlphaBkgContainer = (MDataValue*) ONDataplistfcn -> FindObject("UseNormFactorFromAlphaBkgValue"); UseNormFactorFromAlphaBkg = bool(UseNormFactorFromAlphaBkgContainer -> GetValue()); if (UseNormFactorFromAlphaBkg) { gLog << "fcnSupercuts : UseNormFactorFromAlphaBkg variable set to kTRUE" << endl; } else { gLog << "fcnSupercuts : UseNormFactorFromAlphaBkg variable set to kFALSE" << endl; } // // transfer current parameter values to MSupercuts // // Attention : npar is the number of variable parameters // not the total number of parameters // Double_t fMin, fEdm, fErrdef; Int_t fNpari, fNparx, fIstat; gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); super->SetParameters(TArrayD(fNparx, par)); //$$$$$$$$$$$$$$$$$$$$$ // for testing // TArrayD checkparameters = super->GetParameters(); //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", " // << fNparx << endl; //gLog << "fcnsupercuts : i, par, checkparameters =" << endl; //for (Int_t i=0; iEventloop()) { gLog << "fcnsupercuts : ONDataevtloopfcn->Eventloop() failed" << endl; } // Somehow (??) I can not use the function MTaskList::PrintStatistics... // it produces a segmentation fault... //ONDataTasklistfcn->PrintStatistics(0, kFALSE); MH3* alpha = (MH3*)ONDataplistfcn->FindObject("AlphaFcn", "MH3"); if (!alpha) return; TH1 &alphaHist = alpha->GetHist(); alphaHist.SetName("alpha-fcnSupercuts"); // // plot alpha for OFF data with the current cuts // if(!OFFDataevtloopfcn->Eventloop()) { gLog << "fcnsupercuts : OFFDataevtloopfcn->Eventloop() failed" << endl; } // Somehow (??) I can not use the function MTaskList::PrintStatistics... // it produces a segmentation fault... //OFFDataTasklistfcn->PrintStatistics(0, kFALSE); MH3* alphaOFF = (MH3*) OFFDataplistfcn->FindObject("AlphaOFFFcn", "MH3"); if (!alphaOFF) return; TH1 &alphaHistOFF = alphaOFF->GetHist(); alphaHistOFF.SetName("alphaOFF-fcnSupercuts"); if (UseNormFactorFromAlphaBkg) { // Normalization factor computed using alpha bkg region Double_t NewNormFactor = ComputeNormFactorFromAlphaBkg(&alphaHist, &alphaHistOFF, AlphaBkgMin, AlphaBkgMax); gLog << "Normalization factor computed from alpha plot (after cuts) " << endl << "using counted number of ON and OFF events in alpha region " << endl << "defined by range " << AlphaBkgMin << "-" << AlphaBkgMax << endl << "Normalization factor = " << NewNormFactor << endl; gLog << "Normalization factor used is the one computed in bkg region; " << endl << "i.e. " << NewNormFactor << " instead of " << NormalizationFactorTrain << endl; NormalizationFactorTrain = NewNormFactor; } //------------------------------------------- // set Minuit pointer to zero in order not to destroy the TMinuit // object for optimizing the supercuts gMinuit = NULL; //================================================================= // fit alpha distribution to get the number of excess events and // calculate significance of gamma signal in the alpha plot const Double_t alphasig = AlphaSig; const Double_t alphamin = AlphaBkgMin; // alpha min for bkg region in ON data const Double_t alphamax = AlphaBkgMax; // alpha max for bkg region in ON data const Int_t degreeON = 2; // degree of polynomial used to fit ON data in Bkg region const Int_t degreeOFF = 2; // degree of polynomial used to fit OFF data in ALL region Bool_t drawpoly; Bool_t fitgauss; Bool_t saveplots; if (iflag == 3) {// Even though minimization finished successfully, I will NOT produce // final plots now... i'll do it later via the function // "MFindSupercutsONOFF::OutputOptimizationOnTrainSample()" drawpoly = kFALSE; fitgauss = kFALSE; saveplots = kFALSE; } else { drawpoly = kFALSE; fitgauss = kFALSE; saveplots = kFALSE; } const Bool_t print = kTRUE; MHFindSignificanceONOFF findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); findsig.SetUseFittedQuantities(UseFittedQuantities); // TPostScript* DummyPs = new TPostScript("dummy.ps"); TString DummyPs = ("Dummy"); const Bool_t rc = findsig.FindSigmaONOFF(&alphaHist,&alphaHistOFF, NormalizationFactorTrain, alphamin, alphamax, degreeON, degreeOFF, alphasig, drawpoly, fitgauss, print, saveplots, DummyPs); //DummyPs -> Close(); //delete DummyPs; //================================================================= // reset gMinuit to the MINUIT object for optimizing the supercuts gMinuit = savePointer; //------------------------------------------- if (!rc) { gLog << "fcnSupercuts : FindSigmaONOFF() failed" << endl; f = 1.e10; return; } /* // plot some quantities during the optimization MHCT1Supercuts *plotsuper = (MHCT1Supercuts*)ONDataplistfcn->FindObject("MHCT1Supercuts"); if (plotsuper) plotsuper->Fill(&findsig); */ //------------------------ // get significance const Double_t significance = findsig.GetSignificance(); f = significance>0 ? -significance : 0; gLog << " ///****************************************** ///" << endl; gLog << "Significance (Li&Ma)is now: " << f << endl; gLog << " ///****************************************** ///" << endl; //------------------------ // optimize signal/background ratio //Double_t ratio = findsig.GetNbg()>0.0 ? // findsig.GetNex()/findsig.GetNbg() : 0.0; //f = -ratio; //------------------------------------------- // final calculations //if (iflag == 3) //{ // //} //------------------------------------------------------------- } // -------------------------------------------------------------------------- // // Default constructor. // MFindSupercutsONOFF::MFindSupercutsONOFF(const char *name, const char *title) : fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL) { fName = name ? name : "MFindSupercutsONOFF"; fTitle = title ? title : "Optimizer of the supercuts"; //--------------------------- // camera geometry is needed for conversion mm ==> degree //fCam = new MGeomCamCT1Daniel; // fCam = new MGeomCamCT1; fCam = new MGeomCamMagic; // matrices to contain the training/test samples fMatrixTrain = new MHMatrix("MatrixTrain"); fMatrixTest = new MHMatrix("MatrixTest"); fMatrixTrainOFF = new MHMatrix("MatrixTrainOFF"); fMatrixTestOFF = new MHMatrix("MatrixTestOFF"); // objects of MSupercutsCalcONOFF to which these matrices are attached fCalcHadTrain = new MSupercutsCalcONOFF("SupercutsCalcTrain"); fCalcHadTest = new MSupercutsCalcONOFF("SupercutsCalcTest"); fCalcHadTrainOFF = new MSupercutsCalcONOFF("SupercutsCalcTrainOFF"); fCalcHadTestOFF = new MSupercutsCalcONOFF("SupercutsCalcTestOFF"); // Define columns of matrices fCalcHadTrain->InitMapping(fMatrixTrain); fCalcHadTest->InitMapping(fMatrixTest); fCalcHadTrainOFF->InitMapping(fMatrixTrainOFF); fCalcHadTestOFF->InitMapping(fMatrixTestOFF); // For the time being weights method is not implemented //UseWeights = kFALSE; // Normalization factors are initialized to -1 // Functions determining Number of excess events will not work // with negative norm factors; // This will ensure that Norm factors are computed before running // DetExcessONOFF() function. fNormFactorTrain = -1.0; fNormFactorTest = -1.0; // SigmaLiMa and Nex member variables are initialized to 0 fSigmaLiMaTrain = 0.0; fSigmaLiMaTest = 0.0; fNexTrain = 0.0; fNexTest = 0.0; // Cuts (low and up) in variable MPointingPos.fZd // (at some point real theta) // The default is not cut, i.e. all values (0-1) are taken fThetaMin = 0; // in miliradians // FIXME: change name fThetaMax = 90; // new: in deg; old: (1570) in miliradians // FIXME: change name 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; // By default, bining for alpha plots is set to fNAlphaBins = 54; fAlphaBinLow = -12; fAlphaBinUp = 96; fNormFactorFromAlphaBkg = kTRUE; // By default normalization factor is // computed using bkg region in alpha histograms (after cuts) // fActualCosThetaBinCenter = 0; 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 // 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. // Default value is kTRUE fUseFittedQuantities = kTRUE; // Boolean variable that allows the user to set some limits to the // some of the Minuit parameters. For the time being, only the limits // for the parameters which do NOT depend in size, dist and theta are set, // i.e. static limits. The value of this boolean variable is set in the // constructor of the class. fSetLimitsToSomeMinuitParams = kTRUE; // Limits for the Minuit parameters. For the time being the values are set in the constructor // of the class. One MUST be very careful to set limits such that the expected final values // (optimized values) are far away from the limits. // Values in degrees fMinuitDistUPUpperLimit = 2.0; fMinuitDistUPLowerLimit = 0.5; fMinuitLengthUPUpperLimit = 0.8; fMinuitLengthUPLowerLimit = 0.0; fMinuitWidthUPUpperLimit = 0.5; fMinuitWidthUPLowerLimit = 0.0; fMinuitLeakage1UPUpperLimit = 1.5; fMinuitLeakage1UPLowerLimit = -0.5; fMinuitDistLOWUpperLimit = 1.0; fMinuitDistLOWLowerLimit = 0.0; fMinuitLengthLOWUpperLimit = 0.5; fMinuitLengthLOWLowerLimit = -0.3; fMinuitWidthLOWUpperLimit = 0.4; fMinuitWidthLOWLowerLimit = -0.3; // 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; // Set wether to use or not hillas dist fUseDist = kTRUE; 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 (standard value) fPsFilename = NULL; fPsFilename2 = NULL; //////////////////////////////////////////////////// // TMP // There are quite some problems during the data preprocessing. // For the time being, I will add some cuts to the functions // DefineTrainTestMatrixThetaRange and for OFF, so that I can // make a kind of preprocess on my own. This allows me // to make a very silly preprocess with wolfgangs macro, which // might be free of corrupted data, and then I can do on my own. fSizeCutLow = 0.1; // To prevent empty events fSizeCutUp = 10000000; // Angular cuts are converted to mm, which // are the units of the preprocessed data.... // Angular cuts not yet implemented ... Double_t ConvMMToDeg = 0.00337034; fDistCutLow = 0.4/ConvMMToDeg; fDistCutUp = 1.5/ConvMMToDeg; fLengthCutLow = 0.1/ConvMMToDeg; fLengthCutUp = 1/ConvMMToDeg; fWidthCutLow = 0.07/ConvMMToDeg; fWidthCutUp = 1/ConvMMToDeg; // ENDTMP //////////////////////////////////////////////////// } // -------------------------------------------------------------------------- // // Default destructor. // MFindSupercutsONOFF::~MFindSupercutsONOFF() { *fLog << "destructor of MFindSupercutsONOFF is called" << endl; fPsFilename = NULL; fPsFilename2 = NULL; delete fCam; delete fMatrixTrain; delete fMatrixTest; delete fCalcHadTrain; delete fCalcHadTest; delete fMatrixTrainOFF; delete fMatrixTestOFF; delete fCalcHadTrainOFF; delete fCalcHadTestOFF; *fLog << "destructor of MFindSupercutsONOFF 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 MFindSupercutsONOFF::SetPostScriptFile (TPostScript* PsFile) { fPsFilename = PsFile; *fLog << "MFindSupercutsONOFF : Results (alpha distributions with excess and significances) will be stored in PostScript file " << fPsFilename -> GetName() << endl; } void MFindSupercutsONOFF::SetPostScriptFile2 (TPostScript &PsFile) { fPsFilename2 = new TPostScript (PsFile); *fLog << "MFindSupercutsONOFF : Results (alpha distributions with excess and significances) will be stored in PostScript file " << fPsFilename2 -> GetName() << endl; } void MFindSupercutsONOFF::SetPsFilenameString (const TString filename) { fPsFilenameString = filename; } void MFindSupercutsONOFF::SetSkipOptimization(Bool_t b) { fSkipOptimization = b; if (fSkipOptimization) { *fLog << "MFindSupercutsONOFF :: SetSkipOptimization " << endl << "Variable fSkipOptimization is kTRUE, and therefore " << "the optimization of supercuts is skipped. Hope that's " << "what you want... " << endl; } } void MFindSupercutsONOFF::SetUseInitialSCParams(Bool_t b) { fUseInitialSCParams = b; if (fUseInitialSCParams) { *fLog << "MFindSupercutsONOFF :: SetUseInitialSCParams " << endl << "Variable fUseInitialSCParams is kTRUE. " << endl; if (fSkipOptimization) { *fLog << "The Initial SC Parameters will be applied on the selected data." << endl; } else { *fLog << "However, fSkipOptimization = kFALSE, and therefore, the " << "the supercuts will be optimized. The final cuts that " << "will be applied to the data will NOT be the initial SC parameters." << endl; } } } Bool_t MFindSupercutsONOFF::SetAlphaSig(Double_t alphasig) { // check that alpha is within the limits 0-90 if (alphasig <= 0 || alphasig > 90) { *fLog << "MFindSupercutsONOFF ::SetAlphaSig; " << "value " << alphasig << " is not within the the " << "logical limits of alpha; 0-90" << endl; return kFALSE; } fAlphaSig = alphasig; return kTRUE; } Bool_t MFindSupercutsONOFF::SetAlphaBkgMin(Double_t alpha) { // check that alpha is within the limits 0-90 if (alpha <= 0 || alpha >= 90) { *fLog << "MFindSupercutsONOFF ::SetAlphaBkgMin; " << "value " << alpha << " is not within the the " << "logical limits of alpha; 0-90" << endl; return kFALSE; } fAlphaBkgMin = alpha; return kTRUE; } Bool_t MFindSupercutsONOFF::SetAlphaBkgMax(Double_t alpha) { // check that alpha is within the limits 0-90 if (alpha <= 0 || alpha > 90.001) { *fLog << "MFindSupercutsONOFF ::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 MFindSupercutsONOFF::CheckAlphaSigBkg() { if (fAlphaSig > fAlphaBkgMin) { *fLog << "MFindSupercutsONOFF ::CheckAlphaSigBkg(); " << endl << "fAlphaSig > fAlphaBkgMin, which should not occur..." << endl << "fAlphaSig = " << fAlphaSig << ", fAlphaBkgMin = " << fAlphaBkgMin << endl; return kFALSE; } if (fAlphaBkgMax < fAlphaBkgMin) { *fLog << "MFindSupercutsONOFF ::CheckAlphaSigBkg(); " << endl << "fAlphaBkgMin > fAlphaBkgMax, which should not occur..." << endl << "fAlphaBkgMin = " << fAlphaBkgMin << ", fAlphaBkgMax = " << fAlphaBkgMax << endl; return kFALSE; } return kTRUE; } /* // Function that computes the normalization factor using COUNTED events ON and OFF // in alpha region defined by fAlphaBkgMin and fAlphaBkgMax Double_t MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg(TH1 *histON, TH1 *histOFF) { Double_t NormFactor = 0.0; Double_t ONEvents = 0.0; Double_t OFFEvents = 0.0; const Double_t SmallQuantity = 0.01; Double_t xlo = 0.0; Double_t xup = 0.0; Double_t width = 0.0; Int_t BinCounterOFF = 0; Int_t BinCounterON = 0; // I make a copy of the histograms so that nothing happens to the // histograms used in argument TH1* HistON = (TH1*) histON->Clone(); TH1* HistOFF = (TH1*) histOFF->Clone(); if ( !HistON ) { *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl << "Clone of ON histogram could not be generated" << endl; return 0.0; } if ( !HistOFF ) { *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl << " Clone of OFF histogram could not be generated" << endl; return 0.0; } // Calculate the number of OFF events in the Bkg region // defined by fAlphaBkgMin and fAlphaBkgMax // ___________________________________________________ Int_t nbinsOFF = HistOFF -> GetNbinsX(); Double_t binwidthOFF = HistOFF -> GetBinWidth(1); for (Int_t i=1; i<=nbinsOFF; i++) { xlo = HistOFF->GetBinLowEdge(i); xup = HistOFF->GetBinLowEdge(i+1); // bin must be completely contained in the bkg region if ( xlo >= (fAlphaBkgMin-SmallQuantity) && xup <= (fAlphaBkgMax+SmallQuantity) ) { width = fabs(xup-xlo); if (fabs(width-binwidthOFF) > SmallQuantity) { *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl << " HistOFF has variable binning, which is not allowed" << endl; return 0.0; } BinCounterOFF++; OFFEvents += HistOFF->GetBinContent(i); } } // Calculate the number of ON events in the Bkg region // defined by fAlphaBkgMin and fAlphaBkgMax // ___________________________________________________ Int_t nbinsON = HistON -> GetNbinsX(); Double_t binwidthON = HistON -> GetBinWidth(1); for (Int_t i=1; i<=nbinsON; i++) { xlo = HistON->GetBinLowEdge(i); xup = HistON->GetBinLowEdge(i+1); // bin must be completely contained in the bkg region if ( xlo >= (fAlphaBkgMin-SmallQuantity) && xup <= (fAlphaBkgMax+SmallQuantity) ) { width = fabs(xup-xlo); if (fabs(width-binwidthON) > SmallQuantity) { *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl << " HistON has variable binning, which is not allowed" << endl; return 0.0; } BinCounterON++; ONEvents += HistON->GetBinContent(i); } } // NormFactor is computed if (ONEvents < SmallQuantity || OFFEvents < SmallQuantity) { *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl << "ONEvents or OFFEvents computed in bkg region are < " << SmallQuantity << endl; return 0.0; } NormFactor = ONEvents/OFFEvents; Double_t error = 1/ONEvents + 1/OFFEvents; error = TMath::Sqrt(error); error = error * NormFactor; // tmp info cout << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg;" << endl << "ON Events in bkg region = " << ONEvents << " (" << BinCounterON << " bins)" << " , OFF Events in bkg region = " << OFFEvents << " (" << BinCounterOFF << " bins)" << endl <<"NormFactor computed from bkg region = " << NormFactor << " +/- " << error << endl; // end temp return NormFactor; } */ // Function that set the values of fThetaMin and fThetaMax and also // fThetaRangeString (with value in miliradians); // data members that are used basically to print/plot // information. Bool_t MFindSupercutsONOFF::SetThetaRange(Double_t ThetaMin, Double_t ThetaMax) { // Check that values are reasonable... well ... i guess this was done // in previous functions... fThetaMin = int(ThetaMin*1000.0); fThetaMax = int(ThetaMax*1000.0); fThetaRangeString = ("ThetaRange"); fThetaRangeString += (fThetaMin); fThetaRangeString += ("_"); fThetaRangeString += (fThetaMax); fThetaRangeString += ("mRad"); return kTRUE; } // Function that sets Size range Bool_t MFindSupercutsONOFF::SetSizeRange(Double_t SizeMin, Double_t SizeMax) { fSizeCutLow = SizeMin; fSizeCutUp = SizeMax; *fLog << "MFindSupercutsONOFF::SetSizeRange" << endl << "Data matrices will be filled with events whose MHillas.fSize " << endl << "is in the range " << fSizeCutLow <<"-"<Print("SizeCols"); Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows(); if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MFindSupercutsONOFF::DefineTrainMatrix; no.of generated events (" << howmanygenerated << ") is incompatible with the no.of requested events (" << howmanytrain << ")" << endl; } *fLog << "training matrix was filled" << endl; *fLog << "=============================================" << endl; //-------------------------- // write out training matrix if (filetrain != "") { TFile filetr(filetrain, "RECREATE"); fMatrixTrain->Write(); filetr.Close(); *fLog << "MFindSupercuts::DefineTrainMatrix; Training matrix was written onto file '" << filetrain << "'" << endl; } return kTRUE; } // -------------------------------------------------------------------------- // -------------------------------------------------------------------------- // // Define the matrix for the test sample // // alltogether 'howmanytest' events are read from file 'nametest' // // the events are selected according to a target distribution 'hreftest' // // Bool_t MFindSupercutsONOFF::DefineTestMatrix(const TString &nametest, MH3 &hreftest, const Int_t howmanytest, const TString &filetest) { if (nametest.IsNull() || howmanytest<=0) return kFALSE; *fLog << "=============================================" << endl; *fLog << "fill test matrix from file '" << nametest << "', select " << howmanytest << " events " << endl; if (!fUseOrigDistribution) { *fLog << " according to a distribution given by the MH3 object '" << hreftest.GetName() << "'" << endl; } else { *fLog << " randomly" << endl; } MParList plist; MTaskList tlist; MReadMarsFile read("Events", nametest); read.DisableAutoScheme(); MFEventSelector2 seltest(hreftest); seltest.SetNumMax(howmanytest); seltest.SetName("selectTest"); if (fUseOrigDistribution) { *fLog << "MFindSupercutsONFF; MFEventSelector2::SetUseOrigDistribution(Bool)" << " DOES NOT EXIST NOW..." << endl; // seltest.SetUseOrigDistribution(kTRUE); } MFillH filltest(fMatrixTest); filltest.SetFilter(&seltest); filltest.SetName("fillMatrixTest"); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTest); //****************************** // entries in MTaskList tlist.AddToList(&read); tlist.AddToList(&seltest); tlist.AddToList(&filltest); //****************************** MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTest"); evtloop.SetProgressBar(&bar); if (!evtloop.Eventloop()) return kFALSE; tlist.PrintStatistics(0, kTRUE); fMatrixTest->Print("SizeCols"); const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows(); if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MFindSupercutsONOFF::DefineTestMatrix; no.of generated events (" << howmanygenerated << ") is incompatible with the no.of requested events (" << howmanytest << ")" << endl; } *fLog << "test matrix was filled" << endl; *fLog << "=============================================" << endl; //-------------------------- // write out test matrix if (filetest != "") { TFile filete(filetest, "RECREATE", ""); fMatrixTest->Write(); filete.Close(); *fLog << "MFindSupercutsONOFF::DefineTestMatrix; Test matrix was written onto file '" << filetest << "'" << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // Define the matrices for the training and test sample respectively // // // Bool_t MFindSupercutsONOFF::DefineTrainTestMatrix( const TString &name, MH3 &href, const Int_t howmanytrain, const Int_t howmanytest, const TString &filetrain, const TString &filetest) { *fLog << "=============================================" << endl; *fLog << "fill training and test matrix from file '" << name << "', select " << howmanytrain << " training and " << howmanytest << " test events " << endl; if (!fUseOrigDistribution) { *fLog << " according to a distribution given by the MH3 object '" << href.GetName() << "'" << endl; } else { *fLog << " randomly" << endl; } MParList plist; MTaskList tlist; MReadMarsFile read("Events", name); read.DisableAutoScheme(); MFEventSelector2 selector(href); selector.SetNumMax(howmanytrain+howmanytest); selector.SetName("selectTrainTest"); selector.SetInverted(); if (fUseOrigDistribution) { *fLog << "MFindSupercutsONFF; MFEventSelector2::SetUseOrigDistribution(Bool)" << " DOES NOT EXIST NOW..." << endl; // selector.SetUseOrigDistribution(kTRUE); } MContinue cont(&selector); cont.SetName("ContTrainTest"); Double_t prob = ( (Double_t) howmanytrain ) / ( (Double_t)(howmanytrain+howmanytest) ); MFRandomSplit split(prob); MFillH filltrain(fMatrixTrain); filltrain.SetFilter(&split); filltrain.SetName("fillMatrixTrain"); // consider this event as candidate for a test event // only if event was not accepted as a training event MContinue conttrain(&split); conttrain.SetName("ContTrain"); MFillH filltest(fMatrixTest); filltest.SetName("fillMatrixTest"); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTrain); plist.AddToList(fMatrixTest); //****************************** //****************************** // entries in MTaskList tlist.AddToList(&read); tlist.AddToList(&cont); tlist.AddToList(&split); tlist.AddToList(&filltrain); tlist.AddToList(&conttrain); tlist.AddToList(&filltest); //****************************** MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTrainTest"); evtloop.SetProgressBar(&bar); Int_t maxev = -1; if (!evtloop.Eventloop(maxev)) return kFALSE; tlist.PrintStatistics(0, kTRUE); fMatrixTrain->Print("SizeCols"); const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows(); if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events (" << generatedtrain << ") is incompatible with the no.of requested events (" << howmanytrain << ")" << endl; } fMatrixTest->Print("SizeCols"); const Int_t generatedtest = fMatrixTest->GetM().GetNrows(); if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events (" << generatedtest << ") is incompatible with the no.of requested events (" << howmanytest << ")" << endl; } *fLog << "training and test matrix were filled" << endl; *fLog << "=============================================" << endl; //-------------------------- // write out training matrix if (filetrain != "") { TFile filetr(filetrain, "RECREATE"); fMatrixTrain->Write(); filetr.Close(); *fLog << "MFindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '" << filetrain << "'" << endl; } //-------------------------- // write out test matrix if (filetest != "") { TFile filete(filetest, "RECREATE", ""); fMatrixTest->Write(); filete.Close(); *fLog << "MFindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '" << filetest << "'" << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // Define the matrices for the training and test OFF sample respectively // // // Bool_t MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange( const TString &name, const Double_t whichfractiontrain, const Double_t whichfractiontest, Double_t ThetaMin, Double_t ThetaMax, const TString &filetrain, const TString &filetest) { *fLog << "=============================================" << endl; *fLog << "Fill training and testing OFF matrices from file '" << name << "', select a fraction of " << whichfractiontrain << " events for the training and a fraction of " << endl << whichfractiontest << " for the testing" << endl; MParList plist; MTaskList tlist; MReadMarsFile read("Events", name); read.DisableAutoScheme(); // Cuts in Theta // TString ThetaCutMinString ("ThetaOrig.fVal>"); TString ThetaCutMinString ("MPointingPos.fZd>"); // new! // for magic ThetaCutMinString += ThetaMin; MContinue ThetaCutMin(ThetaCutMinString); ThetaCutMin.SetInverted(); //TString ThetaCutMaxString ("ThetaOrig.fVal<"); TString ThetaCutMaxString ("MPointingPos.fZd<"); // new! // for magic ThetaCutMaxString += ThetaMax; MContinue ThetaCutMax(ThetaCutMaxString); ThetaCutMax.SetInverted(); // TMP // Cuts in Size, TString SizeCutMinString ("MHillas.fSize>"); SizeCutMinString += fSizeCutLow; MContinue SizeCutMin(SizeCutMinString); SizeCutMin.SetInverted(); TString SizeCutMaxString ("MHillas.fSize<"); SizeCutMaxString += fSizeCutUp; MContinue SizeCutMax(SizeCutMaxString); SizeCutMax.SetInverted(); // ENDTMP Double_t prob = whichfractiontrain /(whichfractiontrain+whichfractiontest); MFRandomSplit split(prob); MFillH filltrain(fMatrixTrainOFF); filltrain.SetName("fillMatrixTrainOFF"); filltrain.SetFilter(&split); // consider this event as candidate for a test event // only if event was not accepted as a training event MContinue conttrain(&split); conttrain.SetName("ContTrain"); MFillH filltest(fMatrixTestOFF); filltest.SetName("fillMatrixTestOFF"); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTrainOFF); plist.AddToList(fMatrixTestOFF); //****************************** // entries in MTaskList tlist.AddToList(&read); // tlist.AddToList(&ThetaCutMin); // tlist.AddToList(&ThetaCutMax); //TMP tlist.AddToList(&SizeCutMin); tlist.AddToList(&SizeCutMax); //ENDTMP tlist.AddToList(&split); tlist.AddToList(&filltrain); tlist.AddToList(&conttrain); tlist.AddToList(&filltest); //****************************** MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTrainTestOFF"); evtloop.SetProgressBar(&bar); Int_t maxev = -1; if (!evtloop.Eventloop(maxev)) return kFALSE; tlist.PrintStatistics(0, kTRUE); fMatrixTrainOFF->Print("SizeCols"); fMatrixTestOFF->Print("SizeCols"); *fLog << "train and test matrix OFF were filled" << endl; *fLog << "=============================================" << endl; //-------------------------- // write out training matrix if (filetrain != "") { TFile filetr(filetrain, "RECREATE"); fMatrixTrainOFF->Write(); filetr.Close(); *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange; Training matrix was written onto file '" << filetrain << "'" << endl; } //-------------------------- // write out test matrix if (filetest != "") { TFile filete(filetest, "RECREATE", ""); fMatrixTestOFF->Write(); filete.Close(); *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange; Test matrix was written onto file '" << filetest << "'" << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // Define the matrices for the training and test sample respectively // // // Bool_t MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange( const TString &name, const Double_t whichfractiontrain, const Double_t whichfractiontest, Double_t ThetaMin, Double_t ThetaMax, const TString &filetrain, const TString &filetest) { *fLog << "=============================================" << endl; *fLog << "Fill training and testing ON matrices from file '" << name << "', select a fraction of " << whichfractiontrain << " events for the training and a fraction of " << endl << whichfractiontest << " for the testing" << endl; MParList plist; MTaskList tlist; MReadMarsFile read("Events", name); read.DisableAutoScheme(); // Cuts in Theta //TString ThetaCutMinString ("ThetaOrig.fVal>"); TString ThetaCutMinString ("MPointingPos.fZd>"); // for magic ThetaCutMinString += ThetaMin; MContinue ThetaCutMin(ThetaCutMinString); ThetaCutMin.SetInverted(); //TString ThetaCutMaxString ("ThetaOrig.fVal<"); TString ThetaCutMaxString ("MPointingPos.fZd<"); // for magic ThetaCutMaxString += ThetaMax; MContinue ThetaCutMax(ThetaCutMaxString); ThetaCutMax.SetInverted(); // TMP // Cuts in Size, TString SizeCutMinString ("MHillas.fSize>"); SizeCutMinString += fSizeCutLow; MContinue SizeCutMin(SizeCutMinString); SizeCutMin.SetInverted(); TString SizeCutMaxString ("MHillas.fSize<"); SizeCutMaxString += fSizeCutUp; MContinue SizeCutMax(SizeCutMaxString); SizeCutMax.SetInverted(); // ENDTMP Double_t prob = whichfractiontrain/(whichfractiontrain + whichfractiontest); MFRandomSplit split(prob); MFillH filltrain(fMatrixTrain); filltrain.SetName("fillMatrixTrain"); filltrain.SetFilter(&split); // consider this event as candidate for a test event // only if event was not accepted as a training event MContinue conttrain(&split); conttrain.SetName("ContTrain"); MFillH filltest(fMatrixTest); filltest.SetName("fillMatrixTest"); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTrain); plist.AddToList(fMatrixTest); //****************************** // entries in MTaskList tlist.AddToList(&read); // tlist.AddToList(&ThetaCutMin); // tlist.AddToList(&ThetaCutMax); //TMP tlist.AddToList(&SizeCutMin); tlist.AddToList(&SizeCutMax); //ENDTMP tlist.AddToList(&split); tlist.AddToList(&filltrain); tlist.AddToList(&conttrain); tlist.AddToList(&filltest); //****************************** MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTrainTest"); evtloop.SetProgressBar(&bar); Int_t maxev = -1; if (!evtloop.Eventloop(maxev)) return kFALSE; tlist.PrintStatistics(0, kTRUE); fMatrixTrain->Print("SizeCols"); fMatrixTest->Print("SizeCols"); *fLog << "training and test matrix were filled" << endl; *fLog << "=============================================" << endl; //-------------------------- // write out training matrix if (filetrain != "") { TFile filetr(filetrain, "RECREATE"); fMatrixTrain->Write(); filetr.Close(); *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange; Train matrix was written onto file '" << filetrain << "'" << endl; } //-------------------------- // write out test matrix if (filetest != "") { TFile filete(filetest, "RECREATE", ""); fMatrixTest->Write(); filete.Close(); *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange; Test matrix was written onto file '" << filetest << "'" << endl; } return kTRUE; } /// **********/// // -------------------------------------------------------------------------- // // Read training and test matrices from files // // Bool_t MFindSupercutsONOFF::ReadMatrix(const TString &filetrain, const TString &filetest) { //-------------------------- // read in training matrix TFile filetr(filetrain); fMatrixTrain->Read("MatrixTrain"); fMatrixTrain->Print("SizeCols"); *fLog << "MFindSupercuts::ReadMatrix; Training matrix was read in from file '" << filetrain << "'" << endl; filetr.Close(); //-------------------------- // read in test matrix TFile filete(filetest); fMatrixTest->Read("MatrixTest"); fMatrixTest->Print("SizeCols"); *fLog << "MFindSupercuts::ReadMatrix; Test matrix was read in from file '" << filetest << "'" << endl; filete.Close(); return kTRUE; } // Read training and test matrices OFF from files // // Bool_t MFindSupercutsONOFF::ReadMatrixOFF(const TString &filetrain, const TString &filetest) { //-------------------------- // read in training matrix TFile filetr(filetrain); fMatrixTrainOFF->Read("MatrixTrainOFF"); fMatrixTrainOFF->Print("SizeCols"); *fLog << "MFindSupercutsONOFF::ReadMatrixOFF; Training matrix OFF was read in from file '" << filetrain << "'" << endl; filetr.Close(); //-------------------------- // read in test matrix TFile filete(filetest); fMatrixTestOFF->Read("MatrixTestOFF"); fMatrixTestOFF->Print("SizeCols"); *fLog << "MFindSupercutsONOFF::ReadMatrixONOFF; Test matrix OFF was read in from file '" << filetest << "'" << endl; filete.Close(); return kTRUE; } // Function to compute the normalization factor for Train sample. // The normalization factor is defined as the ratio of OFF/ON events. Bool_t MFindSupercutsONOFF::ComputeNormFactorTrain() { Int_t EventsInTrainMatrixON = fMatrixTrain->GetM().GetNrows(); Int_t EventsInTrainMatrixOFF = fMatrixTrainOFF->GetM().GetNrows(); // Info... *fLog << "MFindSupercutsONOFF::ComputeNormFactorTrain; " << endl << "EventsInTrainMatrixON = " << EventsInTrainMatrixON << " , EventsInTrainMatrixOFF = " << EventsInTrainMatrixOFF << endl; fNormFactorTrain = double(EventsInTrainMatrixON)/double(EventsInTrainMatrixOFF); *fLog << "MFindSupercutsONOFF::ComputeNormFactorTrain; fNormFactorTrain is : " << fNormFactorTrain << endl; return kTRUE; } // Function to compute the normalization factor for Test sample. // The normalization factor is defined as the ratio of OFF/ON events. Bool_t MFindSupercutsONOFF::ComputeNormFactorTest() { Int_t EventsInTestMatrixON = fMatrixTest->GetM().GetNrows(); Int_t EventsInTestMatrixOFF = fMatrixTestOFF->GetM().GetNrows(); // Info... *fLog << "MFindSupercutsONOFF::ComputeNormFactorTest; " << endl << "EventsInTestMatrixON = " << EventsInTestMatrixON << " , EventsInTestMatrixOFF = " << EventsInTestMatrixOFF << endl; fNormFactorTest = double(EventsInTestMatrixON)/double(EventsInTestMatrixOFF); *fLog << "MFindSupercutsONOFF::ComputeNormFactorTest; fNormFactorTest is : " << fNormFactorTest << endl; return kTRUE; } Bool_t MFindSupercutsONOFF::SetGammaEfficiency(Double_t gammaeff) { if (gammaeff < 0.0 || gammaeff >1.0) { *fLog << "MFindSupercutsONOFF::SetGammaEfficiency; " << endl << "The argument of SetGammaEfficiency (" << gammaeff << ") is outside the range 0-1" << endl << "This is not allowed..." << endl; return kFALSE; } fGammaEfficiency = gammaeff; // TEMP INFO cout << "fGammaEfficiency has been set to : " << fGammaEfficiency << endl; // END TEMP return kTRUE; } //------------------------------------------------------------------------ // // Steering program for optimizing the supercuts // --------------------------------------------- // // the criterion for the 'optimum' is // // - a maximum significance of the gamma signal in the alpha plot, // in which the supercuts have been applied // // The various steps are : // // - setup the 2 event loops (ON and OFF) to be executed for each call to fcnSupercuts // - call TMinuit to do the minimization of (-significance) : // the fcnSupercuts function calculates the significance // for the current cut values // for this - the 2 alpha plots (ON and OFF) are produced in the 2 event loops, // in which the SAME cuts have been applied // - the significance of the gamma signal is extracted from the // the 2 alpha plots (ON OFF) using MHFindSignificanceONOFF::FindSigmaONOFF. // // // Needed as input : (to be set by the Set functions) // // - fFilenameParam name of file to which optimum values of the // parameters are written // // - fHadronnessName name of container where MSupercutsCalcONOFF will // put the supercuts hadronness for ON data // - fHadronnessNameOFF name of container where MSupercutsCalcONOFF will // put the supercuts hadronness for OFF data // // - fAlphaSig, fAlphaBkgMin, fAlphaBkgMax ; Signal and Bkg region in alpha histo. // - for the minimization, the starting values of the parameters are taken // - from file parSCinit (if it is != "") // - or from the arrays params and/or steps // //---------------------------------------------------------------------- Bool_t MFindSupercutsONOFF::FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps) { // Setup the event loops which will be executed in the // fcnSupercuts function of MINUIT // // parSCinit is the name of the file containing the initial values // of the parameters; // if parSCinit = "" 'params' and 'steps' are taken as initial values *fLog << "=============================================" << endl; *fLog << "Setup event loop for fcnSupercuts" << endl; if (fHadronnessName.IsNull()) { *fLog << "MFindSupercutsONOFF::FindParams; hadronness name for ON data is not defined... aborting" << endl; return kFALSE; } if (fHadronnessNameOFF.IsNull()) { *fLog << "MFindSupercutsONOFF::FindParams; hadronness name for OFF data is not defined... aborting" << endl; return kFALSE; } if (fMatrixTrain == NULL) { *fLog << "MFindSupercuts::FindParams; training matrix is not defined... aborting" << endl; return kFALSE; } if (fMatrixTrain->GetM().GetNrows() <= 0) { *fLog << "MFindSupercuts::FindParams; training matrix has no entries" << endl; return kFALSE; } if (fMatrixTrainOFF == NULL) { *fLog << "MFindSupercutsONOFF::FindParams; training matrix OFF is not defined... aborting" << endl; return kFALSE; } if (fMatrixTrainOFF->GetM().GetNrows() <= 0) { *fLog << "MFindSupercutsONOFF::FindParams; training matrix OFF has no entries" << endl; return kFALSE; } if (fAlphaDistributionsRootFilename.IsNull()) { *fLog << "MFindSupercutsONOFF::FindParams; fAlphaDistributionsRootFilename is not defined... program aborting..." << endl; return kFALSE; } //-------------------------------- // create container for the supercut parameters // and set them to their initial values MSupercuts super; // take initial values from file parSCinit if (parSCinit != "") { TFile inparam(parSCinit); super.Read("MSupercuts"); inparam.Close(); *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from file " << parSCinit << endl; } // take initial values from 'params' and/or 'steps' else if (params.GetSize() != 0 || steps.GetSize() != 0 ) { if (params.GetSize() != 0) { *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from 'params'" << endl; super.SetParameters(params); } if (steps.GetSize() != 0) { *fLog << "MFindSupercutsONOFF::FindParams; initial step sizes are taken from 'steps'" << endl; super.SetStepsizes(steps); } /* // TMP // Print parameters if (params.GetSize() == steps.GetSize()) { *fLog << "MFindSupercutsONOFF; SC parameters and Setps are: " << endl; for (Int_t z = 0; z < params.GetSize(); z++) { cout << params[z] << setw(20) << steps[z] << endl; } } // ENDTMP */ /* // TMP2 TArrayD paramsBis = super.GetParameters(); TArrayD stepsBis = super.GetStepsizes(); if (paramsBis.GetSize() == stepsBis.GetSize()) { *fLog << "MFindSupercutsONOFF; SC parametersBis and SetpsBis are: " << endl; for (Int_t z = 0; z < paramsBis.GetSize(); z++) { cout << paramsBis[z] << setw(20) << stepsBis[z] << endl; } } // ENDTMP2 */ else { *fLog << "MFindSupercutsONOFF; ERROR: params.GetSize() != steps.GetSize() " << endl; } } else { *fLog << "MFindSupercutsONOFF::FindParams; initial values and step sizes are taken from the MSupercuts constructor" << endl; } // Computation of fNormFactorTrain and creation of a container for // storing this value and include it in the MParList for being // used by function fcnSupercuts if (!ComputeNormFactorTrain()) { *fLog << "Normalization factor for train sample (ON-OFF) could not be computed. Aborting..." << endl; return kFALSE; } *fLog << "MFindSupercutsONOFF::FindParams; fNormFactorTrain = " << fNormFactorTrain << endl; MDataValue NormFactorCont(fNormFactorTrain); NormFactorCont.SetName("NormFactorTrain"); // Value of fAlphaSig, fAlphaBkgMin and fAlphaBkgMax // are stored in MDataValue containers, // then they will be passed to the MParList and will be // accessible to fcnsupercuts MDataValue AlphaSigCont(fAlphaSig); AlphaSigCont.SetName("AlphaSigValue"); MDataValue AlphaBkgMinCont(fAlphaBkgMin); AlphaBkgMinCont.SetName("AlphaBkgMinValue"); MDataValue AlphaBkgMaxCont(fAlphaBkgMax); AlphaBkgMaxCont.SetName("AlphaBkgMaxValue"); // Value of fUseFittedQuantities is stored in container // and passed to the MParList and will be // accessible to fcnsupercuts MDataValue UseFittedQuantitiesCont(fUseFittedQuantities); UseFittedQuantitiesCont.SetName("UseFittedQuantitiesValue"); // Value of fNormFactorFromAlphaBkg is stored in container // and passed to the MParList and will be // accessible to fcnsupercuts MDataValue UseNormFactorFromAlphaBkgCont(fNormFactorFromAlphaBkg); UseNormFactorFromAlphaBkgCont.SetName("UseNormFactorFromAlphaBkgValue"); // A vector of pointers to objects of the class MEvtLoop is defined. // One of the pointed loops will be used to compute ON data alpha plot // and the other for computing OFF data alpha plot. MEvtLoop evtloopfcn[2] = {MEvtLoop("ONDataEvtLoopFCN"), MEvtLoop("OFFDataEvtLoopFCN")}; // ****************************************************************** // EVENT LOOP FOR COMPUTING ALPHA HISTOGRAM FOR ON DATA IS SET // ****************************************************************** // ----------------------------------------------------------------- MParList parlistfcn; MTaskList tasklistfcn; // loop over rows of matrix MMatrixLoop loop(fMatrixTrain); //-------------------------------- // calculate supercuts hadronness fCalcHadTrain->SetHadronnessName(fHadronnessName); // Boolean variable that controls in class MSupercutsCalcONOFF // wether the supercuts are stored or not // is set to the default value, kFALSE. (no storage of supercuts) fCalcHadTrain -> SetStoreAppliedSupercuts(kFALSE); // Set boolean variable that controls wether cuts are // dynamic or static fCalcHadTrain -> SetVariableUseStaticCuts(fUseStaticCuts); if (!fUseStaticCuts) { // Set boolean variable that controls wether the theta variable // is used or not in the computation of the dynamical cuts fCalcHadTrain -> SetVariableNotUseTheta(fNotUseTheta); } // apply the supercuts MF scfilter(fHadronnessName+".fHadronness>0.5"); MContinue supercuts(&scfilter); // plot |alpha| const TString mh3Name = "AlphaFcn"; MBinning binsalpha("Binning"+mh3Name); //binsalpha.SetEdges(54, -12.0, 96.0); binsalpha.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl; // |alpha| is assumed to be in column 7 of the matrix MH3 alpha("MatrixTrain[7]"); alpha.SetName(mh3Name); MFillH fillalpha(&alpha); // book histograms to be filled during the optimization // ! not in the event loop ! MHCT1Supercuts plotsuper; plotsuper.SetupFill(&parlistfcn); //****************************** // entries in MParList (extension of old MParList) parlistfcn.AddToList(&tasklistfcn); parlistfcn.AddToList(&super); parlistfcn.AddToList(&NormFactorCont); parlistfcn.AddToList(&AlphaSigCont); parlistfcn.AddToList(&AlphaBkgMinCont); parlistfcn.AddToList(&AlphaBkgMaxCont); parlistfcn.AddToList(&UseFittedQuantitiesCont); parlistfcn.AddToList(&UseNormFactorFromAlphaBkgCont); parlistfcn.AddToList(fCam); parlistfcn.AddToList(fMatrixTrain); parlistfcn.AddToList(&binsalpha); parlistfcn.AddToList(&alpha); parlistfcn.AddToList(&plotsuper); //****************************** // entries in MTaskList tasklistfcn.AddToList(&loop); tasklistfcn.AddToList(fCalcHadTrain); tasklistfcn.AddToList(&supercuts); tasklistfcn.AddToList(&fillalpha); //****************************** // &evtloopfcn[0] = new MEvtLoop("ONDataEvtLoopFCN"); evtloopfcn[0].SetParList(&parlistfcn); // MEvtLoop evtloopfcn("EvtLoopFCN"); // evtloopfcn.SetParList(&parlistfcn); *fLog << "Event loop for computing ON data alpha histo in fcnSupercuts has been setup" << endl; //---- End of setting of loop for compuing ON data alpha plot --------------- // ****************************************************************** // EVENT LOOP FOR COMPUTING ALPHA HISTOGRAM FOR OFF DATA IS SET // ****************************************************************** // ----------------------------------------------------------------- MParList parlistfcnOFF; MTaskList tasklistfcnOFF; // loop over rows of matrix MMatrixLoop loopOFF(fMatrixTrainOFF); //-------------------------------- // calculate supercuts hadronness fCalcHadTrainOFF->SetHadronnessName(fHadronnessNameOFF); // Boolean variable that controls in class MSupercutsCalcONOFF // wether the supercuts are stored or not // is set to the default value, kFALSE. (no storage of supercuts) fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kFALSE); // Set boolean variable that controls wether cuts are // dynamic or static fCalcHadTrainOFF -> SetVariableUseStaticCuts(fUseStaticCuts); if (!fUseStaticCuts) { // Set boolean variable that controls wether the theta variable // is used or not in the computation of the dynamica cuts fCalcHadTrainOFF -> SetVariableNotUseTheta(fNotUseTheta); } // apply the supercuts MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5"); MContinue supercutsOFF(&scfilterOFF); // plot |alpha| const TString mh3NameOFF = "AlphaOFFFcn"; MBinning binsalphaOFF("Binning"+mh3NameOFF); //binsalphaOFF.SetEdges(54, -12.0, 96.0); binsalphaOFF.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl; // |alpha| is assumed to be in column 7 of the matrix MH3 alphaOFF("MatrixTrainOFF[7]"); alphaOFF.SetName(mh3NameOFF); MFillH fillalphaOFF(&alphaOFF); //****************************** // entries in MParList (extension of old MParList) parlistfcnOFF.AddToList(&tasklistfcnOFF); parlistfcnOFF.AddToList(&super); parlistfcnOFF.AddToList(fCam); parlistfcnOFF.AddToList(fMatrixTrainOFF); parlistfcnOFF.AddToList(&binsalphaOFF); parlistfcnOFF.AddToList(&alphaOFF); //****************************** // entries in MTaskList tasklistfcnOFF.AddToList(&loopOFF); tasklistfcnOFF.AddToList(fCalcHadTrainOFF); tasklistfcnOFF.AddToList(&supercutsOFF); tasklistfcnOFF.AddToList(&fillalphaOFF); //****************************** // &evtloopfcn[1] = new MEvtLoop("OFFDataEvtLoopFCN"); evtloopfcn[1].SetParList(&parlistfcnOFF); // MEvtLoop evtloopfcn("EvtLoopFCN"); // evtloopfcn.SetParList(&parlistfcn); *fLog << "Event loop for computing OFF data alpha histo in fcnSupercuts has been setup" << endl; //---- End of setting of loop for compuing OFF data alpha plot --------------- // address of evtloopfcn is used in CallMinuit() //----------------------------------------------------------------------- // //---------- Start of minimization part -------------------- // // Do the minimization with MINUIT // // Be careful: This is not thread safe // *fLog << "========================================================" << endl; *fLog << "Start minimization for supercuts" << endl; // ------------------------------------------- // prepare call to MINUIT // // get initial values of parameters fVinit = super.GetParameters(); fStep = super.GetStepsizes(); TString name[fVinit.GetSize()]; fStep.Set(fVinit.GetSize()); fLimlo.Set(fVinit.GetSize()); fLimup.Set(fVinit.GetSize()); fFix.Set(fVinit.GetSize()); fNpar = fVinit.GetSize(); for (UInt_t i=0; i= 48) // { // fStep[i] = 0.0; // fFix[i] = 1; // } //} // ------------------------------------------- // call MINUIT Bool_t rc; if(fSkipOptimization) { *fLog << "Parameter optimization is skipped. Using previously optimized parameters." << endl; rc = kTRUE; if(fUseInitialSCParams) { // write Initial parameter values onto root file filenameParam // so that later they can be applied on the data *fLog << "Initial SC parameter values are written onto file '" << fFilenameParam << "' :" << endl; TFile outparam(fFilenameParam, "RECREATE"); super.Write(); outparam.Close(); const TArrayD &check = super.GetParameters(); for (Int_t i=0; iGetM().GetNrows() <= 0) { *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; test matrix ON has no entries" << endl; return kFALSE; } if (fMatrixTestOFF->GetM().GetNrows() <= 0) { *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; test matrix OFF has no entries" << endl; return kFALSE; } // ------------- TEST the supercuts ------------------------------ // *fLog << "Test the supercuts on the test sample ON and OFF" << endl; // ----------------------------------------------------------------- // read optimum parameter values from file filenameParam // into array 'supercutsPar' TFile inparam(fFilenameParam); MSupercuts scin; scin.Read("MSupercuts"); inparam.Close(); *fLog << "Optimum parameter values for supercuts were read from file '"; *fLog << fFilenameParam << "' :" << endl; const TArrayD &supercutsPar = scin.GetParameters(); for (Int_t i=0; i SetHadronnessName(fHadronnessNameOFF); fCalcHadTest -> SetHadronnessName(fHadronnessName); // Settings related to the storage of teh supercuts applied fCalcHadTest -> SetSupercutsAppliedName(fTestONSupercutsAppliedName); fCalcHadTest -> SetStoreAppliedSupercuts(kTRUE); fCalcHadTestOFF -> SetSupercutsAppliedName(fTestOFFSupercutsAppliedName); fCalcHadTestOFF -> SetStoreAppliedSupercuts(kTRUE); // Set boolean variable that controls wether cuts are // dynamic or static fCalcHadTest -> SetVariableUseStaticCuts(fUseStaticCuts); fCalcHadTestOFF -> SetVariableUseStaticCuts(fUseStaticCuts); if (!fUseStaticCuts) { // Set boolean variable that controls wether the theta variable // is used or not in the computation of the dynamical cuts fCalcHadTest -> SetVariableNotUseTheta(fNotUseTheta); fCalcHadTestOFF -> SetVariableNotUseTheta(fNotUseTheta); } // apply the supercuts to OFF data MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5"); MContinue applysupercutsOFF(&scfilterOFF); // apply the supercuts to ON data MF scfilter(fHadronnessName+".fHadronness>0.5"); MContinue applysupercuts(&scfilter); // **************************************************** // Filling OFF alpha distribution MParList parlistOFF; MTaskList tasklistOFF; MMatrixLoop loopOFF(fMatrixTestOFF); // plot alpha before applying the supercuts const TString mh3NameOFFB = "AlphaOFFBefore"; MBinning binsalphaOFFbef("Binning"+mh3NameOFFB); // binsalphaOFFbef.SetEdges(54, -12.0, 96.0); binsalphaOFFbef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); // |alpha| is assumed to be in column 7 of the matrix MH3 alphaOFFbefore("MatrixTestOFF[7]"); alphaOFFbefore.SetName(mh3NameOFFB); MFillH fillalphaOFFbefore(&alphaOFFbefore); fillalphaOFFbefore.SetName("FillAlphaOFFBefore"); // plot alpha OFF after applying the supercuts const TString mh3NameOFFA = "AlphaOFFAfter"; MBinning binsalphaOFFaft("Binning"+mh3NameOFFA); // binsalphaOFFaft.SetEdges(54, -12.0, 96.0); binsalphaOFFaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); MH3 alphaOFFafter("MatrixTestOFF[7]"); alphaOFFafter.SetName(mh3NameOFFA); MFillH fillalphaOFFafter(&alphaOFFafter); fillalphaOFFafter.SetName("FillAlphaOFFAfter"); //****************************** // entries in MParList parlistOFF.AddToList(&tasklistOFF); parlistOFF.AddToList(&supercuts); parlistOFF.AddToList(&supapptestOFF); parlistOFF.AddToList(fCam); parlistOFF.AddToList(fMatrixTestOFF); parlistOFF.AddToList(&binsalphaOFFbef); parlistOFF.AddToList(&alphaOFFbefore); parlistOFF.AddToList(&binsalphaOFFaft); parlistOFF.AddToList(&alphaOFFafter); //****************************** // Entries in MtaskList tasklistOFF.AddToList(&loopOFF); tasklistOFF.AddToList(&fillalphaOFFbefore); tasklistOFF.AddToList(fCalcHadTestOFF); tasklistOFF.AddToList(&applysupercutsOFF); tasklistOFF.AddToList(&fillalphaOFFafter); //****************************** MProgressBar barOFF; MEvtLoop evtloopOFF; evtloopOFF.SetParList(&parlistOFF); evtloopOFF.SetName("EvtLoopTestParamsOFF"); evtloopOFF.SetProgressBar(&barOFF); Int_t maxeventsOFF = -1; if (!evtloopOFF.Eventloop(maxeventsOFF)) return kFALSE; tasklistOFF.PrintStatistics(0, kTRUE); //------------------------------------------- // draw the alpha plots MH3* alphaOFFBefore = (MH3*)parlistOFF.FindObject(mh3NameOFFB, "MH3"); TH1 &alphaOFFHistb = alphaOFFBefore->GetHist(); TString alphaOFFHistbName ("TestAlphaOFFBeforeCuts"); if (fThetaRangeString.IsNull()) { *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; " << " fThetaRangeString is NOT defined..." << endl; alphaOFFHistbName += ("ThetaRangeStringUndefined"); } else { alphaOFFHistbName += (fThetaRangeString); } alphaOFFHistb.SetXTitle("|\\alpha| [\\circ]"); alphaOFFHistb.SetName(alphaOFFHistbName); MH3* alphaOFFAfter = (MH3*)parlistOFF.FindObject(mh3NameOFFA, "MH3"); TH1 &alphaOFFHista = alphaOFFAfter->GetHist(); TString alphaOFFHistaName ("TestAlphaOFFAfterCuts"); if (fThetaRangeString.IsNull()) { *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; " << " fThetaRangeString is NOT defined..." << endl; alphaOFFHistaName += ("ThetaRangeStringUndefined"); } else { alphaOFFHistaName += (fThetaRangeString); } alphaOFFHista.SetXTitle("|\\alpha| [\\circ]"); alphaOFFHista.SetName(alphaOFFHistaName); // ***************************************************** //------------------------------------------- MParList parlist2; MTaskList tasklist2; MMatrixLoop loop(fMatrixTest); // plot alpha before applying the supercuts const TString mh3NameB = "AlphaBefore"; MBinning binsalphabef("Binning"+mh3NameB); //binsalphabef.SetEdges(54, -12.0, 96.0); binsalphabef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); // |alpha| is assumed to be in column 7 of the matrix MH3 alphabefore("MatrixTest[7]"); alphabefore.SetName(mh3NameB); MFillH fillalphabefore(&alphabefore); fillalphabefore.SetName("FillAlphaBefore"); // plot alpha after applying the supercuts const TString mh3NameA = "AlphaAfter"; MBinning binsalphaaft("Binning"+mh3NameA); // binsalphaaft.SetEdges(54, -12.0, 96.0); binsalphaaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); MH3 alphaafter("MatrixTest[7]"); alphaafter.SetName(mh3NameA); MFillH fillalphaafter(&alphaafter); fillalphaafter.SetName("FillAlphaAfter"); //****************************** // entries in MParList parlist2.AddToList(&tasklist2); parlist2.AddToList(&supercuts); parlist2.AddToList(&supapptestON); parlist2.AddToList(fCam); parlist2.AddToList(fMatrixTest); parlist2.AddToList(&binsalphabef); parlist2.AddToList(&alphabefore); parlist2.AddToList(&binsalphaaft); parlist2.AddToList(&alphaafter); //****************************** // entries in MTaskList tasklist2.AddToList(&loop); tasklist2.AddToList(&fillalphabefore); tasklist2.AddToList(fCalcHadTest); tasklist2.AddToList(&applysupercuts); tasklist2.AddToList(&fillalphaafter); //****************************** MProgressBar bar2; MEvtLoop evtloop2; evtloop2.SetParList(&parlist2); evtloop2.SetName("EvtLoopTestParams"); evtloop2.SetProgressBar(&bar2); Int_t maxevents2 = -1; if (!evtloop2.Eventloop(maxevents2)) return kFALSE; tasklist2.PrintStatistics(0, kTRUE); //------------------------------------------- // draw the alpha plots MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3"); TH1 &alphaHist1 = alphaBefore->GetHist(); TString alphaHist1Name ("TestAlphaBeforeCuts"); if (fThetaRangeString.IsNull()) { *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; " << " fThetaRangeString is NOT defined..." << endl; alphaHist1Name += ("ThetaRangeStringUndefined"); } else { alphaHist1Name += (fThetaRangeString); } alphaHist1.SetName(alphaHist1Name); alphaHist1.SetXTitle("|\\alpha| [\\circ]"); MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3"); TH1 &alphaHist2 = alphaAfter->GetHist(); TString alphaHist2Name ("TestAlphaAfterCuts"); if (fThetaRangeString.IsNull()) { *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; " << " fThetaRangeString is NOT defined..." << endl; alphaHist2Name += ("ThetaRangeStringUndefined"); } else { alphaHist2Name += (fThetaRangeString); } alphaHist2.SetXTitle("|\\alpha| [\\circ]"); alphaHist2.SetName(alphaHist2Name); //fPsFilename -> NewPage(); // Canvas is deleted after saving the plot in ps file // This is to allow for GRID analysis TCanvas *c = new TCanvas("TestAlphaBeforeAfterSC", "TestAlphaBeforeAfterSC", 600, 600); c->Divide(2,2); c->cd(1); alphaOFFHistb.DrawCopy(); c->cd(2); alphaOFFHista.DrawCopy(); c->cd(3); alphaHist1.DrawCopy(); c->cd(4); alphaHist2.DrawCopy(); c->Modified(); c-> Update(); // ******************************************** // TMP solutions while the TPostScript thing is not working // PsFileName for storing these plots is derived // from PsFilenameString. cout << "Alpha distributions will be saved in PostScript file " ; if (!fPsFilenameString.IsNull()) { TString filename = (fPsFilenameString); filename += ("AlphaTESTSampleONOFFBeforeAfterCuts.ps"); cout << filename << endl; c -> SaveAs(filename); } // END OF TEMPORAL SOLUTION // ******************************************** // Canvas is deleted after saving the plot in ps file // This is to allow for GRID analysis delete c; //------------------------------------------- // fit alpha distribution to get the number of excess events and // calculate significance of gamma signal in the alpha plot const Double_t alphasig = fAlphaSig; const Double_t alphamin = fAlphaBkgMin; // alpha min for bkg region in ON data const Double_t alphamax = fAlphaBkgMax; const Int_t degree = 2; const Int_t degreeOFF = 2; const Bool_t drawpoly = kTRUE; const Bool_t fitgauss = kTRUE; const Bool_t print = kTRUE; const Bool_t saveplots = kTRUE; // Save plots in Psfile if (!ComputeNormFactorTest()) { *fLog << "Normalization factor for test sample (ON-OFF) couold not be computed. Aborting..." << endl; return kFALSE; } // Normalization factor factor is corrected using the estimated // number of excess events and the total number of ON events if (fTuneNormFactor) { // Run MHFindSignificanceONOFF::FindSigmaONOFF // to get estimation of number of excess gamma events. // This number will be used in the correction of the // normalization factor // No plots will be produced this time... const Bool_t drawpolyTest = kFALSE; const Bool_t fitgaussTest = kFALSE; const Bool_t printTest = kFALSE; const Bool_t saveplotsTest = kFALSE; // Save plots in Psfile // TPostScript* DummyPs = new TPostScript("dummy.ps"); TString DummyPs = ("Dummy"); MHFindSignificanceONOFF findsigTest; findsigTest.SetRebin(kTRUE); findsigTest.SetReduceDegree(kFALSE); findsigTest.SetUseFittedQuantities(fUseFittedQuantities); if (findsigTest.FindSigmaONOFF(&alphaHist2, &alphaOFFHista, fNormFactorTest, alphamin, alphamax, degree, degreeOFF, alphasig, drawpolyTest, fitgaussTest, printTest, saveplotsTest, DummyPs)) { // Update values of Nex and SigmaLiMa for Test sammple fSigmaLiMaTest = double(findsigTest.GetSignificance()); if(fUseFittedQuantities) { fNexTest = double(findsigTest.GetNexONOFFFitted()); } else { fNexTest = double(findsigTest.GetNexONOFF()); } } else { *fLog << "Normalization Factor tuning for Train sample is not possible" << endl; fSigmaLiMaTest = 0.0; fNexTest = 0.0; } // DummyPs -> Close(); //DummyPs = NULL; // delete DummyPs; Int_t EventsInTestMatrixOFF = fMatrixTestOFF->GetM().GetNrows(); Double_t Ngammas = fNexTest/fGammaEfficiency; Double_t GammaFraction = Ngammas/EventsInTestMatrixOFF; *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; " << "fNormFactorTest is corrected with fraction of gammas in ON " << "Data Test sample and the number of OFF events before cuts." << endl << "EventsInTestMatrixOFF = " << EventsInTestMatrixOFF << " Excess events in Test ON sample = " << fNexTest << " fGammaEfficiency = " << fGammaEfficiency << endl << "fNormFactorTest Value before correction is " << fNormFactorTest << endl; fNormFactorTest = fNormFactorTest - GammaFraction; *fLog << "fNormFactorTest Value after correction is " << fNormFactorTest << endl; } if (fNormFactorFromAlphaBkg) { // Normalization factor computed using alpha bkg region Double_t NewNormFactor = ComputeNormFactorFromAlphaBkg(&alphaHist2, &alphaOFFHista, fAlphaBkgMin, fAlphaBkgMax); *fLog << "Normalization factor computed from alpha plot (after cuts) " << endl << "using counted number of ON and OFF events in alpha region " << endl << "defined by range " << fAlphaBkgMin << "-" << fAlphaBkgMax << endl << "Normalization factor = " << NewNormFactor << endl; *fLog << "Normalization factor used is the one computed in bkg region; " << endl << "i.e. " << NewNormFactor << " instead of " << fNormFactorTest << endl; fNormFactorTest = NewNormFactor; } // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF MHFindSignificanceONOFF findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); findsig.SetUseFittedQuantities(fUseFittedQuantities); TString psfilename; // Name of psfile where Alpha plot will be stored if (!fPsFilenameString.IsNull()) { psfilename += (fPsFilenameString); psfilename += ("TESTSample"); } else { psfilename += ("TESTSample"); } findsig.FindSigmaONOFF(&alphaHist2, &alphaOFFHista, fNormFactorTest, alphamin, alphamax, degree, degreeOFF, alphasig, drawpoly, fitgauss, print, saveplots, psfilename); const Double_t significance = findsig.GetSignificance(); const Double_t alphasi = findsig.GetAlphasi(); // Update values of Nex and SigmaLiMa for Test sammple fSigmaLiMaTest = double(findsig.GetSignificance()); if(fUseFittedQuantities) { fNexTest = double(findsig.GetNexONOFFFitted()); } else { fNexTest = double(findsig.GetNexONOFF()); } // Alpha distributions and containers with teh supercuts applied // are written into root file TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE", "Alpha Distributions for several Theta bins"); alphaHist2.Write(); alphaOFFHista.Write(); (supapptestON.GetTreePointer())->Write(); (supapptestOFF.GetTreePointer())->Write(); *fLog << "TTree objects containing supercuts applied to TEST sample are " << "written into root file " << fAlphaDistributionsRootFilename << endl; (supapptestON.GetTreePointer())->Print(); (supapptestOFF.GetTreePointer())->Print(); //rootfile.Close(); // Boolean variable that controls in class MSupercutsCalcONOFF // wether the supercuts are stored or not // is set to the default value, kFALSE. fCalcHadTest -> SetStoreAppliedSupercuts(kFALSE); fCalcHadTestOFF -> SetStoreAppliedSupercuts(kFALSE); *fLog << "Significance of gamma signal after supercuts in Test sample : " << significance << " (for |alpha| < " << alphasi << " degrees)" << endl; //------------------------------------------- *fLog << "Test of supercuts on TEST sample finished" << endl; *fLog << "======================================================" << endl; return kTRUE; } // Function that applies the optimized supercuts on the TEST sample. // Supercuts applied (LengthUp, LengthLow...) // to all the individual events, as well as // the features of the shower images (Length,Width...) are // stored in 2 TTree objects (for ON and OFF data) in the // root file specified by variable "fAlphaDistributionsRootFilename" Bool_t MFindSupercutsONOFF::TestParamsOnTrainSample() { if (fMatrixTrain->GetM().GetNrows() <= 0) { *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; train matrix ON has no entries" << endl; return kFALSE; } if (fMatrixTrainOFF->GetM().GetNrows() <= 0) { *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; train matrix OFF has no entries" << endl; return kFALSE; } // ------------- TEST the supercuts on TRAIN sample -------------------------- // *fLog << "Producing output of the SUPERCUTS ON-OFF optimization " << endl; // ----------------------------------------------------------------- // read optimum parameter values from file filenameParam // into array 'supercutsPar' TFile inparam(fFilenameParam); MSupercuts scin; scin.Read("MSupercuts"); inparam.Close(); *fLog << "Optimum parameter values for supercuts were read from file '"; *fLog << fFilenameParam << "' :" << endl; const TArrayD &supercutsPar = scin.GetParameters(); for (Int_t i=0; i SetHadronnessName(fHadronnessNameOFF); fCalcHadTrain -> SetHadronnessName(fHadronnessName); // Settings related to the storage of teh supercuts applied fCalcHadTrain -> SetSupercutsAppliedName(fTrainONSupercutsAppliedName); fCalcHadTrain -> SetStoreAppliedSupercuts(kTRUE); fCalcHadTrainOFF -> SetSupercutsAppliedName(fTrainOFFSupercutsAppliedName); fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kTRUE); // Set boolean variable that controls wether cuts are // dynamic or static fCalcHadTrain -> SetVariableUseStaticCuts(fUseStaticCuts); fCalcHadTrainOFF -> SetVariableUseStaticCuts(fUseStaticCuts); if (!fUseStaticCuts) { // Set boolean variable that controls wether the theta variable // is used or not in the computation of the dynamical cuts fCalcHadTrain -> SetVariableNotUseTheta(fNotUseTheta); fCalcHadTrainOFF -> SetVariableNotUseTheta(fNotUseTheta); } // apply the supercuts to OFF data MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5"); MContinue applysupercutsOFF(&scfilterOFF); // apply the supercuts to ON data MF scfilter(fHadronnessName+".fHadronness>0.5"); MContinue applysupercuts(&scfilter); // **************************************************** // Filling OFF alpha distribution MParList parlistOFF; MTaskList tasklistOFF; MMatrixLoop loopOFF(fMatrixTrainOFF); // plot alpha before applying the supercuts const TString mh3NameOFFB = "AlphaOFFBefore"; MBinning binsalphaOFFbef("Binning"+mh3NameOFFB); //binsalphaOFFbef.SetEdges(54, -12.0, 96.0); binsalphaOFFbef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); // |alpha| is assumed to be in column 7 of the matrix MH3 alphaOFFbefore("MatrixTrainOFF[7]"); alphaOFFbefore.SetName(mh3NameOFFB); MFillH fillalphaOFFbefore(&alphaOFFbefore); fillalphaOFFbefore.SetName("FillAlphaOFFBefore"); // fill OFF alpha after applying the supercuts const TString mh3NameOFFA = "AlphaOFFAfter"; MBinning binsalphaOFFaft("Binning"+mh3NameOFFA); //binsalphaOFFaft.SetEdges(54, -12.0, 96.0); binsalphaOFFaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); MH3 alphaOFFafter("MatrixTrainOFF[7]"); alphaOFFafter.SetName(mh3NameOFFA); //TH1 &alphaOFFhista = alphaOFFafter.GetHist(); //alphaOFFhista.SetName("alphaOFFAfter-OptimizationParamOutput"); MFillH fillalphaOFFafter(&alphaOFFafter); fillalphaOFFafter.SetName("FillAlphaOFFAfter"); //****************************** // entries in MParList parlistOFF.AddToList(&tasklistOFF); parlistOFF.AddToList(&supercuts); parlistOFF.AddToList(&supapptrainOFF); parlistOFF.AddToList(fCam); parlistOFF.AddToList(fMatrixTrainOFF); parlistOFF.AddToList(&binsalphaOFFbef); parlistOFF.AddToList(&alphaOFFbefore); parlistOFF.AddToList(&binsalphaOFFaft); parlistOFF.AddToList(&alphaOFFafter); //****************************** // Entries in MtaskList tasklistOFF.AddToList(&loopOFF); tasklistOFF.AddToList(&fillalphaOFFbefore); tasklistOFF.AddToList(fCalcHadTrainOFF); tasklistOFF.AddToList(&applysupercutsOFF); tasklistOFF.AddToList(&fillalphaOFFafter); //****************************** MProgressBar barOFF; MEvtLoop evtloopOFF; evtloopOFF.SetParList(&parlistOFF); evtloopOFF.SetName("EvtLoopOptimizationParamOutputOFF"); evtloopOFF.SetProgressBar(&barOFF); Int_t maxeventsOFF = -1; if (!evtloopOFF.Eventloop(maxeventsOFF)) return kFALSE; tasklistOFF.PrintStatistics(0, kTRUE); //------------------------------------------- // draw the alpha plots MH3* alphaOFFBefore = (MH3*)parlistOFF.FindObject(mh3NameOFFB, "MH3"); TH1 &alphaOFFHistb = alphaOFFBefore->GetHist(); alphaOFFHistb.SetXTitle("|\\alpha| [\\circ]"); TString alphaOFFHistbName ("TrainAlphaOFFBeforeCuts"); if (fThetaRangeString.IsNull()) { *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; " << " fThetaRangeString is NOT defined..." << endl; alphaOFFHistbName += ("ThetaRangeStringUndefined"); } else { alphaOFFHistbName += (fThetaRangeString); } alphaOFFHistb.SetName(alphaOFFHistbName); MH3* alphaOFFAfter = (MH3*)parlistOFF.FindObject(mh3NameOFFA, "MH3"); TH1 &alphaOFFHista = alphaOFFAfter->GetHist(); TString alphaOFFHistaName ("TrainAlphaOFFAfterCuts"); if (fThetaRangeString.IsNull()) { *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; " << " fThetaRangeString is NOT defined..." << endl; alphaOFFHistaName += ("ThetaRangeStringUndefined"); } else { alphaOFFHistaName += (fThetaRangeString); } alphaOFFHista.SetXTitle("|\\alpha| [\\circ]"); alphaOFFHista.SetName(alphaOFFHistaName); // ***************************************************** //------------------------------------------- MParList parlist2; MTaskList tasklist2; MMatrixLoop loop(fMatrixTrain); // plot alpha before applying the supercuts const TString mh3NameB = "AlphaBefore"; MBinning binsalphabef("Binning"+mh3NameB); //binsalphabef.SetEdges(54, -12.0, 96.0); binsalphabef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); // |alpha| is assumed to be in column 7 of the matrix MH3 alphabefore("MatrixTrain[7]"); alphabefore.SetName(mh3NameB); MFillH fillalphabefore(&alphabefore); fillalphabefore.SetName("FillAlphaBefore"); // plot alpha after applying the supercuts const TString mh3NameA = "AlphaAfter"; MBinning binsalphaaft("Binning"+mh3NameA); //binsalphaaft.SetEdges(54, -12.0, 96.0); binsalphaaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp); MH3 alphaafter("MatrixTrain[7]"); alphaafter.SetName(mh3NameA); MFillH fillalphaafter(&alphaafter); fillalphaafter.SetName("FillAlphaAfter"); //****************************** // entries in MParList parlist2.AddToList(&tasklist2); parlist2.AddToList(&supercuts); parlist2.AddToList(&supapptrainON); parlist2.AddToList(fCam); parlist2.AddToList(fMatrixTrain); parlist2.AddToList(&binsalphabef); parlist2.AddToList(&alphabefore); parlist2.AddToList(&binsalphaaft); parlist2.AddToList(&alphaafter); //****************************** // entries in MTaskList tasklist2.AddToList(&loop); tasklist2.AddToList(&fillalphabefore); tasklist2.AddToList(fCalcHadTrain); tasklist2.AddToList(&applysupercuts); tasklist2.AddToList(&fillalphaafter); //****************************** MProgressBar bar2; MEvtLoop evtloop2; evtloop2.SetParList(&parlist2); evtloop2.SetName("EvtLoopOptimizationParamOutput"); evtloop2.SetProgressBar(&bar2); Int_t maxevents2 = -1; if (!evtloop2.Eventloop(maxevents2)) return kFALSE; tasklist2.PrintStatistics(0, kTRUE); //------------------------------------------- // draw the alpha plots MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3"); TH1 &alphaHist1 = alphaBefore->GetHist(); TString alphaHist1Name ("TrainAlphaBeforeCuts"); if (fThetaRangeString.IsNull()) { *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; " << " fThetaRangeString is NOT defined..." << endl; alphaHist1Name += ("ThetaRangeStringUndefined"); } else { alphaHist1Name += (fThetaRangeString); } alphaHist1.SetXTitle("|\\alpha| [\\circ]"); alphaHist1.SetName(alphaHist1Name); MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3"); TH1 &alphaHist2 = alphaAfter->GetHist(); TString alphaHist2Name ("TrainAlphaAfterCuts"); if (fThetaRangeString.IsNull()) { *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; " << " fThetaRangeString is NOT defined..." << endl; alphaHist2Name += ("ThetaRangeStringUndefined"); } else { alphaHist2Name += (fThetaRangeString); } alphaHist2.SetXTitle("|\\alpha| [\\circ]"); alphaHist2.SetName(alphaHist2Name); // fPsFilename -> NewPage(); // Canvas is deleted after saving the plot in ps file // This is to allow for GRID analysis TCanvas *c = new TCanvas("TrainAlphaBeforeAfterSC", "TrainAlphaBeforeAfterSC", 600, 600); c->Divide(2,2); c->cd(1); alphaOFFHistb.DrawCopy(); c->cd(2); alphaOFFHista.DrawCopy(); //c->Modified(); //c-> Update(); // fPsFilename -> NewPage(); c->cd(3); alphaHist1.DrawCopy(); c->cd(4); alphaHist2.DrawCopy(); c -> Modified(); c -> Update(); // ******************************************** // TMP solutions while the TPostScript thing is not working // PsFileName for storing these plots is derived // from PsFilenameString. cout << "Alpha distributions will be saved in PostScript file " ; if (!fPsFilenameString.IsNull()) { TString filename = (fPsFilenameString); filename += ("AlphaTRAINSampleONOFFBeforeAfterCuts.ps"); cout << filename << endl; c -> SaveAs(filename); } // END OF TEMPORAL SOLUTION // ******************************************** // fPsFilename -> NewPage(); // Canvas is deleted after saving the plot in ps file // This is to allow for GRID analysis delete c; //------------------------------------------- // fit alpha distribution to get the number of excess events and // calculate significance of gamma signal in the alpha plot const Double_t alphasig = fAlphaSig; const Double_t alphamin = fAlphaBkgMin; // alpha min for bkg region in ON data const Double_t alphamax = fAlphaBkgMax; const Int_t degree = 2; const Int_t degreeOFF = 2; const Bool_t drawpoly = kTRUE; const Bool_t fitgauss = kTRUE; const Bool_t print = kTRUE; Bool_t saveplots = kTRUE; // Save plots in Psfile if (!ComputeNormFactorTrain()) { *fLog << "Normalization factor for train sample (ON-OFF) could not be computed. Aborting..." << endl; return kFALSE; } // Normalization factor factor is corrected using the estimated // number of excess events and the total number of OFF events // fNormFactor = fNormFactor - Ngammas/EventsInMatrixOFF if (fTuneNormFactor) { // Run MHFindSignificanceONOFF::FindSigmaONOFF // to get estimation of number of excess gamma events. // This number will be used in the correction of the // normalization factor // No plots will be produced this time... const Bool_t drawpolyTest = kFALSE; const Bool_t fitgaussTest = kFALSE; const Bool_t printTest = kFALSE; const Bool_t saveplotsTest = kFALSE; // Save plots in Psfile //TPostScript* DummyPs = new TPostScript("dummy.ps"); TString DummyPs = ("Dummy"); MHFindSignificanceONOFF findsigTest; findsigTest.SetRebin(kTRUE); findsigTest.SetReduceDegree(kFALSE); findsigTest.SetUseFittedQuantities(fUseFittedQuantities); if(findsigTest.FindSigmaONOFF(&alphaHist2, &alphaOFFHista, fNormFactorTrain, alphamin, alphamax, degree, degreeOFF, alphasig, drawpolyTest, fitgaussTest, printTest, saveplotsTest, DummyPs)) { // Update values of Nex and SigmaLiMa for Test sammple fSigmaLiMaTrain = double(findsigTest.GetSignificance()); if(fUseFittedQuantities) { fNexTrain = double(findsigTest.GetNexONOFFFitted()); } else { fNexTrain = double(findsigTest.GetNexONOFF()); } } else { *fLog << "Normalization Factor tuning for Train sample is not possible" << endl; fSigmaLiMaTrain = 0.0; fNexTrain = 0.0; } //DummyPs -> Close(); //DummyPs = NULL; //delete DummyPs; Int_t EventsInTrainMatrixOFF = fMatrixTrainOFF->GetM().GetNrows(); Double_t Ngammas = fNexTrain/fGammaEfficiency; Double_t GammaFraction = Ngammas/EventsInTrainMatrixOFF; *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; " << "fNormFactorTrain is corrected with fraction of gammas in ON " << "Data Train sample and the number of OFF events before cuts." << endl << "EventsInTrainMatrixOFF = " << EventsInTrainMatrixOFF << " Excess events in Train ON sample = " << fNexTrain << " fGammaEfficiency = " << fGammaEfficiency << endl << "fNormFactorTrain Value before correction is " << fNormFactorTrain << endl; fNormFactorTrain = fNormFactorTrain - GammaFraction; *fLog << "fNormFactorTrain Value after correction is " << fNormFactorTrain << endl; } if (fNormFactorFromAlphaBkg) { Double_t NewNormFactor = ComputeNormFactorFromAlphaBkg(&alphaHist2, &alphaOFFHista, fAlphaBkgMin, fAlphaBkgMax); *fLog << "Normalization factor computed from alpha plot (after cuts) " << endl << "using counted number of ON and OFF events in alpha region " << endl << "defined by range " << fAlphaBkgMin << "-" << fAlphaBkgMax << endl << "Normalization factor = " << NewNormFactor << endl; *fLog << "Normalization factor used is the one computed in bkg region; " << endl << "i.e. " << NewNormFactor << " instead of " << fNormFactorTrain << endl; fNormFactorTrain = NewNormFactor; } MHFindSignificanceONOFF findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); findsig.SetUseFittedQuantities(fUseFittedQuantities); TString psfilename; // Name of psfile where Alpha plot will be stored if (!fPsFilenameString.IsNull()) { psfilename += (fPsFilenameString); psfilename += ("TRAINSample"); } else { psfilename += ("TRAINSample"); } findsig.FindSigmaONOFF(&alphaHist2, &alphaOFFHista, fNormFactorTrain, alphamin, alphamax, degree, degreeOFF, alphasig, drawpoly, fitgauss, print, saveplots, psfilename); const Double_t significance = findsig.GetSignificance(); const Double_t alphasi = findsig.GetAlphasi(); // Update Train values Nex, SigmaLiMa fSigmaLiMaTrain = double(findsig.GetSignificance()); if(fUseFittedQuantities) { fNexTrain = double(findsig.GetNexONOFFFitted()); } else { fNexTrain = double(findsig.GetNexONOFF()); } // Alpha distributions and containers with the supercuts applied // are written into root file TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE", "Alpha Distributions for several Theta bins"); alphaHist2.Write(); alphaOFFHista.Write(); (supapptrainON.GetTreePointer())->Write(); (supapptrainOFF.GetTreePointer())->Write(); *fLog << "TTree objects containing supercuts applied to TEST sample are " << "written into root file " << fAlphaDistributionsRootFilename << endl; (supapptrainON.GetTreePointer())->Print(); (supapptrainOFF.GetTreePointer())->Print(); //rootfile.Close(); // Boolean variable that controls in class MSupercutsCalcONOFF // wether the supercuts are stored or not // is set to the default value, kFALSE. fCalcHadTrain -> SetStoreAppliedSupercuts(kFALSE); fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kFALSE); *fLog << "Significance of gamma signal after supercuts in train sample : " << significance << " (for |alpha| < " << alphasi << " degrees)" << endl; //------------------------------------------- *fLog << "Test of supercuts on TRAIN sample finished" << endl; *fLog << "======================================================" << endl; return kTRUE; }