/* ======================================================================== *\ ! ! * ! * 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 ! Author(s): Marcos Lopez, 05/2004 ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MMyFindSuperCuts // // // // Class for otimizing the parameters of the supercuts // // // // // // // ///////////////////////////////////////////////////////////////////////////// #include "MMyFindSuperCuts.h" #include // fabs #include #include #include #include #include #include #include "MBinning.h" #include "MContinue.h" #include "MMySuperCuts.h" #include "MMySuperCutsCalc.h" #include "MDataElement.h" #include "MDataMember.h" #include "MEvtLoop.h" #include "MFSelFinal.h" #include "MF.h" #include "MFEventSelector.h" #include "MFEventSelector2.h" #include "MFillH.h" //#include "MGeomCamCT1Daniel.h" #include "MFEventSelector.h" #include "MGeomCamMagic.h" #include "MH3.h" #include "MHSupercuts.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" #include "TVector2.h" #include "MSrcPosCam.h" #include "MHillasSrcCalc.h" #include "MSrcPosCalc.h" #include "MObservatory.h" ClassImp(MMyFindSuperCuts); 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) // // ------------------------- // - Read the Matrix // - Fill a histogram with the alpha (6th column) // - Fit the alpha plot // - Get result of fit (significance) // static void fcnSuperCuts(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) { // // Set by hand the paramters range // if ( par[0]<0 || par[0]>1.8 || par[4]>par[0] || par[4]<0 || par[0]1 || par[3]>1 || par[4]>1 || par[1]>1 || par[6]>1 || par[7]>1 || par[8]<0 || par[8]>1.8 || par[12]>par[0] || par[12]<0 || par[8]1. || par[14]<0 || par[14]>1 || par[15]<0 || par[15]>1 || par[16]<0 || par[16]>1.8 || par[20]<0 || par[20]>1.8 || par[16]0.18 || par[24]>200 || par[28]<-200 || par[32]>1 || par[36]<0 || par[40]>1 ) { gLog << "Parameter out of limit" << endl; f=1e10; cout << f << endl; return; } //------------------------------------------------------------- // 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(); MMySuperCuts *super = (MMySuperCuts*)plistfcn->FindObject("MMySuperCuts"); if (!super) { gLog << "fcnSupercuts : MMySuperCuts object '" << "MMySuperCuts" << "' not found... aborting" << endl; return; } // // transfer current parameter values to MMySuperCuts // // 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 gLog.Separator("Current params"); TArrayD checkparameters = super->GetParameters(); gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", " << fNparx << endl; // gLog << "fcnsupercuts : i, par, checkparameters =" << endl; // for (Int_t i=0; iPrint(); //$$$$$$$$$$$$$$$$$$$$$ // // EventLoop: Fill the alpha plot with these cuts // gLog.SetNullOutput(kTRUE); evtloopfcn->Eventloop(); //tasklistfcn->PrintStatistics(0, kTRUE); // // Calculate significance // 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 = 10.0; const Double_t alphamin = 30.0; const Double_t alphamax = 90.0; const Int_t degree = 0; Bool_t drawpoly; Bool_t fitgauss; TCanvas* c; if (iflag == 3) { c= new TCanvas(); drawpoly = kTRUE; fitgauss = kTRUE; } else { drawpoly = kFALSE; fitgauss = kFALSE; } fitgauss = kTRUE; const Bool_t print = kFALSE; MHFindSignificance findsig; //findsig.SetRebin(kTRUE); findsig.SetRebin(kFALSE); findsig.SetReduceDegree(kFALSE); const Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, alphasig, drawpoly, fitgauss, print); gLog.SetNullOutput(kFALSE); //================================================================= // 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 MHSupercuts *plotsuper = (MHSupercuts*)plistfcn->FindObject("MHSupercuts"); if (plotsuper) plotsuper->Fill(&findsig); // // Print info // gLog << "(Non, Nbg, Nex, S/N, sigma, Significane) = " << findsig.GetNon() << " " << findsig.GetNbg()<< " " << findsig.GetNex() << " " << findsig.GetNex()/ findsig.GetNbg()<< " " << findsig.GetSigmaGauss() << " " << findsig.GetSignificance() << endl; const Double_t sigmagauss = findsig.GetSigmaGauss(); // if (sigmagauss>10) // { // gLog << "fcnSupercuts : Alpha plot too width = " << sigmagauss << endl; // f = 1.e10; // return; // } // // Get significance // const Double_t significance = findsig.GetSignificance(); //const Double_t ratio = findsig.GetNex()/findsig.GetNbg(); // f = significance>0 ? -significance : 0; // f = -TMath::Abs(significance)/sqrt(sigmagauss); // f = - significance*sqrt(findsig.GetNex())/sqrt(sigmagauss); f = - significance*significance/sqrt(sigmagauss); // f = - significance*findsig.GetProb(); cout << f << endl; // // // // optimize signal/background ratio // Double_t ratio = findsig.GetNbg()>0.0 ? // findsig.GetNex()/findsig.GetNbg() : 0.0; // //f = -ratio; // if(significance < 5) // f=1e10; // f=f*ratio*ratio; // cout << f << " " << significance << " " << ratio << endl; } // -------------------------------------------------------------------------- // // Default constructor. // MMyFindSuperCuts::MMyFindSuperCuts(const char *name, const char *title) : fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL),fSizeCutLow(2000.0),fSizeCutUp(10000000), fOptimizationMode(0) { fName = name ? name : "MMyFindSuperCuts"; fTitle = title ? title : "Optimizer of the supercuts"; //--------------------------- // camera geometry is needed for conversion mm ==> degree //fCam = new MGeomCamCT1Daniel; fCam = new MGeomCamMagic; // matrices to contain the training/test samples fMatrixTrain = new MHMatrix("MatrixTrain"); fMatrixTest = new MHMatrix("MatrixTest"); // objects of MMySuperCutsCalc to which these matrices are attached fCalcHadTrain = new MMySuperCutsCalc("SupercutsCalcTrain"); fCalcHadTest = new MMySuperCutsCalc("SupercutsCalcTest"); // Define columns of matrices *fLog << inf << "Init Mapping of Train Matrix" << endl; fCalcHadTrain->InitMapping(fMatrixTrain); *fLog << inf << "Init Mapping of Test Matrix" << endl; fCalcHadTest->InitMapping(fMatrixTest); } // -------------------------------------------------------------------------- // // Default destructor. // MMyFindSuperCuts::~MMyFindSuperCuts() { 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 MMyFindSuperCuts::DefineTrainMatrix( const TString &filenametrain, MH3 &hreftrain, const Int_t howmanytrain, const TString &outfiletrain) { if (filenametrain.IsNull() || howmanytrain <= 0) return kFALSE; *fLog << endl << "***************************************************************" << endl << " Filling training matrix from file: '" << filenametrain << "'" << endl << " Selecting: " << howmanytrain << " events"; if (!hreftrain.GetHist().GetEntries()==0) { *fLog << ", according to a distribution given by the MH3 object '" << hreftrain.GetName() << "'" << endl; } else { *fLog << ", Randomly" << endl; } *fLog << "***************************************************************" << endl; // // ParList // MParList plist; MTaskList tlist; plist.AddToList(&tlist); // entries in MParList plist.AddToList(fCam); plist.AddToList(fMatrixTrain); // // TaskList // MReadTree read("Events", filenametrain); //MReadTree read("Parameters", filenametrain); read.DisableAutoScheme(); // MFEventSelector2 seltrain(hreftrain); // seltrain.SetNumMax(howmanytrain); // seltrain.SetName("selectTrain"); // Random selection of events MFEventSelector seltrain; seltrain.SetNumSelectEvts(howmanytrain); seltrain.SetName("selectTrain"); MFillH filltrain(fMatrixTrain); filltrain.SetFilter(&seltrain); filltrain.SetName("fillMatrixTrain"); // // ReCalculate hillass // MObservatory obs; // plist.AddToList(&obs); // MSrcPosCalc srccalc; // srccalc.SetPositionXY(3.46945e-17, 0.0487805); // MHillasSrcCalc hsrc; // entries in MTaskList tlist.AddToList(&read); tlist.AddToList(&seltrain); // tlist.AddToList(&srccalc); // tlist.AddToList(&hsrc); tlist.AddToList(&filltrain); // // EventLoop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTrain"); evtloop.SetProgressBar(&bar); if (!evtloop.Eventloop()) return kFALSE; // if (!evtloop.PreProcess()) // return kFALSE; // Int_t i=0; // while(tlist.Process()) // { // i++; // if (i>howmanytrain) // break; // } // evtloop.PostProcess(); tlist.PrintStatistics(0, kTRUE); // // Print // fMatrixTrain->Print("SizeCols"); Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows(); if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MMyFindSuperCuts::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 (outfiletrain != "") { TFile filetr(outfiletrain, "RECREATE"); fMatrixTrain->Write(); filetr.Close(); *fLog << "MMyFindSuperCuts::DefineTrainMatrix; Training matrix was written onto file '" << outfiletrain << "'" << endl; } // fMatrixTrain->Print("data"); 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 MMyFindSuperCuts::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 (!hreftest.GetHist().GetEntries()==0) { *fLog << " according to a distribution given by the MH3 object '" << hreftest.GetName() << "'" << endl; } else { *fLog << " randomly" << endl; } // // ParList // MParList plist; MTaskList tlist; // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTest); // // TaskList // MReadTree read("Parameters", nametest); read.DisableAutoScheme(); // MFEventSelector2 seltest(hreftest); // seltest.SetNumMax(howmanytest); // seltest.SetName("selectTest"); MFillH filltest(fMatrixTest); //filltest.SetFilter(&seltest); // entries in MTaskList tlist.AddToList(&read); //tlist.AddToList(&seltest); tlist.AddToList(&filltest); // // EventLoop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTest"); evtloop.SetProgressBar(&bar); if (!evtloop.PreProcess()) return kFALSE; Int_t i=0; while(tlist.Process()) { i++; if (i>100000) break; } evtloop.PostProcess(); // if (!evtloop.Eventloop()) // return kFALSE; tlist.PrintStatistics(0, kTRUE); // // Print // fMatrixTest->Print("SizeCols"); const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows(); if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MMyFindSuperCuts::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 << "MMyFindSuperCuts::DefineTestMatrix; Test matrix was written onto file '" << filetest << "'" << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // Define the matrices for the training and test sample respectively // // // Bool_t MMyFindSuperCuts::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 (!href.GetHist().GetEntries()==0) { *fLog << " according to a distribution given by the MH3 object '" << href.GetName() << "'" << endl; } else { *fLog << " randomly" << endl; } // // ParList // MParList plist; MTaskList tlist; // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTrain); plist.AddToList(fMatrixTest); // // TaskList // MReadMarsFile read("Events", name); read.DisableAutoScheme(); MFEventSelector2 selector(href); selector.SetNumMax(howmanytrain+howmanytest); selector.SetName("selectTrainTest"); selector.SetInverted(); MContinue cont(&selector); cont.SetName("ContTrainTest"); Double_t prob = ( (Double_t) howmanytrain ) / ( (Double_t)(howmanytrain+howmanytest) ); MFEventSelector split; split.SetSelectionRatio(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 MTaskList tlist.AddToList(&read); tlist.AddToList(&cont); tlist.AddToList(&split); tlist.AddToList(&filltrain); tlist.AddToList(&conttrain); tlist.AddToList(&filltest); // // EventLoop // 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 << "MMyFindSuperCuts::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 << "MMyFindSuperCuts::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 << "MMyFindSuperCuts::DefineTrainTestMatrix; Training matrix was written onto file '" << filetrain << "'" << endl; } // // Write out test matrix // if (filetest != "") { TFile filete(filetest, "RECREATE", ""); fMatrixTest->Write(); filete.Close(); *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; Test matrix was written onto file '" << filetest << "'" << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // Read training and test matrices from files // Bool_t MMyFindSuperCuts::ReadMatrix(const TString &filetrain, const TString &filetest) { // // Read in training matrix // TFile filetr(filetrain); fMatrixTrain->Read("MatrixTrain"); fMatrixTrain->Print("SizeCols"); *fLog << "MMyFindSuperCuts::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 <<"MMyFindSuperCuts::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 MMySuperCutsCalc 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 MMyFindSuperCuts::FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps ) { fCalcHadTrain->SetSizeCuts(fSizeCutLow,fSizeCutUp); fCalcHadTest->SetSizeCuts(fSizeCutLow,fSizeCutUp); // 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; // // Check that all the necessary containers already exists // if (fHadronnessName.IsNull()) { *fLog << "MyFindSuperCuts::FindParams; hadronness name is not defined... aborting" << endl; return kFALSE; } if (fMatrixTrain == NULL) { *fLog << "MMyFindSuperCuts::FindParams; training matrix is not defined... aborting" << endl; return kFALSE; } if (fMatrixTrain->GetM().GetNrows() <= 0) { *fLog << "MMyFindSuperCuts::FindParams; training matrix has no entries" << endl; return kFALSE; } // // ParList // MParList parlistfcn; MTaskList tasklistfcn; parlistfcn.AddToList(&tasklistfcn); // // create container for the supercut parameters // and set them to their initial values // MMySuperCuts super; // take initial values from file parSCinit if (parSCinit != "") { TFile inparam(parSCinit); super.Read("MMySuperCuts"); inparam.Close(); *fLog << "MMyFindSuperCuts::FindParams; initial values of parameters are taken from file " << parSCinit << endl; super.Print(); } // take initial values from 'params' and/or 'steps' //else if (params.GetSize() != 0 || steps.GetSize() != 0 ) if (params.GetSize() != 0 || steps.GetSize() != 0 ) { if (params.GetSize() != 0) { *fLog << "MMyFindSuperCuts::FindParams; initial values of parameters are taken from 'params'" << endl; super.SetParameters(params); } if (steps.GetSize() != 0) { *fLog << "MMyFindSuperCuts::FindParams; initial step sizes are taken from 'steps'" << endl; super.SetStepsizes(steps); } } else { *fLog << "MMyFindSuperCuts::FindParams; initial values and step sizes are taken from the MMySuperCuts constructor" << endl; } // Alpha binning const TString mh3Name = "AlphaFcn"; MBinning binsalpha("Binning"+mh3Name); // binsalpha.SetEdges(54, -12.0, 96.0); binsalpha.SetEdges(45, 0.0, 90.0); // Alpha plot |alpha| MH3 alpha("MatrixTrain[6]"); // |alpha| assumed to be in column 6 alpha.SetName(mh3Name); *fLog << warn << "WARNING----->ALPHA IS ASSUMED TO BE IN COLUMN 6!!!!!!!!" << endl; // book histograms to be filled during the optimization // ! not in the event loop ! MHSupercuts plotsuper; plotsuper.SetupFill(&parlistfcn); // Add to ParList parlistfcn.AddToList(&super); parlistfcn.AddToList(fCam); parlistfcn.AddToList(fMatrixTrain); parlistfcn.AddToList(&binsalpha); parlistfcn.AddToList(&alpha); parlistfcn.AddToList(&plotsuper); // // MTaskList // MMatrixLoop loop(fMatrixTrain); // loop over rows of matrix // calculate supercuts hadronness //fCalcHadTrain->SetHadronnessName(fHadronnessName); // apply the supercuts MF scfilter(fHadronnessName+".fHadronness>0.5"); MContinue supercuts(&scfilter); // Fill Alpha Plot MFillH fillalpha(&alpha); // Add to TaskList tasklistfcn.AddToList(&loop); tasklistfcn.AddToList(fCalcHadTrain); tasklistfcn.AddToList(&supercuts); tasklistfcn.AddToList(&fillalpha); // // EventLoop // MEvtLoop evtloopfcn("EvtLoopFCN"); evtloopfcn.SetParList(&parlistfcn); *fLog << "Event loop for fcnSupercuts has been setup" << endl; //----------------------------------------------------------------------- // //---------- 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; // } // // //} // // vary only first 23 parameters // // // for (UInt_t i=0; i= 24) // { // 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 MMySuperCuts // // 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 << "MMyFindSuperCuts::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); MMySuperCuts scin; scin.Read("MMySuperCuts"); 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[6]"); 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[6]"); 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; }