/* ======================================================================== *\ ! ! * ! * 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): Thomas Bretz, 6/2003 ! Author(s): Wolfgang Wittek, 7/2003 ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MCT1FindSupercuts // // // // Class for otimizing the parameters of the supercuts // // // // // // // ///////////////////////////////////////////////////////////////////////////// #include "MCT1FindSupercuts.h" #include // fabs #include #include #include #include #include #include #include "MBinning.h" #include "MContinue.h" #include "MCT1Supercuts.h" #include "MCT1SupercutsCalc.h" #include "MDataElement.h" #include "MDataMember.h" #include "MEvtLoop.h" #include "MFCT1SelFinal.h" #include "MF.h" #include "MFEventSelector.h" #include "MFEventSelector2.h" #include "MFillH.h" //#include "MGeomCamCT1Daniel.h" #include "MFRandomSplit.h" #include "MGeomCamCT1.h" #include "MH3.h" #include "MHCT1Supercuts.h" #include "MHFindSignificance.h" #include "MHMatrix.h" #include "MHOnSubtraction.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(MCT1FindSupercuts); using namespace std; //------------------------------------------------------------------------ // // 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(); MParList *plistfcn = (MParList*)evtloopfcn->GetParList(); //MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList"); MCT1Supercuts *super = (MCT1Supercuts*)plistfcn->FindObject("MCT1Supercuts"); if (!super) { gLog << "fcnSupercuts : MCT1Supercuts object '" << "MCT1Supercuts" << "' not found... aborting" << endl; return; } // // transfer current parameter values to MCT1Supercuts // // 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(); //tasklistfcn->PrintStatistics(0, kTRUE); MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3"); if (!alpha) return; TH1 &alphaHist = alpha->GetHist(); alphaHist.SetName("alpha-fcnSupercuts"); //------------------------------------------- // 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 = 20.0; const Double_t alphamin = 30.0; const Double_t alphamax = 90.0; const Int_t degree = 2; Bool_t drawpoly; Bool_t fitgauss; if (iflag == 3) { drawpoly = kTRUE; fitgauss = kTRUE; } else { drawpoly = kFALSE; fitgauss = kFALSE; } //drawpoly = kFALSE; //fitgauss = kFALSE; const Bool_t print = kTRUE; MHFindSignificance findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); const Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, alphasig, drawpoly, fitgauss, print); //================================================================= // reset gMinuit to the MINUIT object for optimizing the supercuts gMinuit = savePointer; //------------------------------------------- if (!rc) { gLog << "fcnSupercuts : FindSigma() failed" << endl; f = 1.e10; return; } // plot some quantities during the optimization MHCT1Supercuts *plotsuper = (MHCT1Supercuts*)plistfcn->FindObject("MHCT1Supercuts"); if (plotsuper) plotsuper->Fill(&findsig); //------------------------ // get significance const Double_t significance = findsig.GetSignificance(); f = significance>0 ? -significance : 0; //------------------------ // 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. // MCT1FindSupercuts::MCT1FindSupercuts(const char *name, const char *title) : fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL) { fName = name ? name : "MCT1FindSupercuts"; fTitle = title ? title : "Optimizer of the supercuts"; fUseOrigDistribution = kFALSE; //--------------------------- // camera geometry is needed for conversion mm ==> degree //fCam = new MGeomCamCT1Daniel; fCam = new MGeomCamCT1; // matrices to contain the training/test samples fMatrixTrain = new MHMatrix("MatrixTrain"); fMatrixTest = new MHMatrix("MatrixTest"); // objects of MCT1SupercutsCalc to which these matrices are attached fCalcHadTrain = new MCT1SupercutsCalc("SupercutsCalcTrain"); fCalcHadTest = new MCT1SupercutsCalc("SupercutsCalcTest"); // Define columns of matrices fCalcHadTrain->InitMapping(fMatrixTrain); fCalcHadTest->InitMapping(fMatrixTest); } // -------------------------------------------------------------------------- // // Default destructor. // MCT1FindSupercuts::~MCT1FindSupercuts() { delete fCam; delete fMatrixTrain; delete fMatrixTest; delete fCalcHadTrain; delete fCalcHadTest; } // -------------------------------------------------------------------------- // // // void MCT1FindSupercuts::SetUseOrigDistribution(Bool_t b) { fUseOrigDistribution = b; *fLog << "MCT1FindSupercuts : when defining training and test matrices use the original distribution" << endl; } // -------------------------------------------------------------------------- // // Define the matrix 'fMatrixTrain' for the training sample // // alltogether 'howmanytrain' events are read from file 'nametrain'; // the events are selected according to a target distribution 'hreftrain' // // Bool_t MCT1FindSupercuts::DefineTrainMatrix( const TString &nametrain, MH3 &hreftrain, const Int_t howmanytrain, const TString &filetrain) { if (nametrain.IsNull() || howmanytrain <= 0) return kFALSE; *fLog << "=============================================" << endl; *fLog << "fill training matrix from file '" << nametrain << "', select " << howmanytrain << " events " << endl; if (!fUseOrigDistribution) { *fLog << " according to a distribution given by the MH3 object '" << hreftrain.GetName() << "'" << endl; } else { *fLog << " randomly" << endl; } MParList plist; MTaskList tlist; MReadMarsFile read("Events", nametrain); read.DisableAutoScheme(); MFEventSelector2 seltrain(hreftrain); seltrain.SetNumMax(howmanytrain); seltrain.SetName("selectTrain"); if (fUseOrigDistribution) seltrain.SetUseOrigDistribution(kTRUE); MFillH filltrain(fMatrixTrain); filltrain.SetFilter(&seltrain); filltrain.SetName("fillMatrixTrain"); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTrain); //****************************** // entries in MTaskList tlist.AddToList(&read); tlist.AddToList(&seltrain); tlist.AddToList(&filltrain); //****************************** MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTrain"); evtloop.SetProgressBar(&bar); if (!evtloop.Eventloop()) return kFALSE; tlist.PrintStatistics(0, kTRUE); fMatrixTrain->Print("SizeCols"); Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows(); if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MCT1FindSupercuts::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 << "MCT1FindSupercuts::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 MCT1FindSupercuts::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) seltest.SetUseOrigDistribution(kTRUE); MFillH filltest(fMatrixTest); filltest.SetFilter(&seltest); //****************************** // 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 << "MCT1FindSupercuts::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 << "MCT1FindSupercuts::DefineTestMatrix; Test matrix was written onto file '" << filetest << "'" << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // Define the matrices for the training and test sample respectively // // // Bool_t MCT1FindSupercuts::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) 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 << "MCT1FindSupercuts::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 << "MCT1FindSupercuts::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 << "MCT1FindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '" << filetrain << "'" << endl; } //-------------------------- // write out test matrix if (filetest != "") { TFile filete(filetest, "RECREATE", ""); fMatrixTest->Write(); filete.Close(); *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '" << filetest << "'" << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // Read training and test matrices from files // // Bool_t MCT1FindSupercuts::ReadMatrix(const TString &filetrain, const TString &filetest) { //-------------------------- // read in training matrix TFile filetr(filetrain); fMatrixTrain->Read("MatrixTrain"); fMatrixTrain->Print("SizeCols"); *fLog << "MCT1FindSupercuts::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 << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '" << filetest << "'" << endl; filete.Close(); 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 event loop 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 alpha plot is produced in the event loop, // in which the cuts have been applied // - the significance of the gamma signal in the alpha plot // is calculated // // 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 MCT1SupercutsCalc will // put the supercuts hadronness // // - 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 MCT1FindSupercuts::FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps) { // Setup the event loop 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 << "MCT1FindSupercuts::FindParams; hadronness name is not defined... aborting" << endl; return kFALSE; } if (fMatrixTrain == NULL) { *fLog << "MCT1FindSupercuts::FindParams; training matrix is not defined... aborting" << endl; return kFALSE; } if (fMatrixTrain->GetM().GetNrows() <= 0) { *fLog << "MCT1FindSupercuts::FindParams; training matrix has no entries" << endl; return kFALSE; } MParList parlistfcn; MTaskList tasklistfcn; // loop over rows of matrix MMatrixLoop loop(fMatrixTrain); //-------------------------------- // create container for the supercut parameters // and set them to their initial values MCT1Supercuts super; // take initial values from file parSCinit if (parSCinit != "") { TFile inparam(parSCinit); super.Read("MCT1Supercuts"); inparam.Close(); *fLog << "MCT1FindSupercuts::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 << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from 'params'" << endl; super.SetParameters(params); } if (steps.GetSize() != 0) { *fLog << "MCT1FindSupercuts::FindParams; initial step sizes are taken from 'steps'" << endl; super.SetStepsizes(steps); } } else { *fLog << "MCT1FindSupercuts::FindParams; initial values and step sizes are taken from the MCT1Supercits constructor" << endl; } //-------------------------------- // calculate supercuts hadronness fCalcHadTrain->SetHadronnessName(fHadronnessName); // 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); *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(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); //****************************** MEvtLoop evtloopfcn("EvtLoopFCN"); evtloopfcn.SetParList(&parlistfcn); *fLog << "Event loop for fcnSupercuts has been setup" << endl; // 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 TStopwatch clock; clock.Start(); *fLog << "before calling CallMinuit" << endl; MMinuitInterface inter; Bool_t rc = inter.CallMinuit(fcnSupercuts, name, fVinit, fStep, fLimlo, fLimup, fFix, &evtloopfcn, "SIMPLEX", kFALSE); *fLog << "after calling CallMinuit" << endl; *fLog << "Time spent for the minimization in MINUIT : " << endl;; clock.Stop(); clock.Print(); plotsuper.DrawClone(); if (!rc) return kFALSE; *fLog << "Minimization for supercuts finished" << endl; *fLog << "========================================================" << endl; // ----------------------------------------------------------------- // in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put // into MCT1Supercuts // write optimum parameter values onto file filenameParam TFile outparam(fFilenameParam, "RECREATE"); super.Write(); outparam.Close(); *fLog << "Optimum parameter values for supercuts were written onto file '" << fFilenameParam << "' :" << endl; const TArrayD &check = super.GetParameters(); for (Int_t i=0; iGetM().GetNrows() <= 0) { *fLog << "MCT1FindSupercuts::TestParams; test matrix has no entries" << endl; return kFALSE; } // ------------- TEST the supercuts ------------------------------ // *fLog << "Test the supercuts on the test sample" << endl; // ----------------------------------------------------------------- // read optimum parameter values from file filenameParam // into array 'supercutsPar' TFile inparam(fFilenameParam); MCT1Supercuts scin; scin.Read("MCT1Supercuts"); 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; iSetHadronnessName(fHadronnessName); //------------------------------------------- MMatrixLoop loop(fMatrixTest); // plot alpha before applying the supercuts const TString mh3NameB = "AlphaBefore"; MBinning binsalphabef("Binning"+mh3NameB); binsalphabef.SetEdges(54, -12.0, 96.0); // |alpha| is assumed to be in column 7 of the matrix MH3 alphabefore("MatrixTest[7]"); alphabefore.SetName(mh3NameB); TH1 &alphahistb = alphabefore.GetHist(); alphahistb.SetName("alphaBefore-TestParams"); MFillH fillalphabefore(&alphabefore); fillalphabefore.SetName("FillAlphaBefore"); // apply the supercuts MF scfilter(fHadronnessName+".fHadronness>0.5"); MContinue applysupercuts(&scfilter); // plot alpha after applying the supercuts const TString mh3NameA = "AlphaAfter"; MBinning binsalphaaft("Binning"+mh3NameA); binsalphaaft.SetEdges(54, -12.0, 96.0); MH3 alphaafter("MatrixTest[7]"); alphaafter.SetName(mh3NameA); TH1 &alphahista = alphaafter.GetHist(); alphahista.SetName("alphaAfter-TestParams"); MFillH fillalphaafter(&alphaafter); fillalphaafter.SetName("FillAlphaAfter"); //****************************** // entries in MParList parlist2.AddToList(&tasklist2); parlist2.AddToList(&supercuts); 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(); alphaHist1.SetXTitle("|\\alpha| [\\circ]"); MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3"); TH1 &alphaHist2 = alphaAfter->GetHist(); alphaHist2.SetXTitle("|\\alpha| [\\circ]"); alphaHist2.SetName("alpha-TestParams"); TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300); c->Divide(2,1); gROOT->SetSelectedPad(NULL); c->cd(1); alphaHist1.DrawCopy(); c->cd(2); alphaHist2.DrawCopy(); //------------------------------------------- // fit alpha distribution to get the number of excess events and // calculate significance of gamma signal in the alpha plot const Double_t alphasig = 20.0; const Double_t alphamin = 30.0; const Double_t alphamax = 90.0; const Int_t degree = 2; const Bool_t drawpoly = kTRUE; const Bool_t fitgauss = kTRUE; const Bool_t print = kTRUE; MHFindSignificance findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree, alphasig, drawpoly, fitgauss, print); const Double_t significance = findsig.GetSignificance(); const Double_t alphasi = findsig.GetAlphasi(); *fLog << "Significance of gamma signal after supercuts in test sample : " << significance << " (for |alpha| < " << alphasi << " degrees)" << endl; //------------------------------------------- *fLog << "Test of supercuts part finished" << endl; *fLog << "======================================================" << endl; return kTRUE; }