Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2298) +++ 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 2298) +++ 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 2298) +++ 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 2298) +++ 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 \