Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 2299) +++ /trunk/MagicSoft/Mars/Changelog (revision 2300) @@ -1,3 +1,50 @@ -*-*- END OF LINE -*-*- + + + + 2003/08/19: Wolfgang Wittek + + * manalysis/MCT1Supercuts.[h,cc] + - new class + - container for the supercut parameters + + * manalysis/MCT1SupercutsCalc.[h,cc] + - get supercut parameters from container 'MCT1Supercuts' + + * manalysis/MCT1FindSupercuts.[h,cc] + - new class + - optimizes the parameters for the supercuts + + * manalysis/MMinuitInterface.[h,cc] + - new class + - interface for Minuit + + * manalysis/Makefile + AnalysisLinkDef.h + - include MCT1FindSupercuts + MMinuitInterface + + * mhist/MH3.cc + - reset fHist in SetupFill(); + this is necessary if the same MH3 object is used in more than one + eventloop + + * mhist/MHMatrix.cc + - give name to the event loop + + * mhist/MHFindSignificance.[h,cc] + - new class + - calculates the significance of the gamma signal in the alpha plot + + * mhist/MHCT1Supercuts.[h,cc] + - new class + - plots various quantities during the optimization of the supercuts + + * mhist/Makefile + HistLinkDef.h + - MHFindSignificance included + - MHCT1Supercuts included + + 2003/08/01: Thomas Bretz @@ -162,4 +209,6 @@ * mgeom/Makefile, GeomLinkDef.h - add new class + + Index: /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2299) +++ /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2300) @@ -51,7 +51,19 @@ #pragma link C++ class MMcTriggerLvl2Calc+; +#pragma link C++ class MCT1Supercuts+; #pragma link C++ class MCT1SupercutsCalc+; +#pragma link C++ class MCT1FindSupercuts+; +#pragma link C++ class MMinuitInterface+; #pragma link C++ class MFiltercutsCalc+; #endif + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc (revision 2300) +++ /trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc (revision 2300) @@ -0,0 +1,1047 @@ +/* ======================================================================== *\ +! +! * +! * 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 "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; + } + + MHCT1Supercuts *plotsuper = + (MHCT1Supercuts*)plistfcn->FindObject("MHCT1Supercuts"); + if (!plotsuper) + { + gLog << "fcnSupercuts : MHCT1Supercuts object '" << "MHCT1Supercuts" + << "' not found... aborting" << endl; + return; + } + + + //------------------------------------------------------- + // transfer current parameter values to MCT1Supercuts + // + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); + + //gLog << "fcnSupercuts : npar, fNpari, fNparx = " << npar << ", " + // << fNpari << ", " << fNparx << endl; + + TArrayD d(fNparx, par); + super->SetParams(fNparx, d); + + + //------------------------------------------- + // plot alpha with the current cuts + evtloopfcn->Eventloop(); + + TString mh3Name = "AlphaFcn"; + MH3* alpha = (MH3*)plistfcn->FindObject(mh3Name, "MH3"); + + 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 + + Double_t alphasig = 20.0; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + 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; + + Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + 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 + plotsuper->Fill(&findsig); + + + //------------------------ + // get significance + Double_t significance = findsig.GetSignificance(); + f = significance>0.0 ? -significance : 0.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) +{ + fName = name ? name : "MCT1FindSupercuts"; + fTitle = title ? title : "Optimizer of the supercuts"; + + //--------------------------- + // default values + fFilenameTrain = ""; + fFilenameTest = ""; + fFilenameParam = ""; + + fHowManyTrain = 10000; + fHowManyTest = 10000; + + fHadronnessName = ""; + + fMatrixFilter = NULL; + + //--------------------------- + // 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 == "" || 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 ( fabs(howmanygenerated-howmanytrain) > 3.0*sqrt(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 == "" || 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"); + Int_t howmanygenerated = fMatrixTest->GetM().GetNrows(); + if ( fabs(howmanygenerated-howmanytest) > 3.0*sqrt(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); + + Float_t fcorr = (Float_t)read.GetEntries() / + ((Float_t)read.GetEntries()-(Float_t)howmanytrain); + + *fLog << "entries, howmanytrain, fcorr = " << read.GetEntries() + << ", " << howmanytrain << ", " << fcorr << endl; + + //MFEventSelector2 seltest(hreftest); + //seltrain.SetNumMax(howmanytest*fcorr); + MFEventSelector seltest; + seltest.SetName("SelTest"); + seltest.SetNumSelectEvts(howmanytest*fcorr); + + 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); + + Int_t maxevents = -1; + if ( !evtloop.Eventloop(maxevents) ) + return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + fMatrixTrain->Print("SizeCols"); + Int_t generatedtrain = fMatrixTrain->GetM().GetNrows(); + if ( fabs(generatedtrain-howmanytrain) > 3.0*sqrt(howmanytrain) ) + { + *fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events (" + << generatedtrain + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + fMatrixTest->Print("SizeCols"); + Int_t generatedtest = fMatrixTest->GetM().GetNrows(); + if ( fabs(generatedtest-howmanytest) > 3.0*sqrt(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(TString &filetrain, 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(TString &filetrain, 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::GetParams(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 == "") + { + *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; + } + + 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; + super.InitParams(); + + // calculate supercuts hadronness + fCalcHadTrain->SetHadronnessName(fHadronnessName); + + // apply the supercuts + MF scfilter(fHadronnessName+".fHadronness>0.5"); + MContinue supercuts(&scfilter); + + // plot |alpha| + TString mh3Name = "AlphaFcn"; + MBinning binsalpha("Binning"+mh3Name); + binsalpha.SetEdges(54, -12.0, 96.0); + + // |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.SetParList(&parlistfcn); + evtloopfcn.SetName("EvtLoopFCN"); + *fLog << "Event loop for fcnSupercuts has been setup" << endl; + + // address of evtloopfcn is used in CallMinuit() + fObjectFit = &evtloopfcn; + + + //----------------------------------------------------------------------- + // + //---------- 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 + // + + fNpar = 48; + + // get initial values of parameters from MCT1SupercutsCalc + TArrayD vinit(fNpar); + super.GetParams(fNpar, vinit); + + for (UInt_t i=0; iGetM().GetNrows(); + + if ( howmanygenerated <= 0 ) + { + *fLog << "MCT1FindSupercuts::TestParams; test matrix is empty... aborting" + << 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 '" + << fFilenameParam << "' :" << endl; + + TArrayD supercutsPar(72);; + scin.GetParams(72, supercutsPar); + for (Int_t i=0; i<72; i++) + { + *fLog << supercutsPar[i] << ", "; + } + *fLog << endl; + //--------------------------- + + + // ----------------------------------------------------------------- + if ( fHadronnessName == "") + { + *fLog << "MCT1FindSupercuts::TestParams; hadronness name is not defined... aborting" + << endl; + return kFALSE; + } + + + MParList parlist2; + MTaskList tasklist2; + + MCT1Supercuts supercuts; + supercuts.SetParams(72, supercutsPar); + + fCalcHadTest->SetHadronnessName(fHadronnessName); + + + //------------------------------------------- + + MMatrixLoop loop(fMatrixTest); + + // plot alpha before applying the supercuts + 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 + 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", "Alpha", 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 + + Double_t alphasig = 20.0; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 4; + Double_t significance = -99.0; + Bool_t drawpoly = kTRUE; + Bool_t fitgauss = kFALSE; + Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + significance = findsig.GetSignificance(); + 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; +} + + +//=========================================================================== + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h (revision 2300) +++ /trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h (revision 2300) @@ -0,0 +1,118 @@ +#ifndef MARS_MCT1FindSupercuts +#define MARS_MCT1FindSupercuts + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +#include "MFilter.h" +#include "MEvtLoop.h" +#include "MH3.h" +#include "MCT1SupercutsCalc.h" +#include "MGeomCam.h" +#include "MHMatrix.h" + + +class MCT1FindSupercuts : public MParContainer +{ +private: + + TString fFilenameTrain; + TString fFilenameTest; + + Int_t fHowManyTrain; + Int_t fHowManyTest; + + TString fFilenameParam; + + TString fHadronnessName; + + MCT1SupercutsCalc *fCalcHadTrain; + MCT1SupercutsCalc *fCalcHadTest; + + MHMatrix *fMatrixTrain; + MHMatrix *fMatrixTest; + MGeomCam *fCam; + + MEvtLoop *fObjectFit; + + MFilter *fMatrixFilter; + + // to comunicate with MINUIT ----------------- + // attention : dimensions must agree with those in + // MMinuitInterface::CallMinuit() + char fParName [80][100]; + Double_t fVinit[80]; + Double_t fStep[80]; + Double_t fLimlo[80]; + Double_t fLimup[80]; + Int_t fFix[80]; + + UInt_t fNpar; + + TString fMethod; + + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + Int_t fErrMinimize; + //-------------------------------------------- + + +public: + MCT1FindSupercuts(const char *name=NULL, const char *title=NULL); + ~MCT1FindSupercuts(); + + void SetFilenameTraining(const TString &name, const Int_t howmany) + {fFilenameTrain = name; fHowManyTrain = howmany; } + + void SetFilenameTest(const TString &name, const Int_t howmany) + {fFilenameTest = name; fHowManyTest = howmany; } + + void SetFilenameParam(const TString &name) {fFilenameParam = name;} + void SetHadronnessName(const TString &name) {fHadronnessName = name;} + + void SetMatrixFilter(MFilter *filter) {fMatrixFilter = filter;} + + Bool_t DefineTrainMatrix(const TString &name, const Int_t howmany, + MH3 &href); + + Bool_t DefineTestMatrix(const TString &name, const Int_t howmany, + MH3 &href); + + Bool_t DefineTrainTestMatrix(const TString &name, + const Int_t howmanytrain, MH3 &hreftrain, + const Int_t howmanytest, MH3 &hreftest); + + Bool_t ReadMatrix( TString &filetrain, TString &filetest); + Bool_t WriteMatrix(TString &filetrain, TString &filetest); + + Bool_t FindParams(); + Bool_t TestParams(); + + ClassDef(MCT1FindSupercuts, 1) // Class for optimization of the Supercuts +}; + +#endif + + + + + + + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc (revision 2300) +++ /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc (revision 2300) @@ -0,0 +1,346 @@ +/* ======================================================================== *\ +! +! * +! * 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, 08/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MCT1Supercuts // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MCT1Supercuts.h" + +#include +#include + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MCT1Supercuts); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// constructor +// +MCT1Supercuts::MCT1Supercuts(const char *name, const char *title) +{ + fName = name ? name : "MCT1Supercuts"; + fTitle = title ? title : "Container for the supercut parameters"; + + fNpartot = 72; + + fLengthUp.Set(8); + fLengthLo.Set(8); + fWidthUp.Set(8); + fWidthLo.Set(8); + fDistUp.Set(8); + fDistLo.Set(8); + fAsymUp.Set(8); + fAsymLo.Set(8); + fAlphaUp.Set(8); + + // set supercut parameters to their default values + InitParams(); +} + + +// -------------------------------------------------------------------------- +// +// set default values for the supercut parameters +// +void MCT1Supercuts::InitParams() +{ + //--------------------------------- + // supercut parameters + + fLengthUp[0] = 0.315585; + fLengthUp[1] = 0.001455; + fLengthUp[2] = 0.203198; + fLengthUp[3] = 0.005532; + fLengthUp[4] =-0.001670; + fLengthUp[5] =-0.020362; + fLengthUp[6] = 0.007388; + fLengthUp[7] =-0.013463; + + fWidthUp[0] = 0.145412; + fWidthUp[1] =-0.001771; + fWidthUp[2] = 0.054462; + fWidthUp[3] = 0.022280; + fWidthUp[4] =-0.009893; + fWidthUp[5] = 0.056353; + fWidthUp[6] = 0.020711; + fWidthUp[7] =-0.016703; + + fDistUp[0] = 1.787943; + fDistUp[1] = 0.; + fDistUp[2] = 2.942310; + fDistUp[3] = 0.199815; + fDistUp[4] = 0.; + fDistUp[5] = 0.249909; + fDistUp[6] = 0.189697; + fDistUp[7] = 0.; + + fLengthLo[0] = 0.151530; + fLengthLo[1] = 0.028323; + fLengthLo[2] = 0.510707; + fLengthLo[3] = 0.053089; + fLengthLo[4] = 0.013708; + fLengthLo[5] = 2.357993; + fLengthLo[6] = 0.000080; + fLengthLo[7] =-0.007157; + + fWidthLo[0] = 0.089187; + fWidthLo[1] =-0.006430; + fWidthLo[2] = 0.074442; + fWidthLo[3] = 0.003738; + fWidthLo[4] =-0.004256; + fWidthLo[5] =-0.014101; + fWidthLo[6] = 0.006126; + fWidthLo[7] =-0.002849; + + fDistLo[0] = 0.589406; + fDistLo[1] = 0.; + fDistLo[2] =-0.083964; + fDistLo[3] =-0.007975; + fDistLo[4] = 0.; + fDistLo[5] = 0.045374; + fDistLo[6] =-0.001750; + fDistLo[7] = 0.; + + fAsymUp[0] = 0.061267; + fAsymUp[1] = 0.014462; + fAsymUp[2] = 0.014327; + fAsymUp[3] = 0.014540; + fAsymUp[4] = 0.013391; + fAsymUp[5] = 0.012319; + fAsymUp[6] = 0.010444; + fAsymUp[7] = 0.008328; + + fAsymLo[0] =-0.012055; + fAsymLo[1] = 0.009157; + fAsymLo[2] = 0.005441; + fAsymLo[3] = 0.000399; + fAsymLo[4] = 0.001433; + fAsymLo[5] =-0.002050; + fAsymLo[6] =-0.000104; + fAsymLo[7] =-0.001188; + + fAlphaUp[0] = 13.123440; + fAlphaUp[1] = 0.; + fAlphaUp[2] = 0.; + fAlphaUp[3] = 0.; + fAlphaUp[4] = 0.; + fAlphaUp[5] = 0.; + fAlphaUp[6] = 0.; + fAlphaUp[7] = 0.; + //--------------------------------- +} + + +// -------------------------------------------------------------------------- +// +// Set the parameter values from the array 'd' +// +// + +void MCT1Supercuts::SetParams(UInt_t npar, TArrayD &d) +{ + UInt_t ncutpar = fLengthUp.GetSize(); + UInt_t k0 = 0; + + for (UInt_t j=0; j= npar) return; + fLengthUp[j] = d[k]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + fWidthUp[j] = d[k]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + fDistUp[j] = d[k]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + fLengthLo[j] = d[k]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + fWidthLo[j] = d[k]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + fDistLo[j] = d[k]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + fAsymUp[j] = d[k]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + fAsymLo[j] = d[k]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + fAlphaUp[j] = d[k]; + } + k0 += ncutpar; +} + +// -------------------------------------------------------------------------- +// +// Get the parameter values and put them into the array 'd' +// +// + +void MCT1Supercuts::GetParams(UInt_t npar, TArrayD &d) +{ + UInt_t ncutpar = fLengthUp.GetSize(); + UInt_t k0 = 0; + + for (UInt_t j=0; j= npar) return; + d[k] = fLengthUp[j]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + d[k] = fWidthUp[j]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + d[k] = fDistUp[j]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + d[k] = fLengthLo[j]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + d[k] = fWidthLo[j]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + d[k] = fDistLo[j]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + d[k] = fAsymUp[j]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + d[k] = fAsymLo[j]; + } + k0 += ncutpar; + + for (UInt_t j=0; j= npar) return; + d[k] = fAlphaUp[j]; + } +} +//========================================================================== + + + + + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h (revision 2300) +++ /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h (revision 2300) @@ -0,0 +1,75 @@ +#ifndef MARS_MCT1Supercuts +#define MARS_MCT1Supercuts + + +#ifndef ROOT_TArrayD +#include +#endif + +#ifndef ROOT_TString +#include +#endif + +#ifndef MARS_MParContainer +#include +#endif + + + +class MCT1Supercuts : public MParContainer +{ +private: + + UInt_t fNpartot; + + //--------------------------------- + // supercut parameters + + TArrayD fLengthUp; + TArrayD fLengthLo; + TArrayD fWidthUp; + TArrayD fWidthLo; + TArrayD fDistUp; + TArrayD fDistLo; + TArrayD fAsymUp; + TArrayD fAsymLo; + TArrayD fAlphaUp; + //--------------------------------- + + +public: + MCT1Supercuts(const char *name=NULL, const char *title=NULL); + + void InitParams(); + + void SetParams(UInt_t npar, TArrayD &d); + void GetParams(UInt_t npar, TArrayD &d); + + TArrayD &GetLengthUp() { return fLengthUp; } + TArrayD &GetLengthLo() { return fLengthLo; } + TArrayD &GetWidthUp() { return fWidthUp; } + TArrayD &GetWidthLo() { return fWidthLo; } + TArrayD &GetDistUp() { return fDistUp; } + TArrayD &GetDistLo() { return fDistLo; } + TArrayD &GetAsymUp() { return fAsymUp; } + TArrayD &GetAsymLo() { return fAsymLo; } + TArrayD &GetAlphaUp() { return fAlphaUp; } + + ClassDef(MCT1Supercuts, 1) // A container for the Supercut parameters +}; + +#endif + + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2299) +++ /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2300) @@ -33,4 +33,7 @@ #include +#include "TFile.h" +#include "TArrayD.h" + #include "MParList.h" #include "MHillasExt.h" @@ -41,4 +44,5 @@ #include "MHadronness.h" #include "MHMatrix.h" +#include "MCT1Supercuts.h" #include "MLog.h" @@ -49,233 +53,18 @@ using namespace std; -void MCT1SupercutsCalc::InitParams() -{ - fLengthUp.Set(8); - fLengthLo.Set(8); - fWidthUp.Set(8); - fWidthLo.Set(8); - fDistUp.Set(8); - fDistLo.Set(8); - fAsymUp.Set(8); - fAsymLo.Set(8); - fAlphaUp.Set(8); - - //--------------------------------- - // cut parameters - fLengthUp[0] = 0.315585; - fLengthUp[1] = 0.001455; - fLengthUp[2] = 0.203198; - fLengthUp[3] = 0.005532; - fLengthUp[4] =-0.001670; - fLengthUp[5] =-0.020362; - fLengthUp[6] = 0.007388; - fLengthUp[7] =-0.013463; - - fWidthUp[0] = 0.145412; - fWidthUp[1] =-0.001771; - fWidthUp[2] = 0.054462; - fWidthUp[3] = 0.022280; - fWidthUp[4] =-0.009893; - fWidthUp[5] = 0.056353; - fWidthUp[6] = 0.020711; - fWidthUp[7] =-0.016703; - - fDistUp[0] = 1.787943; - fDistUp[1] = 0.; - fDistUp[2] = 2.942310; - fDistUp[3] = 0.199815; - fDistUp[4] = 0.; - fDistUp[5] = 0.249909; - fDistUp[6] = 0.189697; - fDistUp[7] = 0.; - - fLengthLo[0] = 0.151530; - fLengthLo[1] = 0.028323; - fLengthLo[2] = 0.510707; - fLengthLo[3] = 0.053089; - fLengthLo[4] = 0.013708; - fLengthLo[5] = 2.357993; - fLengthLo[6] = 0.000080; - fLengthLo[7] =-0.007157; - - fWidthLo[0] = 0.089187; - fWidthLo[1] =-0.006430; - fWidthLo[2] = 0.074442; - fWidthLo[3] = 0.003738; - fWidthLo[4] =-0.004256; - fWidthLo[5] =-0.014101; - fWidthLo[6] = 0.006126; - fWidthLo[7] =-0.002849; - - fDistLo[0] = 0.589406; - fDistLo[1] = 0.; - fDistLo[2] =-0.083964; - fDistLo[3] =-0.007975; - fDistLo[4] = 0.; - fDistLo[5] = 0.045374; - fDistLo[6] =-0.001750; - fDistLo[7] = 0.; - - fAsymUp[0] = 0.061267; - fAsymUp[1] = 0.014462; - fAsymUp[2] = 0.014327; - fAsymUp[3] = 0.014540; - fAsymUp[4] = 0.013391; - fAsymUp[5] = 0.012319; - fAsymUp[6] = 0.010444; - fAsymUp[7] = 0.008328; - - fAsymLo[0] =-0.012055; - fAsymLo[1] = 0.009157; - fAsymLo[2] = 0.005441; - fAsymLo[3] = 0.000399; - fAsymLo[4] = 0.001433; - fAsymLo[5] =-0.002050; - fAsymLo[6] =-0.000104; - fAsymLo[7] =-0.001188; - - fAlphaUp[0] = 13.123440; - fAlphaUp[1] = 0.; - fAlphaUp[2] = 0.; - fAlphaUp[3] = 0.; - fAlphaUp[4] = 0.; - fAlphaUp[5] = 0.; - fAlphaUp[6] = 0.; - fAlphaUp[7] = 0.; - //--------------------------------- -} - -// -------------------------------------------------------------------------- -// -// Set the parameter values from vector 'par' -// -// Attention : it is assumed that there are (9*ncutpar) values -// -void MCT1SupercutsCalc::SetParams(Double_t *par) -{ - UInt_t ncutpar = fLengthUp.GetSize(); - UInt_t k0 = 0; - - TArrayD lup(ncutpar, par + k0); - fLengthUp = lup; - k0 += ncutpar; - - TArrayD wup(ncutpar, par + k0); - fWidthUp = wup; - k0 += ncutpar; - - TArrayD dup(ncutpar, par + k0); - fDistUp = dup; - k0 += ncutpar; - - TArrayD llo(ncutpar, par + k0); - fLengthLo = llo; - k0 += ncutpar; - - TArrayD wlo(ncutpar, par + k0); - fWidthLo = wlo; - k0 += ncutpar; - - TArrayD dlo(ncutpar, par + k0); - fDistLo = dlo; - k0 += ncutpar; - - TArrayD aup(ncutpar, par + k0); - fAsymUp = aup; - k0 += ncutpar; - - TArrayD alo(ncutpar, par + k0); - fAsymLo = alo; - k0 += ncutpar; - - TArrayD alphaup(ncutpar, par + k0); - fAlphaUp = alphaup; -} - -// -------------------------------------------------------------------------- -// -// Get the parameter values -// -// Attention : it is assumed that there are (9*ncutpar) values -// -void MCT1SupercutsCalc::GetParams(Double_t *par) -{ - UInt_t ncutpar = fLengthUp.GetSize(); - UInt_t k0 = 0; - - for (UInt_t j=0; jFindObject(fSuperName, "MCT1Supercuts"); + if (!fSuper) + { + *fLog << err << fSuperName << " [MCT1Supercuts] not found... aborting." << endl; + return kFALSE; + } + + if (fMatrix) return kTRUE; + //----------------------------------------------------------- fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); if (!fHil) @@ -325,4 +123,5 @@ } + return kTRUE; } @@ -332,9 +131,5 @@ // Calculation of upper and lower limits // -Double_t MCT1SupercutsCalc::CtsMCut( -#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00) -const -#endif - TArrayD &a, Double_t ls, Double_t ct, +Double_t MCT1SupercutsCalc::CtsMCut(TArrayD &a, Double_t ls, Double_t ct, Double_t ls2, Double_t dd2) { @@ -398,5 +193,5 @@ fMap[5] = fMatrix->AddColumn("MHillas.fMeanY"); fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist"); - fMap[7] = fMatrix->AddColumn("MHillasSrc.fAlpha"); + fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)"); } @@ -437,16 +232,16 @@ const Double_t width = width0 * fMm2Deg; - if (newdist < 1.05 && - newdist < CtsMCut (fDistUp, dmls, dmcza, dmls2, dd2) && - newdist > CtsMCut (fDistLo, dmls, dmcza, dmls2, dd2) && - dist < 1.05 && - length < CtsMCut (fLengthUp, dmls, dmcza, dmls2, dd2) && - length > CtsMCut (fLengthLo, dmls, dmcza, dmls2, dd2) && - width < CtsMCut (fWidthUp, dmls, dmcza, dmls2, dd2) && - width > CtsMCut (fWidthLo, dmls, dmcza, dmls2, dd2) && - //asym < CtsMCut (asymup, dmls, dmcza, dmls2, dd2) && - //asym > CtsMCut (asymlow, dmls, dmcza, dmls2, dd2) && - dist < CtsMCut (fDistUp, dmls, dmcza, dmls2, dd2) && - dist > CtsMCut (fDistLo, dmls, dmcza, dmls2, dd2) ) + if (newdist < 1.05 && + newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && + newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && + dist < 1.05 && + length < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) && + length > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) && + width < CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) && + width > CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) && + //asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) && + //asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) && + dist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && + dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) ) fHadronness->SetHadronness(0.25); else @@ -458,2 +253,6 @@ } //========================================================================== + + + + Index: /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2299) +++ /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2300) @@ -18,4 +18,5 @@ class MHadronness; class MHMatrix; +class MCT1Supercuts; @@ -23,12 +24,14 @@ { private: - MHillas *fHil; - MHillasSrc *fHilSrc; - MMcEvt *fMcEvt; - MHadronness *fHadronness; //! output container for hadronness + MHillas *fHil; + MHillasSrc *fHilSrc; + MMcEvt *fMcEvt; + MHadronness *fHadronness; //! output container for hadronness + MCT1Supercuts *fSuper; // container for supercut parameters TString fHadronnessName; // name of container to store hadronness TString fHilName; TString fHilSrcName; + TString fSuperName; // name of container for supercut parameters Double_t fMm2Deg; @@ -37,30 +40,10 @@ MHMatrix *fMatrix; - //--------------------------------- - // cut parameters - TArrayD fLengthUp; - TArrayD fLengthLo; - TArrayD fWidthUp; - TArrayD fWidthLo; - TArrayD fDistUp; - TArrayD fDistLo; - TArrayD fAsymUp; - TArrayD fAsymLo; - TArrayD fAlphaUp; - //--------------------------------- - - void InitParams(); - Int_t PreProcess(MParList *pList); Int_t Process(); - void Set(TArrayD &a, const TArrayD &b) { if (a.GetSize()==b.GetSize()) a=b; } Double_t GetVal(Int_t i) const; - Double_t CtsMCut( -#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00) -const -#endif - TArrayD &a, Double_t ls, Double_t ct, + Double_t CtsMCut(TArrayD &a, Double_t ls, Double_t ct, Double_t ls2, Double_t dd2); @@ -70,7 +53,4 @@ const char *name=NULL, const char *title=NULL); - void SetParams(Double_t *par); - void GetParams(Double_t *par); - void SetHadronnessName(const TString name) { fHadronnessName = name; } TString GetHadronnessName() const { return fHadronnessName; } @@ -79,23 +59,4 @@ void StopMapping() { InitMapping(NULL); } - void SetLengthUp(const TArrayD &d) { Set(fLengthUp, d); } - void SetLengthLo(const TArrayD &d) { Set(fLengthLo, d); } - void SetWidthUp (const TArrayD &d) { Set(fWidthUp, d); } - void SetWidthLo (const TArrayD &d) { Set(fWidthLo, d); } - void SetDistUp (const TArrayD &d) { Set(fDistUp, d); } - void SetDistLo (const TArrayD &d) { Set(fDistLo, d); } - void SetAsymUp (const TArrayD &d) { Set(fAsymUp, d); } - void SetAsymLo (const TArrayD &d) { Set(fAsymLo, d); } - void SetAlphaUp (const TArrayD &d) { Set(fAlphaUp, d); } - - const TArrayD &GetLengthUp() const { return fLengthUp; } - const TArrayD &GetLengthLo() const { return fLengthLo; } - const TArrayD &GetWidthUp() const { return fWidthUp; } - const TArrayD &GetWithLo() const { return fWidthLo; } - const TArrayD &GetDistUp() const { return fDistUp; } - const TArrayD &GetDistLo() const { return fDistLo; } - const TArrayD &GetAsymUp() const { return fAsymUp; } - const TArrayD &GetAsymLo() const { return fAsymLo; } - const TArrayD &GetAlphaUp() const { return fAlphaUp; } ClassDef(MCT1SupercutsCalc, 0) // A class to evaluate the Supercuts @@ -111,2 +72,8 @@ + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc (revision 2300) +++ /trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc (revision 2300) @@ -0,0 +1,325 @@ +/* ======================================================================== *\ +! +! * +! * 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 7/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MMinuitInterface // +// // +// Class for interfacing with Minuit // +// // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MMinuitInterface.h" + +#include // fabs + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" +#include "MParContainer.h" + + +ClassImp(MMinuitInterface); + +using namespace std; + + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MMinuitInterface::MMinuitInterface(const char *name, const char *title) +{ + fName = name ? name : "MMinuitInterface"; + fTitle = title ? title : "Interface for Minuit"; + +} + +// -------------------------------------------------------------------------- +// +// Default destructor. +// +MMinuitInterface::~MMinuitInterface() +{ +} + + +// ----------------------------------------------------------------------- +// +// Interface to MINUIT +// +// +Bool_t MMinuitInterface::CallMinuit( + void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t), + UInt_t npar, char name[80][100], + Double_t vinit[80], Double_t step[80], + Double_t limlo[80], Double_t limup[80], Int_t fix[80], + TObject *objectfit , TString method, Bool_t nulloutput) +{ + TString str(method); + + // + // Be carefull: This is not thread safe + // + UInt_t maxpar = 100; + + if (npar > maxpar) + { + *fLog << "CallMinuit : too many parameters, npar = " << npar + << ", maxpar = " << maxpar << endl; + return kFALSE; + } + + + + //.............................................. + // Set the maximum number of parameters + TMinuit &minuit = *(new TMinuit(maxpar)); + minuit.SetName("MinuitSupercuts"); + //TMinuit minuit(maxpar); + //gMinuit = new TMinuit(maxpar); + //TMinuit &minuit = *gMinuit; + + + //.............................................. + // Set the print level + // -1 no output except SHOW comands + // 0 minimum output + // 1 normal output (default) + // 2 additional ouput giving intermediate results + // 3 maximum output, showing progress of minimizations + // + //Int_t printLevel = -1; + Int_t printLevel = 1; + minuit.SetPrintLevel(printLevel); + + //.............................................. + // Printout for warnings + // SET WAR print warnings + // SET NOW suppress warnings + Int_t errWarn; + Double_t tmpwar = 0; + minuit.mnexcm("SET NOW", &tmpwar, 0, errWarn); + //minuit.mnexcm("SET WAR", &tmpwar, 0, errWarn); + + //.............................................. + // Set the address of the minimization function + minuit.SetFCN(fcn); + + + //.............................................. + // Store address of object to be used in fcn + minuit.SetObjectFit(objectfit ); + + + //.............................................. + // Set starting values and step sizes for parameters + for (UInt_t i=0; i 0) + { + Int_t parNo = i; + minuit.FixParameter(parNo); + step[parNo] = 0.0; + } + } + + //.............................................. + //NumPars = minuit.GetNumPars(); + //*fLog << "CallMinuit : number of free parameters = " + // << NumPars << endl; + + //.............................................. + // Set maximum number of iterations (default = 500) + //Int_t maxiter = 100000; + Int_t maxiter = 10; + minuit.SetMaxIterations(maxiter); + + //.............................................. + // minimization by the method of Migrad + if (str.Contains("Migrad", TString::kIgnoreCase)) + { + //*fLog << "call MIGRAD" << endl; + if (nulloutput) + fLog->SetNullOutput(kTRUE); + Double_t tmp = 0; + minuit.mnexcm("MIGRAD", &tmp, 0, fErrMinimize); + if (nulloutput) + fLog->SetNullOutput(kFALSE); + //*fLog << "return from MIGRAD" << endl; + } + + //.............................................. + // same minimization as by Migrad + // but switches to the SIMPLEX method if MIGRAD fails to converge + if (str.Contains("Minimize", TString::kIgnoreCase)) + { + *fLog << "call MINIMIZE" << endl; + Double_t tmp = 0; + minuit.mnexcm("MINIMIZE", &tmp, 0, fErrMinimize); + *fLog << "return from MINIMIZE" << endl; + } + + //.............................................. + // minimization by the SIMPLEX method + if (str.Contains("Simplex", TString::kIgnoreCase)) + { + *fLog << "call SIMPLEX" << endl; + if (nulloutput) + fLog->SetNullOutput(kTRUE); + Double_t tmp = 0; + minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize); + if (nulloutput) + fLog->SetNullOutput(kFALSE); + //*fLog << "return from SIMPLEX" << endl; + } + + //.............................................. + // check quality of minimization + // istat = 0 covariance matrix not calculated + // 1 diagonal approximation only (not accurate) + // 2 full matrix, but forced positive-definite + // 3 full accurate covariance matrix + // (indication of normal convergence) + minuit.mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); + + //if (fErrMinimize != 0 || fIstat < 3) + if (fErrMinimize != 0) + { + *fLog << "CallMinuit : Minimization failed" << endl; + *fLog << " fMin = " << fMin << ", fEdm = " << fEdm + << ", fErrdef = " << fErrdef << ", fIstat = " << fIstat + << ", fErrMinimize = " << fErrMinimize << endl; + return kFALSE; + } + + //*fLog << "CallMinuit : Minimization was successful" << endl; + //*fLog << " fMin = " << fMin << ", fEdm = " << fEdm + // << ", fErrdef = " << fErrdef << ", fIstat = " << fIstat + // << ", fErrMinimize = " << fErrMinimize << endl; + + + //.............................................. + // minimization by the method of Migrad + if (str.Contains("Hesse", TString::kIgnoreCase)) + { + //*fLog << "call HESSE" << endl; + Double_t tmp = 0; + minuit.mnexcm("HESSE", &tmp, 0, fErrMinimize); + //*fLog << "return from HESSE" << endl; + } + + //.............................................. + // Minos error analysis + if (str.Contains("Minos", TString::kIgnoreCase)) + { + //*fLog << "call MINOS" << endl; + Double_t tmp = 0; + minuit.mnexcm("MINOS", &tmp, 0, fErrMinimize); + //*fLog << "return from MINOS" << endl; + } + + //.............................................. + // Print current status of minimization + // if nkode = 0 only function value + // 1 parameter values, errors, limits + // 2 values, errors, step sizes, internal values + // 3 values, errors, step sizes, 1st derivatives + // 4 values, paraboloc errors, MINOS errors + + //Int_t nkode = 4; + //minuit.mnprin(nkode, fmin); + + //.............................................. + // call fcn with IFLAG = 3 (final calculation : calculate p(chi2)) + // iflag = 1 initial calculations only + // 2 calculate 1st derivatives and function + // 3 calculate function only + // 4 calculate function + final calculations + const char *command = "CALL"; + Double_t iflag = 3; + Int_t errfcn3; + minuit.mnexcm(command, &iflag, 1, errfcn3); + + //*fLog << "Exit CallMinuit" << endl; + + return kTRUE; +} + + +//=========================================================================== + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MMinuitInterface.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MMinuitInterface.h (revision 2300) +++ /trunk/MagicSoft/Mars/manalysis/MMinuitInterface.h (revision 2300) @@ -0,0 +1,68 @@ +#ifndef MARS_MMinuitInterface +#define MARS_MMinuitInterface + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + + + +class MMinuitInterface : public MTask +{ +private: + + UInt_t fNpar; + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + Int_t fErrMinimize; + + +public: + MMinuitInterface(const char *name=NULL, const char *title=NULL); + ~MMinuitInterface(); + +Bool_t MMinuitInterface::CallMinuit( + void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t), + UInt_t npar, char name[80][100], + Double_t vinit[80], Double_t step[80], + Double_t limlo[80], Double_t limup[80], Int_t fix[80], + TObject *fObjectFit, TString method, Bool_t nulloutput); + + ClassDef(MMinuitInterface, 1) // Class for interfacing with Minuit +}; + +#endif + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- /trunk/MagicSoft/Mars/manalysis/Makefile (revision 2299) +++ /trunk/MagicSoft/Mars/manalysis/Makefile (revision 2300) @@ -24,5 +24,5 @@ INCLUDES = -I. -I../mbase -I../mmc -I../mraw -I../mgeom -I../mfilter \ -I../mdata -I../mhist -I../mgui -I../mimage -I../mhistmc \ - -I../mfileio + -I../mfileio -I../mmain @@ -60,5 +60,8 @@ MMcTriggerLvl2.cc \ MMcTriggerLvl2Calc.cc \ + MCT1Supercuts.cc \ MCT1SupercutsCalc.cc \ + MCT1FindSupercuts.cc \ + MMinuitInterface.cc \ MFiltercutsCalc.cc # MCT1PadONOFF.cc \ Index: /trunk/MagicSoft/Mars/mhist/HistLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 2299) +++ /trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 2300) @@ -39,4 +39,7 @@ #pragma link C++ class MHSigmaTheta+; #pragma link C++ class MHSigmaPixel+; +#pragma link C++ class MHOnSubtraction+; +#pragma link C++ class MHFindSignificance+; +#pragma link C++ class MHCT1Supercuts+; #pragma link C++ class MHCompProb+; @@ -47,2 +50,11 @@ #endif + + + + + + + + + Index: /trunk/MagicSoft/Mars/mhist/MH3.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MH3.cc (revision 2299) +++ /trunk/MagicSoft/Mars/mhist/MH3.cc (revision 2300) @@ -241,4 +241,6 @@ Bool_t MH3::SetupFill(const MParList *plist) { + // reset histogram (necessary if the same eventloop is run more than once) + fHist->Reset(); TString bname("Binning"); Index: /trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc (revision 2300) +++ /trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc (revision 2300) @@ -0,0 +1,246 @@ +/* ======================================================================== *\ +! +! * +! * 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 2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////// +// +// MHCT1Supercuts +// +// This class contains histograms for the supercuts +// +/////////////////////////////////////////////////////////////////////// +#include "MHCT1Supercuts.h" + +#include + +#include +#include +#include +#include + +#include "MParList.h" +#include "MHFindSignificance.h" + +#include "MLog.h" +#include "MLogManip.h" + + + +ClassImp(MHCT1Supercuts); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Setup four histograms for Alpha, and Dist +// +MHCT1Supercuts::MHCT1Supercuts(const char *name, const char *title) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHCT1Supercuts"; + fTitle = title ? title : "Container for supercuts"; + + + fDegree = new TH1F("Degree", "Degree of polynomial", 5, -0.5, 4.5); + fProb = new TH1F("Prob", "chi2 probability", 40, 0, 1.0); + fNdf = new TH1F("NDF", "NDF of polynomial fit", 60, -0.5, 59.5); + fGamma = new TH1F("gamma", "gamma", 40, 0.0, 8.0); + fNexNon = new TH1F("NexNon", "Nex / Non", 50, 0.0, 1.0); + fSigLiMa= new TH1F("Significance", "significance of gamma signal", + 60, 0.0, 120.0); + fSigtoBackg= new TH2F("SigtoBackg", "Significance vs signal/backg ratio", + 50, 0.0, 10.0, 60, 0.0, 120.0); + fSigDegree = new TH2F("SigDegree", "Significance vs Degree of polynomial", + 5, -0.5, 4.5, 60, 0.0, 120.0); + fSigNbins = new TH2F("SigNbins", "Significance vs number of bins", + 40, -0.5, 79.5, 60, 0.0, 120.0); + + fDegree->SetDirectory(NULL); + fProb->SetDirectory(NULL); + fNdf->SetDirectory(NULL); + fGamma->SetDirectory(NULL); + fNexNon->SetDirectory(NULL); + fSigLiMa->SetDirectory(NULL); + fSigtoBackg->SetDirectory(NULL); + fSigDegree->SetDirectory(NULL); + fSigNbins->SetDirectory(NULL); + + fDegree->SetXTitle("order of polynomial"); + fProb->SetXTitle("chi2 probability of polynomial fit"); + fNdf->SetXTitle("NDF of polynomial fit"); + fGamma->SetXTitle("gamma"); + fNexNon->SetXTitle("Nex / Non"); + fSigLiMa->SetXTitle("significance"); + + fSigtoBackg->SetXTitle("signa./background ratio"); + fSigtoBackg->SetYTitle("significance"); + + fSigDegree->SetXTitle("order of polynomial"); + fSigDegree->SetYTitle("significance"); + + fSigNbins->SetXTitle("number of bins"); + fSigNbins->SetYTitle("significance"); + + fDegree->SetYTitle("Counts"); + fProb->SetYTitle("Counts"); + fNdf->SetYTitle("Counts"); + fGamma->SetYTitle("Counts"); + fNexNon->SetYTitle("Counts"); + fSigLiMa->SetYTitle("Counts"); + + fSigtoBackg->SetZTitle("Counts"); + fSigDegree->SetZTitle("Counts"); + fSigNbins->SetZTitle("Counts"); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms +// +MHCT1Supercuts::~MHCT1Supercuts() +{ + delete fDegree; + delete fProb; + delete fNdf; + delete fGamma; + delete fNexNon; + delete fSigLiMa; + + delete fSigtoBackg; + delete fSigDegree; + delete fSigNbins; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers (if there are any) +// + +Bool_t MHCT1Supercuts::SetupFill(const MParList *plist) +{ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms from the 'MHFindSignificance' container +// + +Bool_t MHCT1Supercuts::Fill(const MParContainer *par, const Stat_t w) +{ + if (!par) + { + *fLog << err << "MHCT1Supercuts::Fill: Pointer (!=NULL) expected." << endl; + return kFALSE; + } + + MHFindSignificance &h = *(MHFindSignificance*)par; + + fDegree ->Fill(h.GetDegree( ), w); + fProb ->Fill(h.GetProb(), w); + fNdf ->Fill(h.GetNdf(), w); + fGamma ->Fill(h.GetGamma(), w); + + Double_t ratio = h.GetNon()>0.0 ? h.GetNex()/h.GetNon() : 0.0; + fNexNon ->Fill(ratio, w); + + fSigLiMa ->Fill(h.GetSigLiMa(), w); + + Double_t sigtobackg = h.GetNbg()!=0.0 ? h.GetNex() / h.GetNbg() : 0.0; + fSigtoBackg->Fill(sigtobackg, h.GetSigLiMa(), w); + + fSigDegree ->Fill(h.GetDegree(), h.GetSigLiMa(), w); + fSigNbins ->Fill(h.GetMbins(), h.GetSigLiMa(), w); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Creates a new canvas and draws the two histograms into it. +// Be careful: The histograms belongs to this object and won't get deleted +// together with the canvas. +// +void MHCT1Supercuts::Draw(Option_t *) +{ + //TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + //pad->SetBorderMode(0); + //AppendPad(""); + + TCanvas *pad = new TCanvas("Supercuts", "Supercut plots", 900, 900); + gROOT->SetSelectedPad(NULL); + + pad->Divide(3, 3); + + pad->cd(1); + gPad->SetBorderMode(0); + fDegree->Draw(); + + pad->cd(2); + gPad->SetBorderMode(0); + fProb->Draw(); + + pad->cd(3); + gPad->SetBorderMode(0); + fNdf->Draw(); + + pad->cd(4); + gPad->SetBorderMode(0); + fGamma->Draw(); + + pad->cd(5); + gPad->SetBorderMode(0); + fNexNon->Draw(); + + pad->cd(6); + gPad->SetBorderMode(0); + fSigLiMa->Draw(); + + pad->cd(7); + gPad->SetBorderMode(0); + fSigtoBackg->Draw(); + + pad->cd(8); + gPad->SetBorderMode(0); + fSigDegree->Draw(); + + pad->cd(9); + gPad->SetBorderMode(0); + fSigNbins->Draw(); + + pad->Modified(); + pad->Update(); +} +//========================================================================= + + + + + + + + + Index: /trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.h (revision 2300) +++ /trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.h (revision 2300) @@ -0,0 +1,44 @@ +#ifndef MARS_MHCT1Supercuts +#define MARS_MHCT1Supercuts + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1F; +class TH2F; + +class MHCT1Supercuts : public MH +{ +private: + TH1F *fDegree; // order of polynomial for background fit + TH1F *fProb; // chi2 probability of polynomial fit + TH1F *fNdf; // NDF of polynomial fit + TH1F *fGamma; // Nbg = gamma * Noff + TH1F *fNexNon; // no.of excess events / no.of events in signal region + TH1F *fSigLiMa; // significance of gamma signal + + TH2F *fSigtoBackg; // significance vs signal to background ratio (Nex/Nbg) + TH2F *fSigDegree; // significance vs order of polynomial + TH2F *fSigNbins; // significance vs number of bins + + +public: + MHCT1Supercuts(const char *name=NULL, const char *title=NULL); + ~MHCT1Supercuts(); + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + + void Draw(Option_t *opt=NULL); + + + ClassDef(MHCT1Supercuts, 1) // Container which holds histograms for the supercuts +}; + +#endif + + + + + Index: /trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc (revision 2300) +++ /trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc (revision 2300) @@ -0,0 +1,2117 @@ +/* ======================================================================== *\ +! +! * +! * 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, July 2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHFindSignificance +// +// determines the significance of a gamma signal in an |alpha| plot +// +// +// Input : TH1 histogram of |alpha| : with 0 < |alpha| < 90 degrees +// alphamin, alphamax : defining the background region +// alphasig : defining the signal region for which +// the significance is calculated +// degree : the degree of the polynomial to be fitted to the background +// ( a0 + a1*x + a2*x**2 + a3*x**3 + ...) +// +// Output : +// +// - polynomial which describes the background in the background region +// - the number of events in the signal region (Non) +// the number of background events in the signal region (Nbg) +// - the number of excess events in the signal region (Nex = Non - Nbg) +// - thew effective number of background events (Noff), and gamma : +// Nbg = gamma * Noff +// - the significance of the gamma signal according to Li & Ma +// +// +// call member function 'FindSigma' +// to fit the background and to determine the significance +// +// call the member function 'SigmaVsAlpha' +// to determine the significance as a function of alphasig +// +///////////////////////////////////////////////////////////////////////////// +#include "MHFindSignificance.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" +#include "MMinuitInterface.h" + + +ClassImp(MHFindSignificance); + +using namespace std; + +const TString MHFindSignificance::gsDefName = "MHFindSignificance"; +const TString MHFindSignificance::gsDefTitle = "Find Significance in alpha plot"; + + + +// -------------------------------------------------------------------------- +// +// fcnpoly +// +// calculates the chi2 for the fit of the polynomial function 'poly' +// to the histogram 'fhist' +// +// it is called by CallMinuit() (which is called in FitPolynomial()) +// +// bins of fhist with huge errors are ignored in the calculation of the chi2 +// (the huge errors were set in 'FitPolynomial()') +// + +static void fcnpoly(Int_t &npar, Double_t *gin, Double_t &f, + Double_t *par, Int_t iflag) +{ + TH1 *fhist = (TH1*)gMinuit->GetObjectFit(); + TF1 *fpoly = fhist->GetFunction("Poly"); + + + //------------------------------------------- + + Double_t chi2 = 0.0; + + Int_t nbins = fhist->GetNbinsX(); + Int_t mbins = 0; + for (Int_t i=1; i<=nbins; i++) + { + Double_t content = fhist->GetBinContent(i); + Double_t error = fhist->GetBinError(i); + Double_t center = fhist->GetBinCenter(i); + + //----------------------------- + // ignore unwanted points + if (error > 1.e19) + continue; + + if (content <= 0.0) + { + gLog << "fcnpoly : bin with zero content; i, content, error = " + << i << ", " << content << ", " << error << endl; + continue; + } + + if (error <= 0.0) + { + gLog << "fcnpoly : bin with zero error; i, content, error = " + << i << ", " << content << ", " << error << endl; + continue; + } + + //----------------------------- + mbins++; + + Double_t fu; + fu = fpoly->EvalPar(¢er, par); + + // the fitted function must not be negative + if (fu <= 0.0) + { + chi2 = 1.e10; + break; + } + + Double_t temp = (content - fu) / error; + chi2 += temp*temp; + } + + //------------------------------------------- + + f = chi2; + + //------------------------------------------- + // final calculations + //if (iflag == 3) + //{ + //} + + //------------------------------------------------------------- +} + + +// -------------------------------------------------------------------------- +// +// fcnpolygauss +// +// calculates the chi2 for the fit of the (polynomial+Gauss) function +// 'PolyGauss' to the histogram 'fhist' +// +// it is called by CallMinuit() (which is called in FitGaussPoly()) +// +// bins of fhist with huge errors are ignored in the calculation of the chi2 +// (the huge errors were set in 'FitGaussPoly()') +// + +static void fcnpolygauss(Int_t &npar, Double_t *gin, Double_t &f, + Double_t *par, Int_t iflag) +{ + TH1 *fhist = (TH1*)gMinuit->GetObjectFit(); + TF1 *fpolygauss = fhist->GetFunction("PolyGauss"); + + + //------------------------------------------- + + Double_t chi2 = 0.0; + + Int_t nbins = fhist->GetNbinsX(); + Int_t mbins = 0; + for (Int_t i=1; i<=nbins; i++) + { + Double_t content = fhist->GetBinContent(i); + Double_t error = fhist->GetBinError(i); + Double_t center = fhist->GetBinCenter(i); + + //----------------------------- + // ignore unwanted points + if (error > 1.e19) + continue; + + if (content <= 0.0) + { + gLog << "fcnpolygauss : bin with zero content; i, content, error = " + << i << ", " << content << ", " << error << endl; + continue; + } + + if (error <= 0.0) + { + gLog << "fcnpolygauss : bin with zero error; i, content, error = " + << i << ", " << content << ", " << error << endl; + continue; + } + + //----------------------------- + mbins++; + + Double_t fu; + fu = fpolygauss->EvalPar(¢er, par); + + // the fitted function must not be negative + if (fu <= 0.0) + { + chi2 = 1.e10; + break; + } + + Double_t temp = (content - fu) / error; + chi2 += temp*temp; + } + + //------------------------------------------- + + f = chi2; + + //------------------------------------------- + // final calculations + //if (iflag == 3) + //{ + //} + + //------------------------------------------------------------- +} + + + +// -------------------------------------------------------------------------- +// +// Constructor +// +MHFindSignificance::MHFindSignificance(const char *name, const char *title) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + fSigVsAlpha = NULL; + + fPoly = NULL; + fGPoly = NULL; + fGBackg = NULL; + + fHist = NULL; + fHistOrig = NULL; + + // allow rebinning of the alpha plot + fRebin = kTRUE; + + // allow reducing the degree of the polynomial + fReduceDegree = kTRUE; + + fCanvas = NULL; +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +// =====> it is not clear why one obtains sometimes a segmentation violation +// when the destructor is active <======================= +// +// therefore the 'return'statement +// + +MHFindSignificance::~MHFindSignificance() +{ + return; + + *fLog << "destructor of MHFindSignificance is called" << endl; + + //delete fHist; + + delete fSigVsAlpha; + delete fPoly; + delete fGPoly; + delete fGBackg; + delete fCanvas; +} + +// -------------------------------------------------------------------------- +// +// Set flag fRebin +// +// if flag is kTRUE rebinning of the alpha plot is allowed +// +// +void MHFindSignificance::SetRebin(Bool_t b) +{ + fRebin = b; + + *fLog << "MHFindSignificance::SetRebin; flag fRebin set to " + << (b? "kTRUE" : "kFALSE") << endl; +} + +// -------------------------------------------------------------------------- +// +// Set flag fReduceDegree +// +// if flag is kTRUE reducing of the degree of the polynomial is allowed +// +// +void MHFindSignificance::SetReduceDegree(Bool_t b) +{ + fReduceDegree = b; + + *fLog << "MHFindSignificance::SetReduceDegree; flag fReduceDegree set to " + << (b? "kTRUE" : "kFALSE") << endl; +} + +// -------------------------------------------------------------------------- +// +// FindSigma +// +// calls FitPolynomial to fit the background in the background region +// calls DetExcess to determine the number of excess events +// using an extrapolation of the polynomial +// into the signal region +// calls SigmaLiMa to determine the significance of the gamma signal +// in the range |alpha| < alphasig +// calls FitGaussPoly to fit a (polynomial+Gauss) function in the +// whole |alpha| region +// +// +Bool_t MHFindSignificance::FindSigma(TH1 *fhist, Double_t alphamin, + Double_t alphamax, Int_t degree, Double_t alphasig, + Bool_t drawpoly, Bool_t fitgauss, Bool_t print) +{ + *fLog << "MHFindSignificance::FindSigma;" << endl; + + fHistOrig = fhist; + + fHist = (TH1*)fHistOrig->Clone("|alpha| plot"); + if ( !fHist ) + { + *fLog << "MHFindSignificance::FindSigma; Clone of histogram could not be generated" + << endl; + return kFALSE; + } + + fHist->Sumw2(); + fHist->SetNameTitle("Alpha", "alpha plot"); + fHist->SetXTitle("|alpha| [\\circ]"); + fHist->SetYTitle("Counts"); + + fAlphamin = alphamin; + fAlphamax = alphamax; + fAlphammm = (alphamin+alphamax)/2.0; + fDegree = degree; + fAlphasig = alphasig; + + fDraw = drawpoly; + fFitGauss = fitgauss; + + + //-------------------------------------------- + // fit a polynomial in the background region + + //*fLog << "MHFindSignificance::FindSigma; calling FitPolynomial()" << endl; + if ( !FitPolynomial() ) + { + *fLog << "MHFindSignificance::FindSigma; FitPolynomial failed" + << endl; + return kFALSE; + } + + + //-------------------------------------------- + // calculate the number of excess events in the signal region + + //*fLog << "MHFindSignificance::FindSigma; calling DetExcess()" << endl; + if ( !DetExcess() ) + { + *fLog << "MHFindSignificance::FindSigma; DetExcess failed" + << endl; + return kFALSE; + } + + + //-------------------------------------------- + // calculate the significance of the excess + + //*fLog << "MHFindSignificance::FindSigma; calling SigmaLiMa()" << endl; + Double_t siglima = 0.0; + if ( !SigmaLiMa(fNon, fNoff, fGamma, &siglima) ) + { + *fLog << "MHFindSignificance::FindSigma; SigmaLiMa failed" + << endl; + return kFALSE; + } + fSigLiMa = siglima; + + //-------------------------------------------- + // calculate the error of the number of excess events + + fdNex = fNex / fSigLiMa; + + + //-------------------------------------------- + + //*fLog << "MHFindSignificance::FindSigma; calling PrintPoly()" << endl; + if (print) + PrintPoly(); + + + //-------------------------------------------- + // fit a (polynomial + Gauss) function + + if (fFitGauss) + { + //-------------------------------------------------- + // delete objects from this fit + // in order to have independent starting conditions for the next fit + + delete gMinuit; + gMinuit = NULL; + //-------------------------------------------------- + + //*fLog << "MHFindSignificance::FindSigma; calling FitGaussPoly()" + // << endl; + if ( !FitGaussPoly() ) + { + *fLog << "MHFindSignificance::FindSigma; FitGaussPoly failed" + << endl; + return kFALSE; + } + + if (print) + { + //*fLog << "MHFindSignificance::FindSigma; calling PrintPolyGauss()" + // << endl; + PrintPolyGauss(); + } + } + + //-------------------------------------------------- + // draw the histogram if requested + + if (fDraw) + { + //*fLog << "MHFindSignificance::FindSigma; calling DrawFit()" << endl; + if ( !DrawFit() ) + { + *fLog << "MHFindSignificance::FindSigma; DrawFit failed" + << endl; + return kFALSE; + } + } + + + //-------------------------------------------------- + // delete objects from this fit + // in order to have independent starting conditions for the next fit + + delete gMinuit; + gMinuit = NULL; + //-------------------------------------------------- + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// SigmaVsAlpha (like FindSigma. However, alphasig is scanned and +// the significance is plotted versus alphasig) +// +// calls FitPolynomial to fit the background in the background region +// +// scan alphasig; for a given alphasig : +// calls DetExcess to determine the number of excess events +// calls SigmaLiMa to determine the significance of the gamma signal +// in the range fAlphalow < |alpha| < alphasig +// + +Bool_t MHFindSignificance::SigmaVsAlpha(TH1 *fhist, Double_t alphamin, + Double_t alphamax, Int_t degree, Bool_t print) +{ + //*fLog << "MHFindSignificance::SigmaVsAlpha;" << endl; + + fHistOrig = fhist; + + fHist = (TH1*)fHistOrig->Clone("|alpha| plot"); + fHist->Sumw2(); + fHist->SetNameTitle("Alpha", "alpha plot"); + fHist->SetXTitle("|alpha| [\\circ]"); + fHist->SetYTitle("Counts"); + + + fAlphamin = alphamin; + fAlphamax = alphamax; + fAlphammm = (alphamin+alphamax)/2.0; + fDegree = degree; + + + //-------------------------------------------- + // fit a polynomial in the background region + + //*fLog << "MHFindSignificance::SigmaVsAlpha calling FitPolynomial()" + // << endl; + if ( !FitPolynomial() ) + { + *fLog << "MHFindSignificance::SigmaVsAlpha; FitPolynomial() failed" + << endl; + return kFALSE; + } + + + //-------------------------------------------- + // loop over different signal regions + + Int_t nsteps = 15; + + fSigVsAlpha = new TH1D("SigVsAlpha","Sigma vs Alpha", nsteps, 0.0, alphamin); + fSigVsAlpha->SetXTitle("upper edge of signal region in |alpha| [\\circ]"); + fSigVsAlpha->SetYTitle("Significance of gamma signal"); + + for (Int_t i=1; i<=nsteps; i++) + { + fAlphasig = fSigVsAlpha->GetBinCenter(i); + + if ( !DetExcess() ) + { + *fLog << "MHFindSignificance::SigmaVsAlpha; DetExcess() failed" << endl; + continue; + } + + + *fLog << "before SigmaLiMa : fNon, fNoff, fGamma = " << fNon << ", " + << fNoff << ", " << fGamma << endl; + + Double_t siglima = 0.0; + if ( !SigmaLiMa(fNon, fNoff, fGamma, &siglima) ) + { + *fLog << "MHFindSignificance::SigmaVsAlpha; SigmaLiMa() failed" << endl; + continue; + } + + fdNex = fNex / siglima; + fSigVsAlpha->SetBinContent(i, siglima); + + if (print) + PrintPoly(); + } + + //-------------------------------------------- + // plot significance versus alphasig + + TCanvas *ccc = new TCanvas("SigVsAlpha", "Sigma vs Alpha", 600, 600); + + gROOT->SetSelectedPad(NULL); + + ccc->cd(); + fSigVsAlpha->DrawCopy(); + + ccc->Modified(); + ccc->Update(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// FitPolynomial +// +// - create a clone 'fHist' of the |alpha| distribution 'fHistOrig' +// - fit a polynomial of degree 'fDegree' to the alpha distribution +// 'fHist' in the region alphamin < |alpha| < alphamax +// +// in pathological cases the histogram is rebinned before fitting +// (this is done only if fRebin is kTRUE) +// +// if the highest coefficient of the polynomial is compatible with zero +// the fit is repeated with a polynomial of lower degree +// (this is done only if fReduceDegree is kTRUE) +// +// +Bool_t MHFindSignificance::FitPolynomial() +{ + //-------------------------------------------------- + // check the histogram : + // - calculate initial values of the parameters + // - check for bins with zero entries + // - set minimum errors + // - save the original errors + // - set errors huge outside the fit range + // (in 'fcnpoly' points with huge errors will be ignored) + + + Double_t dummy = 1.e20; + + Double_t mean; + Double_t rms; + Double_t nclose; + Double_t nfar; + Double_t a2init = 0.0; + TArrayD saveError; + + Int_t nbins; + Int_t nrebin = 1; + + //---------------- start while loop for rebinning ----------------- + while(1) + { + + fNzero = 0; + fMbins = 0; + fMlow = 0; + fNbgtot = 0.0; + + fAlphami = 10000.0; + fAlphamm = 10000.0; + fAlphama = -10000.0; + + mean = 0.0; + rms = 0.0; + nclose = 0.0; + nfar = 0.0; + + nbins = fHist->GetNbinsX(); + saveError.Set(nbins); + + for (Int_t i=1; i<=nbins; i++) + { + saveError[i-1] = fHist->GetBinError(i); + + // bin should be completely contained in the fit range + // (fAlphamin, fAlphamax) + Double_t xlo = fHist->GetBinLowEdge(i); + Double_t xup = fHist->GetBinLowEdge(i+1); + + if ( xlo >= fAlphamin-fEps && xlo <= fAlphamax+fEps && + xup >= fAlphamin-fEps && xup <= fAlphamax+fEps ) + { + fMbins++; + + if ( xlo < fAlphami ) + fAlphami = xlo; + + if ( xup > fAlphama ) + fAlphama = xup; + + Double_t content = fHist->GetBinContent(i); + fNbgtot += content; + + mean += content; + rms += content*content; + + // count events in low-alpha and high-alpha region + if ( xlo >= fAlphammm-fEps && xup >= fAlphammm-fEps) + { + nfar += content; + if ( xlo < fAlphamm ) + fAlphamm = xlo; + if ( xup < fAlphamm ) + fAlphamm = xup; + } + else + { + nclose += content; + if ( xlo > fAlphamm ) + fAlphamm = xlo; + if ( xup > fAlphamm ) + fAlphamm = xup; + } + + // count bins with zero entry + if (content <= 0.0) + fNzero++; + + // set minimum error + if (content < 9.0) + { + fMlow += 1; + fHist->SetBinError(i, 3.0); + } + + //*fLog << "Take : i, content, error = " << i << ", " + // << fHist->GetBinContent(i) << ", " + // << fHist->GetBinError(i) << endl; + + continue; + } + // bin is not completely contained in the fit range : set error huge + + fHist->SetBinError(i, dummy); + + //*fLog << "Omit : i, content, error = " << i << ", " + // << fHist->GetBinContent(i) << ", " << fHist->GetBinError(i) + // << endl; + + } + + // mean of entries/bin in the fit range + if (fMbins > 0) + { + mean /= ((Double_t) fMbins); + rms /= ((Double_t) fMbins); + } + + rms = sqrt( rms - mean*mean ); + + // if there are no events in the background region + // there is no reason for rebinning + // and this is the condition for assuming a constant background (= 0) + if (mean <= 0.0) + break; + + Double_t helpmi = fAlphami*fAlphami*fAlphami; + Double_t helpmm = fAlphamm*fAlphamm*fAlphamm; + Double_t helpma = fAlphama*fAlphama*fAlphama; + Double_t help = (helpma-helpmm) * (fAlphamm-fAlphami) + - (helpmm-helpmi) * (fAlphama-fAlphamm); + if (help != 0.0) + a2init = ( (fAlphamm-fAlphami)*nfar - (fAlphama-fAlphamm)*nclose ) + * 1.5 * fHist->GetBinWidth(1) / help; + else + a2init = 0.0; + + + //-------------------------------------------- + // rebin the histogram + // - if a bin has no entries + // - or if there are too many bins with too few entries + // - or if the new bin width would exceed half the size of the + // signal region + + if ( !fRebin || + ( fNzero <= 0 && (Double_t)fMlow<0.05*(Double_t)fMbins ) || + (Double_t)(nrebin+1)/(Double_t)nrebin * fHist->GetBinWidth(1) + > fAlphasig/2.0 ) + { + //*fLog << "before break" << endl; + break; + } + + nrebin += 1; + //if (fHist) + delete fHist; + fHist = NULL; + + *fLog << "MHFindSignificance::FitPolynomial; rebin the |alpha| plot, grouping " + << nrebin << " bins together" << endl; + + // TH1::Rebin doesn't work properly + //fHist = fHistOrig->Rebin(nrebin, "Rebinned"); + // use private routine RebinHistogram() + fHist = new TH1F; + fHist->Sumw2(); + fHist->SetNameTitle("Rebinned", "Rebinned alpha plot"); + + + // do rebinning such that x0 remains a lower bin edge + Double_t x0 = 0.0; + if ( !RebinHistogram(x0, nrebin) ) + { + *fLog << "MHFindSignificance::FitPolynomial; RebinHistgram() failed" + << endl; + return kFALSE; + } + + fHist->SetXTitle("|alpha| [\\circ]"); + fHist->SetYTitle("Counts"); + + } + //---------------- end of while loop for rebinning ----------------- + + + // if there are still too many bins with too few entries don't fit + // and assume a constant background + + fConstantBackg = kFALSE; + if ( fNzero > 0 || (Double_t)fMlow>0.05*(Double_t)fMbins ) + { + *fLog << "MHFindSignificance::FitPolynomial; polynomial fit not possible, fNzero, fMlow, fMbins = " + << fNzero << ", " << fMlow << ", " << fMbins << endl; + *fLog << " assume a constant background" << endl; + + fConstantBackg = kTRUE; + fDegree = 0; + + TString funcname = "Poly"; + Double_t xmin = 0.0; + Double_t xmax = 90.0; + + TString formula = "[0]"; + + fPoly = new TF1(funcname, formula, xmin, xmax); + TList *funclist = fHist->GetListOfFunctions(); + funclist->Add(fPoly); + + //-------------------- + Int_t nparfree = 1; + fChisq = 0.0; + fNdf = fMbins - nparfree; + fProb = 0.0; + fIstat = 0; + + fValues.Set(1); + fErrors.Set(1); + + Double_t val, err; + val = mean; + err = sqrt( mean / (Double_t)fMbins ); + + fPoly->SetParameter(0, val); + fPoly->SetParError (0, err); + + fValues[0] = val; + fErrors[0] = err; + + fEma[0][0] = err*err; + fCorr[0][0] = 1.0; + //-------------------- + + //-------------------------------------------------- + // reset the errors of the points in the histogram + for (Int_t i=1; i<=nbins; i++) + { + fHist->SetBinError(i, saveError[i-1]); + } + + + return kTRUE; + } + + + //=========== start loop for reducing the degree ================== + // of the polynomial + while (1) + { + + //-------------------------------------------------- + // prepare fit of a polynomial : (a0 + a1*x + a2*x**2 + a3*x**3 + ...) + + TString funcname = "Poly"; + Double_t xmin = 0.0; + Double_t xmax = 90.0; + + TString formula = "[0]"; + TString bra1 = "+["; + TString bra2 = "]"; + TString xpower = "*x"; + TString newpower = "*x"; + for (Int_t i=1; i<=fDegree; i++) + { + formula += bra1; + formula += i; + formula += bra2; + formula += xpower; + + xpower += newpower; + } + + //*fLog << "FitPolynomial : formula = " << formula << endl; + + fPoly = new TF1(funcname, formula, xmin, xmax); + TList *funclist = fHist->GetListOfFunctions(); + funclist->Add(fPoly); + + //------------------------ + // attention : the dimensions must agree with those in CallMinuit() + + char parname[80][100]; + Double_t vinit[80]; + Double_t step[80]; + Double_t limlo[80]; + Double_t limup[80]; + Int_t fix[80]; + + UInt_t npar = fDegree+1; + + for (UInt_t j=0; jGetEntries(); + + + // use the subsequernt loop if you want to apply the + // constraint : uneven derivatives (at alpha=0) = zero + for (UInt_t j=1; jmnstat(fmin, fedm, errdef, npari, nparx, fIstat); + + //*fLog << "MHFindSignificance::FitPolynomial; fmin, fedm, errdef, npari, nparx, fIstat = " + // << fmin << ", " << fedm << ", " << errdef << ", " << npari + // << ", " << nparx << ", " << fIstat << endl; + + + //------------------- + // store the results + + Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0; + fChisq = fmin; + fNdf = fMbins - nparfree; + fProb = TMath::Prob(fChisq, fNdf); + + + // get fitted parameter values and errors + fValues.Set(fDegree+1); + fErrors.Set(fDegree+1); + + for (Int_t j=0; j<=fDegree; j++) + { + Double_t val, err; + if (gMinuit) + gMinuit->GetParameter(j, val, err); + + fPoly->SetParameter(j, val); + fPoly->SetParError(j, err); + + fValues[j] = val; + fErrors[j] = err; + } + + //-------------------------------------------------- + // if the highest coefficient (j0) of the polynomial + // is consistent with zero reduce the degree of the polynomial + + Int_t j0 = 0; + for (Int_t j=fDegree; j>1; j--) + { + // ignore fixed parameters + if (fErrors[j] == 0.0) + continue; + + // this is the highest coefficient + j0 = j; + break; + } + + if ( !fReduceDegree || j0 == 0 || fabs(fValues[j0]) > fErrors[j0] ) + break; + + // reduce the degree of the polynomial + *fLog << "MHFindSignificance::FitPolynomial; reduce the degree of the polynomial from " + << fDegree << " to " << (j0-2) << endl; + fDegree = j0 - 2; + + funclist->Remove(fPoly); + //if (fPoly) + delete fPoly; + fPoly = NULL; + + // delete the Minuit object in order to have independent starting + // conditions for the next minimization + //if (gMinuit) + delete gMinuit; + gMinuit = NULL; + } + //=========== end of loop for reducing the degree ================== + // of the polynomial + + + //-------------------------------------------------- + // get the error matrix of the fitted parameters + + + if (fIstat >= 1) + { + // error matrix was calculated + if (gMinuit) + gMinuit->mnemat(&fEmat[0][0], fNdim); + + // copy covariance matrix into a matrix which includes also the fixed + // parameters + TString name; + Double_t bnd1, bnd2, val, err; + Int_t jvarbl; + Int_t kvarbl; + for (Int_t j=0; j<=fDegree; j++) + { + if (gMinuit) + gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl); + for (Int_t k=0; k<=fDegree; k++) + { + if (gMinuit) + gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl); + if (jvarbl == 0 || kvarbl == 0) + fEma[j][k] = 0.0; + else + fEma[j][k] = fEmat[jvarbl-1][kvarbl-1]; + } + } + } + else + { + // error matrix was not calculated, construct it + *fLog << "MHFindSignificance::FitPolynomial; error matrix not defined" + << endl; + for (Int_t j=0; j<=fDegree; j++) + { + for (Int_t k=0; k<=fDegree; k++) + { + fEma[j][k] = 0.0; + } + fEma[j][j] = fErrors[j]*fErrors[j]; + } + } + + + //-------------------------------------------------- + // calculate correlation matrix + for (Int_t j=0; j<=fDegree; j++) + { + for (Int_t k=0; k<=fDegree; k++) + { + if (fEma[j][j]*fEma[k][k] != 0.0) + fCorr[j][k] = fEma[j][k] / sqrt( fEma[j][j]*fEma[k][k] ); + else + fCorr[j][k] = 0.0; + } + } + + + //-------------------------------------------------- + // reset the errors of the points in the histogram + for (Int_t i=1; i<=nbins; i++) + { + fHist->SetBinError(i, saveError[i-1]); + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// ReBinHistogram +// +// rebin the histogram 'fHistOrig' by grouping 'nrebin' bins together +// put the result into the histogram 'fHist' +// the rebinning is made such that 'x0' remains a lower bound of a bin +// + +Bool_t MHFindSignificance::RebinHistogram(Double_t x0, Int_t nrebin) +{ + //----------------------------------------- + // search bin i0 which has x0 as lower edge + + Int_t i0 = -1; + Int_t nbold = fHistOrig->GetNbinsX(); + for (Int_t i=1; i<=nbold; i++) + { + if ( fabs(fHistOrig->GetBinLowEdge(i) - x0) < 1.e-4 ) + { + i0 = i; + break; + } + } + + if (i0 == -1) + { + i0 = 1; + *fLog << "MHFindsignificance::Rebin; no bin found with " << x0 + << " as lower edge, start rebinning with bin 1" << endl; + } + + Int_t istart = i0 - nrebin * ( (i0-1)/nrebin ); + + //----------------------------------------- + // get new bin edges + + Int_t nbnew = (nbold-istart+1) / nrebin; + Double_t xmin = fHistOrig->GetBinLowEdge(istart); + Double_t xmax = xmin + + ((Double_t)nbnew) * ((Double_t)nrebin) * fHistOrig->GetBinWidth(1); + fHist->SetBins(nbnew, xmin, xmax); + + *fLog << "MHFindSignificance::ReBin; x0, i0, nbold, nbnew, xmin, xmax = " + << x0 << ", " << i0 << ", " << nbold << ", " << nbnew << ", " + << xmin << ", " << xmax << endl; + + //----------------------------------------- + // get new bin entries + + for (Int_t i=1; i<=nbnew; i++) + { + Int_t j = nrebin*(i-1) + istart; + + Double_t content = 0.0; + Double_t error2 = 0.0; + for (Int_t k=0; kGetBinContent(j+k); + error2 += fHistOrig->GetBinError(j+k) * fHistOrig->GetBinError(j+k); + } + fHist->SetBinContent(i, content); + fHist->SetBinError (i, sqrt(error2)); + } + fHist->SetEntries( fHistOrig->GetEntries() ); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// FitGaussPoly +// +// fits a (Gauss + polynomial function) to the alpha distribution 'fhist' +// +// +Bool_t MHFindSignificance::FitGaussPoly() +{ + *fLog << "Entry FitGaussPoly" << endl; + + //-------------------------------------------------- + // check the histogram : + // - calculate initial values of the parameters + // - check for bins with zero entries + // - set minimum errors + // - save the original errors + // - set errors huge outside the fit range + // (in 'fcnpoly' points with huge errors will be ignored) + + + Double_t dummy = 1.e20; + + fGNzero = 0; + fGMbins = 0; + + //------------------------------------------ + // if a constant background has been assumed (due to low statistics) + // fit only in the signal region + if ( !fConstantBackg ) + { + fAlphalow = 0.0; + fAlphahig = fAlphamax; + } + else + { + fAlphalow = 0.0; + fAlphahig = 2.0*fAlphasig>25.0 ? 25.0 : 2.0*fAlphasig; + } + //------------------------------------------ + + + fAlphalo = 10000.0; + fAlphahi = -10000.0; + + + Int_t nbins = fHist->GetNbinsX(); + TArrayD saveError(nbins); + + for (Int_t i=1; i<=nbins; i++) + { + saveError[i-1] = fHist->GetBinError(i); + + // bin should be completely contained in the fit range + // (fAlphalow, fAlphahig) + Double_t xlo = fHist->GetBinLowEdge(i); + Double_t xup = fHist->GetBinLowEdge(i+1); + + if ( xlo >= fAlphalow-fEps && xlo <= fAlphahig+fEps && + xup >= fAlphalow-fEps && xup <= fAlphahig+fEps ) + { + fGMbins++; + + if ( xlo < fAlphalo ) + fAlphalo = xlo; + + if ( xup > fAlphahi ) + fAlphahi = xup; + + Double_t content = fHist->GetBinContent(i); + + + // count bins with zero entry + if (content <= 0.0) + fGNzero++; + + // set minimum error + if (content < 9.0) + fHist->SetBinError(i, 3.0); + + //*fLog << "Take : i, content, error = " << i << ", " + // << fHist->GetBinContent(i) << ", " + // << fHist->GetBinError(i) << endl; + + continue; + } + // bin is not completely contained in the fit range : set error huge + + fHist->SetBinError(i, dummy); + + //*fLog << "Omit : i, content, error = " << i << ", " + // << fHist->GetBinContent(i) << ", " << fHist->GetBinError(i) + // << endl; + + } + + + // if a bin has no entries don't fit + if (fGNzero > 0) + { + *fLog << "MHFindSignificance::FitGaussPoly; out of " << fGMbins + << " bins there are " << fGNzero + << " bins with zero entry" << endl; + + fGPoly = NULL; + return kFALSE; + } + + + //-------------------------------------------------- + // prepare fit of a (polynomial+Gauss) : + // (a0 + a1*x + a2*x**2 + a3*x**3 + ...) + A*exp( -0.5*((x-x0)/sigma)**2 ) + + TString funcname = "PolyGauss"; + Double_t xmin = 0.0; + Double_t xmax = 90.0; + + TString bra1 = "["; + TString bra2 = "]"; + TString xpower = "*x"; + TString newpower = "*x"; + + TString formulaBackg = "[0]"; + for (Int_t i=1; i<=fDegree; i++) + { + formulaBackg += "+"; + formulaBackg += bra1; + formulaBackg += i; + formulaBackg += bra2; + formulaBackg += xpower; + + xpower += newpower; + } + + TString formulaGauss = bra1; + formulaGauss += fDegree+1; + formulaGauss += bra2; + formulaGauss += "/"; + formulaGauss += bra1; + formulaGauss += fDegree+3; + formulaGauss += bra2; + formulaGauss += "*exp(-0.5*((x-"; + formulaGauss += bra1; + formulaGauss += fDegree+2; + formulaGauss += bra2; + formulaGauss += ")/"; + formulaGauss += bra1; + formulaGauss += fDegree+3; + formulaGauss += bra2; + formulaGauss += ")^2)"; + + TString formula = formulaBackg; + formula += "+"; + formula += formulaGauss; + + *fLog << "FitGaussPoly : formulaBackg = " << formulaBackg << endl; + *fLog << "FitGaussPoly : formulaGauss = " << formulaGauss << endl; + *fLog << "FitGaussPoly : formula = " << formula << endl; + + fGPoly = new TF1(funcname, formula, xmin, xmax); + TList *funclist = fHist->GetListOfFunctions(); + funclist->Add(fGPoly); + + fGBackg = new TF1("Backg", formulaBackg, xmin, xmax); + funclist->Add(fGBackg); + + //------------------------ + // attention : the dimensions must agree with those in CallMinuit() + + char parname[80][100]; + Double_t vinit[80]; + Double_t step[80]; + Double_t limlo[80]; + Double_t limup[80]; + Int_t fix[80]; + + Int_t npar = fDegree+1 + 3; + + // take as initial values for the polynomial + // the result from the polynomial fit + for (Int_t j=0; j<=fDegree; j++) + { + vinit[j] = fPoly->GetParameter(j); + } + + + + Double_t sigma = 8.0; + vinit[fDegree+1] = 2.0 * fNex * fHist->GetBinWidth(1) / + sqrt( 2.0*TMath::Pi() ); + vinit[fDegree+2] = 0.0; + vinit[fDegree+3] = sigma; + + *fLog << "FitGaussPoly : starting value for Gauss-amplitude = " + << vinit[fDegree+1] << endl; + + for (Int_t j=0; jGetEntries(); + + // limit the sigma of the Gauss function + limlo[fDegree+3] = 0.0; + limup[fDegree+3] = 20.0; + + + // use the subsequernt loop if you want to apply the + // constraint : uneven derivatives (at alpha=0) = zero + for (Int_t j=1; j<=fDegree; j += 2) + { + vinit[j] = 0.0; + step[j] = 0.0; + fix[j] = 1; + } + + // fix position of Gauss function + vinit[fDegree+2] = 0.0; + step[fDegree+2] = 0.0; + fix[fDegree+2] = 1; + + // if a constant background has been assumed (due to low statistics) + // fix the background + if (fConstantBackg) + { + step[0] = 0.0; + fix[0] = 1; + } + + TString method = "Migrad"; + TObject *objectfit = fHist; + Bool_t nulloutput = kFALSE; + + MMinuitInterface inter; + Bool_t rc = inter.CallMinuit( fcnpolygauss, npar, parname, vinit, step, + limlo, limup, fix, objectfit, method, + nulloutput); + + if (rc != 0) + { + // *fLog << "MHFindSignificance::FitGaussPoly; (polynomial+Gauss) fit failed" + // << endl; + // return kFALSE; + } + + + //------------------- + // get status of the minimization + Double_t fmin; + Double_t fedm; + Double_t errdef; + Int_t npari; + Int_t nparx; + + if (gMinuit) + gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fGIstat); + + *fLog << "MHFindSignificance::FitGaussPoly; fmin, fedm, errdef, npari, nparx, fGIstat = " + << fmin << ", " << fedm << ", " << errdef << ", " << npari + << ", " << nparx << ", " << fGIstat << endl; + + + //------------------- + // store the results + + Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0; + fGChisq = fmin; + fGNdf = fGMbins - nparfree; + fGProb = TMath::Prob(fGChisq, fGNdf); + + + // get fitted parameter values and errors + fGValues.Set(npar); + fGErrors.Set(npar); + + for (Int_t j=0; jGetParameter(j, val, err); + + fGPoly->SetParameter(j, val); + fGPoly->SetParError(j, err); + + fGValues[j] = val; + fGErrors[j] = err; + + if (j <=fDegree) + { + fGBackg->SetParameter(j, val); + fGBackg->SetParError(j, err); + } + } + + fSigmaGauss = fGValues[fDegree+3]; + fdSigmaGauss = fGErrors[fDegree+3]; + // fitted total number of excess events + fNexGauss = fGValues[fDegree+1] * sqrt(2.0*TMath::Pi()) / + ( 2.0 * fHist->GetBinWidth(1) ); + fdNexGauss = fNexGauss * fGErrors[fDegree+1]/fGValues[fDegree+1]; + + //-------------------------------------------------- + // get the error matrix of the fitted parameters + + + if (fGIstat >= 1) + { + // error matrix was calculated + if (gMinuit) + gMinuit->mnemat(&fGEmat[0][0], fGNdim); + + // copy covariance matrix into a matrix which includes also the fixed + // parameters + TString name; + Double_t bnd1, bnd2, val, err; + Int_t jvarbl; + Int_t kvarbl; + for (Int_t j=0; jmnpout(j, name, val, err, bnd1, bnd2, jvarbl); + for (Int_t k=0; kmnpout(k, name, val, err, bnd1, bnd2, kvarbl); + if (jvarbl == 0 || kvarbl == 0) + fGEma[j][k] = 0.0; + else + fGEma[j][k] = fGEmat[jvarbl-1][kvarbl-1]; + } + } + } + else + { + // error matrix was not calculated, construct it + *fLog << "MHFindSignificance::FitPolynomial; error matrix not defined" + << endl; + for (Int_t j=0; jSetBinError(i, saveError[i-1]); + } + + return kTRUE; + +} + +// -------------------------------------------------------------------------- +// +// DetExcess +// +// using the result of the polynomial fit (fValues), DetExcess determines +// +// - the total number of events in the signal region (fNon) +// - the number of backgound events in the signal region (fNbg) +// - the number of excess events (fNex) +// - the effective number of background events (fNoff), and fGamma : +// fNbg = fGamma * fNoff; fdNbg = fGamma * sqrt(fNoff); +// +// It assumed that the polynomial is defined as +// a0 + a1*x + a2*x**2 + a3*x**3 + .. +// +// and that the alpha distribution has the range 0 < alpha < 90 degrees +// + +Bool_t MHFindSignificance::DetExcess() +{ + //*fLog << "MHFindSignificance::DetExcess;" << endl; + + //-------------------------------------------- + // calculate the total number of events (fNon) in the signal region + + fNon = 0.0; + fdNon = 0.0; + + Double_t alphaup = -1000.0; + Double_t binwidth = fHist->GetBinWidth(1); + + Int_t nbins = fHist->GetNbinsX(); + for (Int_t i=1; i<=nbins; i++) + { + Double_t xlo = fHist->GetBinLowEdge(i); + Double_t xup = fHist->GetBinLowEdge(i+1); + + // bin must be completely contained in the signal region + if ( xlo <= (fAlphasig+fEps) && xup <= (fAlphasig+fEps) ) + { + Double_t width = fabs(xup-xlo); + if (fabs(width-binwidth) > fEps) + { + *fLog << "MHFindSignificance::DetExcess; alpha plot has variable binning, which is not allowed" + << endl; + return kFALSE; + } + + if (xup > alphaup) + alphaup = xup; + + fNon += fHist->GetBinContent(i); + fdNon += fHist->GetBinError(i) * fHist->GetBinError(i); + } + } + fdNon = sqrt(fdNon); + + // the actual signal range is : + if (alphaup == -1000.0) + return kFALSE; + + fAlphasi = alphaup; + + //*fLog << "fAlphasi, fNon, fdNon, binwidth, fDegree = " << fAlphasi << ", " + // << fNon << ", " << fdNon << ", " << binwidth << ", " + // << fDegree << endl; + + //-------------------------------------------- + // calculate the number of background events (fNbg) in the signal region + // and its error (fdNbg) + + Double_t fac = 1.0/binwidth; + + fNbg = 0.0; + Double_t altothejplus1 = fAlphasi; + for (Int_t j=0; j<=fDegree; j++) + { + fNbg += fValues[j] * altothejplus1 / ((Double_t)(j+1)); + altothejplus1 *= fAlphasi; + } + fNbg *= fac; + + // derivative of Nbg + Double_t facj; + Double_t fack; + + Double_t sum = 0.0; + altothejplus1 = fAlphasi; + for (Int_t j=0; j<=fDegree; j++) + { + facj = altothejplus1 / ((Double_t)(j+1)); + + Double_t altothekplus1 = fAlphasi; + for (Int_t k=0; k<=fDegree; k++) + { + fack = altothekplus1 / ((Double_t)(k+1)); + + sum += facj * fack * fEma[j][k]; + altothekplus1 *= fAlphasi; + } + altothejplus1 *= fAlphasi; + } + sum *= fac*fac; + + if (sum < 0.0) + { + *fLog << "MHFindsignificance::DetExcess; error squared is negative" + << endl; + return kFALSE; + } + + fdNbg = sqrt(sum); + + + //-------------------------------------------- + // AS A CHECK : + // calculate the number of background events (fNbgtotFitted) in the + // background region, and its error (fdNbgtotFitted) + // expect fdnbg to be approximately equal to sqrt(fNbgtotFitted) + + Double_t fNmi = 0.0; + altothejplus1 = fAlphami; + for (Int_t j=0; j<=fDegree; j++) + { + fNmi += fValues[j] * altothejplus1 / ((Double_t)(j+1)); + altothejplus1 *= fAlphami; + } + fNmi *= fac; + + Double_t fNma = 0.0; + altothejplus1 = fAlphama; + for (Int_t j=0; j<=fDegree; j++) + { + fNma += fValues[j] * altothejplus1 / ((Double_t)(j+1)); + altothejplus1 *= fAlphama; + } + fNma *= fac; + + fNbgtotFitted = fNma - fNmi; + + //---------------------- + + sum = 0.0; + Double_t altothejma = fAlphama; + Double_t altothejmi = fAlphami; + for (Int_t j=0; j<=fDegree; j++) + { + facj = (altothejma-altothejmi) / ((Double_t)(j+1)); + + Double_t altothekma = fAlphama; + Double_t altothekmi = fAlphami; + for (Int_t k=0; k<=fDegree; k++) + { + fack = (altothekma-altothekmi) / ((Double_t)(k+1)); + + sum += facj * fack * fEma[j][k]; + altothekma *= fAlphama; + altothekmi *= fAlphami; + } + altothejma *= fAlphama; + altothejmi *= fAlphami; + } + sum *= fac*fac; + + fdNbgtotFitted = sqrt(sum); + if ( fabs(fdNbgtotFitted - sqrt(fNbgtotFitted)) > 0.2 * sqrt(fNbgtotFitted) ) + { + *fLog << "MHFindSignificance::DetExcess; error of calculated number of background events (in the background region) does not agree with the expectation :" + << endl; + *fLog << " fNbgtotFitted, fdNbgtotFitted = " + << fNbgtotFitted << ", " << fdNbgtotFitted + << ", expected : " << sqrt(fNbgtotFitted) << endl; + } + + + //-------------------------------------------- + // calculate the number of excess events in the signal region + + fNex = fNon - fNbg; + + //-------------------------------------------- + // calculate the effective number of background events (fNoff) , and fGamma : + // fNbg = fGamma * fNoff; dfNbg = fGamma * sqrt(fNoff); + + if (fNbg < 0.0) + { + *fLog << "MHFindSignificamce::DetExcess; number of background events is negative, fNbg, fdNbg = " + << fNbg << ", " << fdNbg << endl; + + fGamma = 1.0; + fNoff = 0.0; + return kFALSE; + } + + if (fNbg > 0.0) + { + fGamma = fdNbg*fdNbg / fNbg; + fNoff = fNbg*fNbg / (fdNbg*fdNbg); + } + else + { + fGamma = 1.0; + fNoff = 0.0; + } + + //*fLog << "Exit DetExcess()" << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// SigmaLiMa +// +// calculates the significance according to Li & Ma +// ApJ 272 (1983) 317 +// +Bool_t MHFindSignificance::SigmaLiMa(Double_t non, Double_t noff, + Double_t gamma, Double_t *siglima) +{ + if (gamma <= 0.0 || non <= 0.0 || noff <= 0.0) + { + *siglima = 0.0; + return kFALSE; + } + + Double_t help1 = non * log( (1.0+gamma)*non / (gamma*(non+noff)) ); + Double_t help2 = noff * log( (1.0+gamma)*noff / ( non+noff ) ); + *siglima = sqrt( 2.0 * (help1+help2) ); + + Double_t nex = non - gamma*noff; + if (nex < 0.0) + *siglima = - *siglima; + + //*fLog << "MHFindSignificance::SigmaLiMa; non, noff, gamma, *siglima = " + // << non << ", " << noff << ", " << gamma << ", " << *siglima << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +Bool_t MHFindSignificance::DrawFit(const Option_t *opt) +{ + *fLog << "entry DrawFit" << endl; + + + //TCanvas *fCanvas = new TCanvas("Alpha", "Alpha plot", 600, 600); + fCanvas = new TCanvas("Alpha", "Alpha plot", 600, 600); + + //gStyle->SetOptFit(1011); + gROOT->SetSelectedPad(NULL); + + fCanvas->cd(); + if (fHist == NULL) + *fLog << "MHFindSignificance::DrawFit; fHist = NULL" << endl; + + if (fHist) + { + fHist->DrawCopy(); + } + + TF1 *fpoly = fHist->GetFunction("Poly"); + if (fpoly == NULL) + *fLog << "MHFindSignificance::DrawFit; fpoly = NULL" << endl; + + if (fpoly) + { + // 2, 1 is red and solid + fpoly->SetLineColor(2); + fpoly->SetLineStyle(1); + fpoly->SetLineWidth(2); + fpoly->DrawCopy("same"); + } + + if (fFitGauss) + { + TF1 *fpolygauss = fHist->GetFunction("PolyGauss"); + if (fpolygauss == NULL) + *fLog << "MHFindSignificance::DrawFit; fpolygauss = NULL" << endl; + + if (fpolygauss) + { + // 4, 1 is blue and solid + fpolygauss->SetLineColor(4); + fpolygauss->SetLineStyle(1); + fpolygauss->SetLineWidth(4); + fpolygauss->DrawCopy("same"); + } + + TF1 *fbackg = fHist->GetFunction("Backg"); + if (fbackg == NULL) + *fLog << "MHFindSignificance::DrawFit; fbackg = NULL" << endl; + + if (fbackg) + { + // 6, 4 is pink and dotted + fbackg->SetLineColor(4); + fbackg->SetLineStyle(4); + fbackg->SetLineWidth(4); + fbackg->DrawCopy("same"); + } + } + + *fLog << "DrawFit : before text" << endl; + + //------------------------------- + // print results onto the figure + TPaveText *pt = new TPaveText(0.30, 0.35, 0.70, 0.90, "NDC"); + char tx[100]; + + sprintf(tx, "Results of polynomial fit (order %2d) :", fDegree); + TText *t1 = pt->AddText(tx); + t1->SetTextSize(0.03); + t1->SetTextColor(2); + + sprintf(tx, " (%6.2f< |alpha| <%6.2f [\\circ])", fAlphami, fAlphama); + pt->AddText(tx); + + sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f", + fChisq, fNdf, fProb); + pt->AddText(tx); + + sprintf(tx, " Nbgtot(fit) = %8.1f #pm %8.1f", + fNbgtotFitted, fdNbgtotFitted); + pt->AddText(tx); + + sprintf(tx, " Nbgtot(meas) = %8.1f", fNbgtot); + pt->AddText(tx); + + + //sprintf(tx, " "); + //pt->AddText(tx); + + //-------------- + sprintf(tx, "Results for |alpha|< %6.2f [\\circ] :", fAlphasi); + TText *t6 = pt->AddText(tx); + t6->SetTextSize(0.03); + t6->SetTextColor(8); + + sprintf(tx, " Non = %8.1f #pm %8.1f", fNon, fdNon); + pt->AddText(tx); + + sprintf(tx, " Nex = %8.1f #pm %8.1f", fNex, fdNex); + pt->AddText(tx); + + sprintf(tx, " Nbg = %8.1f #pm %8.1f, gamma = %6.1f", + fNbg, fdNbg, fGamma); + pt->AddText(tx); + + Double_t ratio = fNbg>0.0 ? fNex/fNbg : 0.0; + sprintf(tx, " Significance = %6.2f, Nex/Nbg = %6.2f", + fSigLiMa, ratio); + pt->AddText(tx); + + //sprintf(tx, " "); + //pt->AddText(tx); + + //-------------- + if (fFitGauss) + { + sprintf(tx, "Results of (polynomial+Gauss) fit :"); + TText *t7 = pt->AddText(tx); + t7->SetTextSize(0.03); + t7->SetTextColor(4); + + sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f", + fGChisq, fGNdf, fGProb); + pt->AddText(tx); + + sprintf(tx, " Sigma of Gauss = %8.1f #pm %8.1f [\\circ]", + fSigmaGauss, fdSigmaGauss); + pt->AddText(tx); + + sprintf(tx, " total no.of excess events = %8.1f #pm %8.1f", + fNexGauss, fdNexGauss); + pt->AddText(tx); + } + //-------------- + + pt->SetFillStyle(0); + pt->SetBorderSize(0); + pt->SetTextAlign(12); + + *fLog << "DrawFit : before pt->Draw()" << endl; + + pt->Draw(); + + fCanvas->Modified(); + fCanvas->Update(); + + *fLog << "exit DrawFit" << endl; + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// Print the results of the polynomial fit to the alpha distribution +// +// +void MHFindSignificance::PrintPoly(Option_t *o) +{ + *fLog << "---------------------------" << endl; + *fLog << "MHFindSignificance::PrintPoly :" << endl; + + *fLog << "fAlphami, fAlphama, fDegree, fAlphasi = " + << fAlphami << ", " << fAlphama << ", " << fDegree << ", " + << fAlphasi << endl; + + *fLog << "fMbins, fNzero, fIstat = " << fMbins << ", " + << fNzero << ", " << fIstat << endl; + + *fLog << "fChisq, fNdf, fProb = " << fChisq << ", " + << fNdf << ", " << fProb << endl; + + *fLog << "fNon, fNbg, fdNbg, fNbgtot, fNbgtotFitted, fdNbgtotFitted = " + << fNon << ", " << fNbg << ", " << fdNbg << ", " << fNbgtot + << ", " << fNbgtotFitted << ", " << fdNbgtotFitted << endl; + + Double_t sigtoback = fNbg>0.0 ? fNex/fNbg : 0.0; + *fLog << "fNex, fdNex, fGamma, fNoff, fSigLiMa, sigtoback = " + << fNex << ", " << fdNex << ", " << fGamma << ", " << fNoff + << ", " << fSigLiMa << ", " << sigtoback << endl; + + //------------------------------------ + // get errors + + Double_t eplus; + Double_t eminus; + Double_t eparab; + Double_t gcc; + Double_t errdiag; + + + if ( !fConstantBackg ) + { + *fLog << "parameter value error eplus eminus eparab errdiag gcc" + << endl; + + for (Int_t j=0; j<=fDegree; j++) + { + if (gMinuit) + gMinuit->mnerrs(j, eplus, eminus, eparab, gcc); + errdiag = sqrt(fEma[j][j]); + *fLog << j << " " << fValues[j] << " " << fErrors[j] << " " + << eplus << " " << eminus << " " << eparab << " " + << errdiag << " " << gcc << endl; + } + } + else + { + *fLog << "parameter value error errdiag " + << endl; + + for (Int_t j=0; j<=fDegree; j++) + { + errdiag = sqrt(fEma[j][j]); + *fLog << j << " " << fValues[j] << " " << fErrors[j] << " " + << errdiag << endl; + } + } + + //---------------------------------------- + *fLog << "Covariance matrix :" << endl; + for (Int_t j=0; j<=fDegree; j++) + { + *fLog << "j = " << j << " : "; + for (Int_t k=0; k<=fDegree; k++) + { + *fLog << fEma[j][k] << " "; + } + *fLog << endl; + } + + *fLog << "Correlation matrix :" << endl; + for (Int_t j=0; j<=fDegree; j++) + { + *fLog << "j = " << j << " : "; + for (Int_t k=0; k<=fDegree; k++) + { + *fLog << fCorr[j][k] << " "; + } + *fLog << endl; + } + + *fLog << "---------------------------" << endl; +} + +// -------------------------------------------------------------------------- +// +// Print the results of the (polynomial+Gauss) fit to the alpha distribution +// +// +void MHFindSignificance::PrintPolyGauss(Option_t *o) +{ + *fLog << "---------------------------" << endl; + *fLog << "MHFindSignificance::PrintPolyGauss :" << endl; + + *fLog << "fAlphalo, fAlphahi = " + << fAlphalo << ", " << fAlphahi << endl; + + *fLog << "fGMbins, fGNzero, fGIstat = " << fGMbins << ", " + << fGNzero << ", " << fGIstat << endl; + + *fLog << "fGChisq, fGNdf, fGProb = " << fGChisq << ", " + << fGNdf << ", " << fGProb << endl; + + + //------------------------------------ + // get errors + + Double_t eplus; + Double_t eminus; + Double_t eparab; + Double_t gcc; + Double_t errdiag; + + *fLog << "parameter value error eplus eminus eparab errdiag gcc" + << endl; + for (Int_t j=0; j<=(fDegree+3); j++) + { + if (gMinuit) + gMinuit->mnerrs(j, eplus, eminus, eparab, gcc); + errdiag = sqrt(fGEma[j][j]); + *fLog << j << " " << fGValues[j] << " " << fGErrors[j] << " " + << eplus << " " << eminus << " " << eparab << " " + << errdiag << " " << gcc << endl; + } + + + *fLog << "Covariance matrix :" << endl; + for (Int_t j=0; j<=(fDegree+3); j++) + { + *fLog << "j = " << j << " : "; + for (Int_t k=0; k<=(fDegree+3); k++) + { + *fLog << fGEma[j][k] << " "; + } + *fLog << endl; + } + + *fLog << "Correlation matrix :" << endl; + for (Int_t j=0; j<=(fDegree+3); j++) + { + *fLog << "j = " << j << " : "; + for (Int_t k=0; k<=(fDegree+3); k++) + { + *fLog << fGCorr[j][k] << " "; + } + *fLog << endl; + } + + *fLog << "---------------------------" << endl; +} + +//============================================================================ + Index: /trunk/MagicSoft/Mars/mhist/MHFindSignificance.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHFindSignificance.h (revision 2300) +++ /trunk/MagicSoft/Mars/mhist/MHFindSignificance.h (revision 2300) @@ -0,0 +1,198 @@ +#ifndef MARS_MHFindSignificance +#define MARS_MHFindSignificance + +#ifdef MARS_MLogManip +#error Please make ensure that MLogManip.h are included _after_ MHFindSignificance.h +#endif + + +#ifndef MARS_MH +#include "MH.h" +#endif + +#include +#include +#include + +class TArrayD; +class TF1; + + +class MHFindSignificance : public MH +{ +private: + + TH1 *fHistOrig; // original plot of |alpha| (0.0 to 90.0 degrees) + TH1 *fHist; // copy of fHistOrig or rebinned histogram + + TH1D *fSigVsAlpha; + + Double_t fAlphamin; // requested lower limit of fit range + Double_t fAlphammm; // center of fit range + Double_t fAlphamax; // requested lower limit of fit range + + Double_t fAlphami; // actual lower limit of fit range + Double_t fAlphamm; // actual center of fit range + Double_t fAlphama; // actual upper limit of fit range + + Double_t fAlphasig; // requested signal range + Double_t fAlphasi; // actual signal range + + Double_t fAlphalow; // requested lower edge of signal range + Double_t fAlphalo; // actual lower edge of signal range + + Double_t fAlphahig; // requested upper edge of background range + Double_t fAlphahi; // actual upper edge of background range + + // number of events in signal region + Double_t fNon; // total number of events in signal region + Double_t fNbg; // number of background events in signal region + Double_t fNex; // number of excess events in signal region + + Double_t fdNon; + Double_t fdNbg; + Double_t fdNex; + + // number of events in background region + Double_t fNbgtot; // total number of events in background region + Double_t fNbgtotFitted; // fitted total no. of events in background region + Double_t fdNbgtotFitted; // fitted error of this number + + // effective number of background events + Double_t fNoff; + Double_t fGamma; // Nbg = Non - gamma * Noff + + Double_t fSigLiMa; // significance of gamma signal according to Li & Ma + + const static Double_t fEps = 1.e-4; // tolerance for floating point comparisons + + Bool_t fDraw; // if true : draw plots + Bool_t fFitGauss; // if true : do the (polynomial+Gauss fit) + Bool_t fRebin; // if true : allow rebinning of the alpha plot + Bool_t fReduceDegree; // if true : allow reducing of the order of the polynomial + + Bool_t fConstantBackg; // if set true if background fit is not possible + // due to low statistics + + TCanvas *fCanvas; + + Double_t fNexGauss; // total number of excess events + // (from fitted Gauss function) + Double_t fdNexGauss; // error of the total number of excess events + + Double_t fSigmaGauss; // sigma of fitted Gauss function + Double_t fdSigmaGauss; // error of this sigma + + //-------------------- + TF1 *fPoly; // polynomial function + Int_t fFitBad; // if != 0 fit failed + Int_t fDegree; // degree of polynomial to be fitted to the background + Int_t fNdf; // number of degrees of freedom of polynomial fit + Double_t fChisq; // chi squared of polynomial fit + Double_t fProb; // chi squared probability ofg polynomial fit + + TArrayD fValues; + TArrayD fErrors; + + const static Int_t fNdim = 6; + Double_t fEmat[fNdim][fNdim]; + Double_t fEma [fNdim][fNdim]; + Double_t fCorr[fNdim][fNdim]; + + Int_t fMbins; // number of bins in the fit range + Int_t fMlow; // number of bins in the fit range with too few entries + Int_t fNzero; // number of bins in the fit range with zero entry + Int_t fIstat; + + //-------------------- + + //-------------------- + TF1 *fGPoly; // (Gauss+polynomial) function + TF1 *fGBackg; // polynomial part of (Gauss+polynomial) function + Int_t fGFitBad; // if != 0 fit failed + Int_t fGDegree; // degree of polynomial to be fitted to the background + Int_t fGNdf; // number of degrees of freedom of polynomial fit + Double_t fGChisq; // chi squared of polynomial fit + Double_t fGProb; // chi squared probability ofg polynomial fit + + TArrayD fGValues; + TArrayD fGErrors; + + const static Int_t fGNdim = 9; + Double_t fGEmat[fGNdim][fGNdim]; + Double_t fGEma[fGNdim][fGNdim]; + Double_t fGCorr[fGNdim][fGNdim]; + + Int_t fGMbins; // number of bins in the fit range + Int_t fGNzero; // numnber of bins in the fit range with zero entry + Int_t fGIstat; + + //-------------------- + + static const TString gsDefName; //! Default Name + static const TString gsDefTitle; //! Default Title + + Bool_t DetExcess(); + Bool_t FitPolynomial(); + Bool_t FitGaussPoly(); + Bool_t RebinHistogram(Double_t x0, Int_t nrebin); + +public: + MHFindSignificance(const char *name=NULL, const char *title=NULL); + ~MHFindSignificance(); + + Bool_t FindSigma(TH1 *fhist, Double_t alphamin, Double_t alphamax, + Int_t degree, Double_t alphasig, + Bool_t drawpoly, Bool_t fitgauss, Bool_t print); + + Bool_t SigmaLiMa(Double_t non, Double_t noff, Double_t gamma, + Double_t *siglima); + + Bool_t SigmaVsAlpha(TH1 *fhist, Double_t alphamin, Double_t alphamax, + Int_t degree, Bool_t print); + + Double_t GetSignificance() { return fSigLiMa; } + + Bool_t DrawFit(Option_t *opt=NULL); + + Float_t GetDegree() const { return fDegree; } + Float_t GetProb() const { return fProb; } + Float_t GetNdf() const { return fNdf; } + Float_t GetGamma() const { return fGamma; } + Float_t GetNon() const { return fNon; } + Float_t GetNex() const { return fNex; } + Float_t GetNbg() const { return fNbg; } + Float_t GetSigLiMa() const { return fSigLiMa; } + Float_t GetMbins() const { return fMbins; } + Float_t GetAlphasi() const { return fAlphasi; } + + void SetRebin(Bool_t b=kTRUE); + void SetReduceDegree(Bool_t b=kTRUE); + + void PrintPoly(Option_t *opt=NULL); + void PrintPolyGauss(Option_t *opt=NULL); + + ClassDef(MHFindSignificance, 1) // Determine significance from alpha plot +}; + +#endif + + + + + + + + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/mhist/MHMatrix.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 2299) +++ /trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 2300) @@ -635,5 +635,7 @@ MEvtLoop evtloop; evtloop.SetParList(plist); + evtloop.SetName("EvtLoopMatrix"); //evtloop.SetProgressBar(&bar); + if (!evtloop.Eventloop()) Index: /trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mhist/Makefile (revision 2299) +++ /trunk/MagicSoft/Mars/mhist/Makefile (revision 2300) @@ -60,7 +60,10 @@ MHSigmaTheta.cc \ MHTriggerLvl0.cc \ + MHOnSubtraction.cc \ + MHFindSignificance.cc \ + MHCT1Supercuts.cc \ MHCamera.cc # MHCurrents.cc \ -# MHOnSubtraction.cc \ + SRCS = $(SRCFILES)