/* ======================================================================== *\ ! ! * ! * 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 ! 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 "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) { //------------------------------------------------------------- // 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(); 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 : i, par, checkparameters =" << endl; //for (Int_t i=0; iEventloop(); MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3"); if (!alpha) return; TH1 &alphaHist = alpha->GetHist(); //alphaHist.SetXTitle("|\\alpha| [\\circ]"); //------------------------------------------- // 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 = 4; 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"; //--------------------------- // camera geometry is needed for conversion mm ==> degree fCam = new MGeomCamCT1Daniel; // matrices to contain the 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; } // -------------------------------------------------------------------------- // // 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, const Int_t howmanytrain, MH3 &hreftrain) { if (nametrain.IsNull() || howmanytrain <= 0) return kFALSE; *fLog << "=============================================" << endl; *fLog << "fill training matrix from file '" << nametrain << "', select " << howmanytrain << " events according to a distribution given by the MH3 object '" << hreftrain.GetName() << "'" << endl; MParList plist; MTaskList tlist; MReadMarsFile read("Events", nametrain); read.DisableAutoScheme(); //MFEventSelector2 seltrain(hreftrain); //seltrain.SetNumMax(howmanytrain); MFEventSelector seltrain; seltrain.SetNumSelectEvts(howmanytrain); MFillH filltrain(fMatrixTrain); filltrain.SetFilter(&seltrain); //****************************** // 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; 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, const Int_t howmanytest, MH3 &hreftest) { if (nametest.IsNull() || howmanytest<=0) return kFALSE; *fLog << "=============================================" << endl; *fLog << "fill test matrix from file '" << nametest << "', select " << howmanytest << " events according to a distribution given by the MH3 object '" << hreftest.GetName() << "'" << endl; MParList plist; MTaskList tlist; MReadMarsFile read("Events", nametest); read.DisableAutoScheme(); //MFEventSelector2 seltest(hreftest); //seltest.SetNumMax(howmanytest); MFEventSelector seltest; seltest.SetNumSelectEvts(howmanytest); 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; return kTRUE; } // -------------------------------------------------------------------------- // // Define the matrices for the training and test sample respectively // // // Bool_t MCT1FindSupercuts::DefineTrainTestMatrix( const TString &name, const Int_t howmanytrain, MH3 &hreftrain, const Int_t howmanytest, MH3 &hreftest) { *fLog << "=============================================" << endl; *fLog << "fill training matrix from file '" << name << "', select " << howmanytrain << " events for the training according to a distribution given by the MH3 object '" << hreftrain.GetName() << "'" << endl; *fLog << "fill test matrix from the same file '" << name << "', select " << howmanytest << " events for the training according to a distribution given by the MH3 object '" << hreftest.GetName() << "'" << endl; MParList plist; MTaskList tlist; MReadMarsFile read("Events", name); read.DisableAutoScheme(); //MFEventSelector2 seltrain(hreftrain); //seltrain.SetNumMax(howmanytrain); MFEventSelector seltrain; seltrain.SetName("SelTrain"); seltrain.SetNumSelectEvts(howmanytrain); MFillH filltrain(fMatrixTrain); filltrain.SetName("fillMatrixTrain"); filltrain.SetFilter(&seltrain); // consider this event as candidate for a test event // only if event was not accepted as a training event MContinue cont(&seltrain); const Float_t fcorr = (Float_t)read.GetEntries() / (read.GetEntries()-howmanytrain); *fLog << "entries, howmanytrain, fcorr = " << read.GetEntries() << ", " << howmanytrain << ", " << fcorr << endl; //MFEventSelector2 seltest(hreftest); //seltrain.SetNumMax(howmanytest*fcorr); MFEventSelector seltest; seltest.SetName("SelTest"); seltest.SetNumSelectEvts((Int_t)(fcorr*howmanytest)); MFillH filltest(fMatrixTest); filltest.SetName("fillMatrixTest"); filltest.SetFilter(&seltest); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTrain); plist.AddToList(fMatrixTest); //****************************** // entries in MTaskList tlist.AddToList(&read); tlist.AddToList(&seltrain); tlist.AddToList(&filltrain); tlist.AddToList(&cont); tlist.AddToList(&seltest); tlist.AddToList(&filltest); //****************************** MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTrainTest"); evtloop.SetProgressBar(&bar); if (!evtloop.Eventloop()) 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::DefineTrainMatrix; 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::DefineTestMatrix; 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; return kTRUE; } // -------------------------------------------------------------------------- // // Read training and test matrices from files // // // $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$ // 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; //-------------------------- // 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; return kTRUE; } // -------------------------------------------------------------------------- // // Write training and test matrices onto files // // // $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$ // // Bool_t MCT1FindSupercuts::WriteMatrix(const TString &filetrain, const TString &filetest) { //-------------------------- // write out training matrix TFile filetr(filetrain, "RECREATE"); *fLog << "nach TFile" << endl; fMatrixTrain->Write(); *fLog << "MCT1FindSupercuts::WriteMatrix; Training matrix was written onto file '" << filetrain << "'" << endl; //-------------------------- // write out test matrix TFile filete(filetest, "RECREATE"); fMatrixTest->Print("SizeCols"); fMatrixTest->Write(); *fLog << "MCT1FindSupercuts::WriteMatrix; Test matrix was written onto file '" << filetest << "'" << endl; 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 as // MCT1Supercuts::GetParameters(fVinit) // //---------------------------------------------------------------------- Bool_t MCT1FindSupercuts::FindParams() { // Setup the event loop which will be executed in the // fcnSupercuts function of MINUIT // *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 default values MCT1Supercuts super; // 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 from MCT1Supercuts fVinit = super.GetParameters(); 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); 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); 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]"); 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 = 4; const Bool_t drawpoly = kTRUE; const Bool_t fitgauss = kFALSE; 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; }