Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMyFindSuperCuts.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMyFindSuperCuts.cc (revision 6310) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMyFindSuperCuts.cc (revision 6310) @@ -0,0 +1,1475 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 6/2003 +! Author(s): Wolfgang Wittek, 7/2003 +! Author(s): Marcos Lopez, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MMyFindSuperCuts // +// // +// Class for otimizing the parameters of the supercuts // +// // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MMyFindSuperCuts.h" + +#include // fabs + +#include +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MContinue.h" +#include "MMySuperCuts.h" +#include "MMySuperCutsCalc.h" +#include "MDataElement.h" +#include "MDataMember.h" + +#include "MEvtLoop.h" +#include "MFSelFinal.h" +#include "MF.h" +#include "MFEventSelector.h" +#include "MFEventSelector2.h" +#include "MFillH.h" +//#include "MGeomCamCT1Daniel.h" +#include "MFEventSelector.h" +#include "MGeomCamMagic.h" +#include "MH3.h" +#include "MHSupercuts.h" +#include "MHFindSignificance.h" +#include "MHMatrix.h" +#include "MHOnSubtraction.h" + +#include "MLog.h" +#include "MLogManip.h" +#include "MMatrixLoop.h" +#include "MMinuitInterface.h" +#include "MParList.h" +#include "MProgressBar.h" +#include "MReadMarsFile.h" +#include "MReadTree.h" +#include "MTaskList.h" + +#include "TVector2.h" +#include "MSrcPosCam.h" +#include "MHillasSrcCalc.h" +#include "MSrcPosCalc.h" +#include "MObservatory.h" + +ClassImp(MMyFindSuperCuts); + +using namespace std; + + +//------------------------------------------------------------------------ +// +// fcnSupercuts +// +// - calculates the quantity to be minimized (using TMinuit) +// +// - the quantity to be minimized is (-1)*significance of the gamma signal +// in the alpha distribution (after cuts) +// +// - the parameters to be varied in the minimization are the cut parameters +// (par) +// +// ------------------------- +// - Read the Matrix +// - Fill a histogram with the alpha (6th column) +// - Fit the alpha plot +// - Get result of fit (significance) +// +static void fcnSuperCuts(Int_t &npar, Double_t *gin, Double_t &f, + Double_t *par, Int_t iflag) +{ + + // + // Set by hand the paramters range + // + if ( par[0]<0 || par[0]>1.8 || par[4]>par[0] || par[4]<0 || par[0]1 || par[3]>1 || par[4]>1 || par[1]>1 || par[6]>1 || par[7]>1 || + par[8]<0 || par[8]>1.8 || par[12]>par[0] || par[12]<0 || par[8]1. || par[14]<0 || par[14]>1 || par[15]<0 || par[15]>1 || + par[16]<0 || par[16]>1.8 || par[20]<0 || par[20]>1.8 || par[16]0.18 || + par[24]>200 || par[28]<-200 || + par[32]>1 || par[36]<0 || + par[40]>1 + ) + { + gLog << "Parameter out of limit" << endl; + f=1e10; + cout << f << endl; + return; + } + + //------------------------------------------------------------- + // save pointer to the MINUIT object for optimizing the supercuts + // because it will be overwritten + // when fitting the alpha distribution in MHFindSignificance + TMinuit *savePointer = gMinuit; + //------------------------------------------------------------- + + + MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit(); + + MParList *plistfcn = (MParList*)evtloopfcn->GetParList(); + + MMySuperCuts *super = (MMySuperCuts*)plistfcn->FindObject("MMySuperCuts"); + if (!super) + { + gLog << "fcnSupercuts : MMySuperCuts object '" << "MMySuperCuts" + << "' not found... aborting" << endl; + return; + } + + // + // transfer current parameter values to MMySuperCuts + // + // Attention : npar is the number of variable parameters + // not the total number of parameters + // + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); + + super->SetParameters(TArrayD(fNparx, par)); + + + //$$$$$$$$$$$$$$$$$$$$$ + // for testing + gLog.Separator("Current params"); + TArrayD checkparameters = super->GetParameters(); + gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", " + << fNparx << endl; + // gLog << "fcnsupercuts : i, par, checkparameters =" << endl; + // for (Int_t i=0; iPrint(); + + //$$$$$$$$$$$$$$$$$$$$$ + + + // + // EventLoop: Fill the alpha plot with these cuts + // + gLog.SetNullOutput(kTRUE); + evtloopfcn->Eventloop(); + + //tasklistfcn->PrintStatistics(0, kTRUE); + + + + // + // Calculate significance + // + MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3"); + if (!alpha) + return; + + TH1 &alphaHist = alpha->GetHist(); + alphaHist.SetName("alpha-fcnSupercuts"); + + //------------------------------------------- + // set Minuit pointer to zero in order not to destroy the TMinuit + // object for optimizing the supercuts + gMinuit = NULL; + + + //================================================================= + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + // + const Double_t alphasig = 10.0; + const Double_t alphamin = 30.0; + const Double_t alphamax = 90.0; + const Int_t degree = 0; + + Bool_t drawpoly; + Bool_t fitgauss; + TCanvas* c; + if (iflag == 3) + { + c= new TCanvas(); + drawpoly = kTRUE; + fitgauss = kTRUE; + } + else + { + drawpoly = kFALSE; + fitgauss = kFALSE; + } + fitgauss = kTRUE; + const Bool_t print = kFALSE; + + MHFindSignificance findsig; + //findsig.SetRebin(kTRUE); + findsig.SetRebin(kFALSE); + findsig.SetReduceDegree(kFALSE); + + + + const Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + gLog.SetNullOutput(kFALSE); + + //================================================================= + + // reset gMinuit to the MINUIT object for optimizing the supercuts + gMinuit = savePointer; + //------------------------------------------- + + if (!rc) + { + gLog << "fcnSupercuts : FindSigma() failed" << endl; + f = 1.e10; + return; + } + + // plot some quantities during the optimization + MHSupercuts *plotsuper = (MHSupercuts*)plistfcn->FindObject("MHSupercuts"); + if (plotsuper) + plotsuper->Fill(&findsig); + + + // + // Print info + // + gLog << "(Non, Nbg, Nex, S/N, sigma, Significane) = " + << findsig.GetNon() << " " + << findsig.GetNbg()<< " " + << findsig.GetNex() << " " + << findsig.GetNex()/ findsig.GetNbg()<< " " + << findsig.GetSigmaGauss() << " " + << findsig.GetSignificance() << endl; + + const Double_t sigmagauss = findsig.GetSigmaGauss(); + // if (sigmagauss>10) +// { +// gLog << "fcnSupercuts : Alpha plot too width = " << sigmagauss << endl; +// f = 1.e10; +// return; +// } + + + // + // Get significance + // + const Double_t significance = findsig.GetSignificance(); + //const Double_t ratio = findsig.GetNex()/findsig.GetNbg(); + + // f = significance>0 ? -significance : 0; + // f = -TMath::Abs(significance)/sqrt(sigmagauss); + // f = - significance*sqrt(findsig.GetNex())/sqrt(sigmagauss); + //f = - significance*significance/sqrt(sigmagauss); + + f = - significance*findsig.GetProb(); + cout << f << endl; + +// // +// // optimize signal/background ratio +// Double_t ratio = findsig.GetNbg()>0.0 ? +// findsig.GetNex()/findsig.GetNbg() : 0.0; + +// //f = -ratio; +// if(significance < 5) +// f=1e10; + +// f=f*ratio*ratio; +// cout << f << " " << significance << " " << ratio << endl; + +} + + + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MMyFindSuperCuts::MMyFindSuperCuts(const char *name, const char *title) + : fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL),fSizeCutLow(2000.0),fSizeCutUp(10000000), fOptimizationMode(0) +{ + fName = name ? name : "MMyFindSuperCuts"; + fTitle = title ? title : "Optimizer of the supercuts"; + + //--------------------------- + // camera geometry is needed for conversion mm ==> degree + //fCam = new MGeomCamCT1Daniel; + fCam = new MGeomCamMagic; + + // matrices to contain the training/test samples + fMatrixTrain = new MHMatrix("MatrixTrain"); + fMatrixTest = new MHMatrix("MatrixTest"); + + // objects of MMySuperCutsCalc to which these matrices are attached + fCalcHadTrain = new MMySuperCutsCalc("SupercutsCalcTrain"); + fCalcHadTest = new MMySuperCutsCalc("SupercutsCalcTest"); + + + // Define columns of matrices + *fLog << inf << "Init Mapping of Train Matrix" << endl; + fCalcHadTrain->InitMapping(fMatrixTrain); + *fLog << inf << "Init Mapping of Test Matrix" << endl; + fCalcHadTest->InitMapping(fMatrixTest); +} + + +// -------------------------------------------------------------------------- +// +// Default destructor. +// +MMyFindSuperCuts::~MMyFindSuperCuts() +{ + delete fCam; + delete fMatrixTrain; + delete fMatrixTest; + delete fCalcHadTrain; + delete fCalcHadTest; +} + + + +// -------------------------------------------------------------------------- +// +// Define the matrix 'fMatrixTrain' for the training sample +// +// alltogether 'howmanytrain' events are read from file 'nametrain'; +// the events are selected according to a target distribution 'hreftrain' +// +// +Bool_t MMyFindSuperCuts::DefineTrainMatrix( + const TString &filenametrain, MH3 &hreftrain, + const Int_t howmanytrain, const TString &outfiletrain) +{ + if (filenametrain.IsNull() || howmanytrain <= 0) + return kFALSE; + + *fLog << endl + << "***************************************************************" + << endl + << " Filling training matrix from file: '" << filenametrain << "'" + << endl + << " Selecting: " << howmanytrain << " events"; + + if (!hreftrain.GetHist().GetEntries()==0) + { + *fLog << ", according to a distribution given by the MH3 object '" + << hreftrain.GetName() << "'" << endl; + } + else + { + *fLog << ", Randomly" << endl; + } + + *fLog << "***************************************************************" + << endl; + + + // + // ParList + // + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + + // entries in MParList + plist.AddToList(fCam); + plist.AddToList(fMatrixTrain); + + + // + // TaskList + // + MReadTree read("Events", filenametrain); + //MReadTree read("Parameters", filenametrain); + read.DisableAutoScheme(); + + // MFEventSelector2 seltrain(hreftrain); +// seltrain.SetNumMax(howmanytrain); +// seltrain.SetName("selectTrain"); + + // Random selection of events + MFEventSelector seltrain; + seltrain.SetNumSelectEvts(howmanytrain); + seltrain.SetName("selectTrain"); + + MFillH filltrain(fMatrixTrain); + filltrain.SetFilter(&seltrain); + filltrain.SetName("fillMatrixTrain"); + +// // ReCalculate hillass +// MObservatory obs; +// plist.AddToList(&obs); + +// MSrcPosCalc srccalc; +// srccalc.SetPositionXY(3.46945e-17, 0.0487805); +// MHillasSrcCalc hsrc; + + + + // entries in MTaskList + tlist.AddToList(&read); + tlist.AddToList(&seltrain); + // tlist.AddToList(&srccalc); +// tlist.AddToList(&hsrc); + tlist.AddToList(&filltrain); + + + // + // EventLoop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTrain"); + evtloop.SetProgressBar(&bar); + + if (!evtloop.Eventloop()) + return kFALSE; + + // if (!evtloop.PreProcess()) +// return kFALSE; + +// Int_t i=0; +// while(tlist.Process()) +// { +// i++; +// if (i>howmanytrain) +// break; +// } + +// evtloop.PostProcess(); + + tlist.PrintStatistics(0, kTRUE); + + + // + // Print + // + fMatrixTrain->Print("SizeCols"); + Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows(); + if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + { + *fLog << "MMyFindSuperCuts::DefineTrainMatrix; no.of generated events (" + << howmanygenerated + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + *fLog << "training matrix was filled" << endl; + *fLog << "=============================================" << endl; + + + + // + // Write out training matrix + // + if (outfiletrain != "") + { + TFile filetr(outfiletrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MMyFindSuperCuts::DefineTrainMatrix; Training matrix was written onto file '" + << outfiletrain << "'" << endl; + } + + // fMatrixTrain->Print("data"); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Define the matrix for the test sample +// +// alltogether 'howmanytest' events are read from file 'nametest' +// +// the events are selected according to a target distribution 'hreftest' +// +// +Bool_t MMyFindSuperCuts::DefineTestMatrix( + const TString &nametest, MH3 &hreftest, + const Int_t howmanytest, const TString &filetest) +{ + if (nametest.IsNull() || howmanytest<=0) + return kFALSE; + + *fLog << "=============================================" << endl; + *fLog << "fill test matrix from file '" << nametest + << "', select " << howmanytest + << " events " << endl; + if (!hreftest.GetHist().GetEntries()==0) + { + *fLog << " according to a distribution given by the MH3 object '" + << hreftest.GetName() << "'" << endl; + } + else + { + *fLog << " randomly" << endl; + } + + + // + // ParList + // + MParList plist; + MTaskList tlist; + + // entries in MParList + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTest); + + + // + // TaskList + // + MReadTree read("Parameters", nametest); + read.DisableAutoScheme(); + + // MFEventSelector2 seltest(hreftest); +// seltest.SetNumMax(howmanytest); +// seltest.SetName("selectTest"); + + MFillH filltest(fMatrixTest); + //filltest.SetFilter(&seltest); + + // entries in MTaskList + tlist.AddToList(&read); + //tlist.AddToList(&seltest); + tlist.AddToList(&filltest); + + + // + // EventLoop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTest"); + evtloop.SetProgressBar(&bar); + + + if (!evtloop.PreProcess()) + return kFALSE; + + Int_t i=0; + while(tlist.Process()) + { + i++; + if (i>100000) + break; + } + + evtloop.PostProcess(); + + +// if (!evtloop.Eventloop()) +// return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + + // + // Print + // + fMatrixTest->Print("SizeCols"); + const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows(); + if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest)) + { + *fLog << "MMyFindSuperCuts::DefineTestMatrix; no.of generated events (" + << howmanygenerated + << ") is incompatible with the no.of requested events (" + << howmanytest << ")" << endl; + } + + *fLog << "test matrix was filled" << endl; + *fLog << "=============================================" << endl; + + + + // + // Write out test matrix + // + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MMyFindSuperCuts::DefineTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// Define the matrices for the training and test sample respectively +// +// +// +Bool_t MMyFindSuperCuts::DefineTrainTestMatrix( + const TString &name, MH3 &href, + const Int_t howmanytrain, const Int_t howmanytest, + const TString &filetrain, const TString &filetest) +{ + *fLog << "=============================================" << endl; + *fLog << "fill training and test matrix from file '" << name + << "', select " << howmanytrain + << " training and " << howmanytest << " test events " << endl; + if (!href.GetHist().GetEntries()==0) + { + *fLog << " according to a distribution given by the MH3 object '" + << href.GetName() << "'" << endl; + } + else + { + *fLog << " randomly" << endl; + } + + + // + // ParList + // + MParList plist; + MTaskList tlist; + + // entries in MParList + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTrain); + plist.AddToList(fMatrixTest); + + + + // + // TaskList + // + MReadMarsFile read("Events", name); + read.DisableAutoScheme(); + + MFEventSelector2 selector(href); + selector.SetNumMax(howmanytrain+howmanytest); + selector.SetName("selectTrainTest"); + selector.SetInverted(); + + MContinue cont(&selector); + cont.SetName("ContTrainTest"); + + Double_t prob = ( (Double_t) howmanytrain ) + / ( (Double_t)(howmanytrain+howmanytest) ); + MFEventSelector split; + split.SetSelectionRatio(prob); + + MFillH filltrain(fMatrixTrain); + filltrain.SetFilter(&split); + filltrain.SetName("fillMatrixTrain"); + + + // consider this event as candidate for a test event + // only if event was not accepted as a training event + + MContinue conttrain(&split); + conttrain.SetName("ContTrain"); + + MFillH filltest(fMatrixTest); + filltest.SetName("fillMatrixTest"); + + + // entries in MTaskList + tlist.AddToList(&read); + tlist.AddToList(&cont); + tlist.AddToList(&split); + tlist.AddToList(&filltrain); + tlist.AddToList(&conttrain); + tlist.AddToList(&filltest); + + + // + // EventLoop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTrainTest"); + evtloop.SetProgressBar(&bar); + + Int_t maxev = -1; + if (!evtloop.Eventloop(maxev)) + return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + fMatrixTrain->Print("SizeCols"); + const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows(); + if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + { + *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; no.of generated events (" + << generatedtrain + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + fMatrixTest->Print("SizeCols"); + const Int_t generatedtest = fMatrixTest->GetM().GetNrows(); + if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest)) + { + *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; no.of generated events (" + << generatedtest + << ") is incompatible with the no.of requested events (" + << howmanytest << ")" << endl; + } + + + *fLog << "training and test matrix were filled" << endl; + *fLog << "=============================================" << endl; + + + + // + // Write out training matrix + // + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + // + // Write out test matrix + // + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// Read training and test matrices from files +// +Bool_t MMyFindSuperCuts::ReadMatrix(const TString &filetrain, const TString &filetest) +{ + // + // Read in training matrix + // + TFile filetr(filetrain); + fMatrixTrain->Read("MatrixTrain"); + fMatrixTrain->Print("SizeCols"); + + *fLog << "MMyFindSuperCuts::ReadMatrix; Training matrix was read in from file '" << filetrain << "'" << endl; + filetr.Close(); + + + // + // Read in test matrix + // + TFile filete(filetest); + fMatrixTest->Read("MatrixTest"); + fMatrixTest->Print("SizeCols"); + + *fLog <<"MMyFindSuperCuts::ReadMatrix; Test matrix was read in from file '" + << filetest << "'" << endl; + filete.Close(); + + + return kTRUE; +} + + + +//------------------------------------------------------------------------ +// +// Steering program for optimizing the supercuts +// --------------------------------------------- +// +// the criterion for the 'optimum' is +// +// - a maximum significance of the gamma signal in the alpha plot, +// in which the supercuts have been applied +// +// The various steps are : +// +// - setup the event loop to be executed for each call to fcnSupercuts +// - call TMinuit to do the minimization of (-significance) : +// the fcnSupercuts function calculates the significance +// for the current cut values +// for this - the alpha plot is produced in the event loop, +// in which the cuts have been applied +// - the significance of the gamma signal in the alpha plot +// is calculated +// +// Needed as input : (to be set by the Set functions) +// +// - fFilenameParam name of file to which optimum values of the +// parameters are written +// +// - fHadronnessName name of container where MMySuperCutsCalc will +// put the supercuts hadronness +// +// - for the minimization, the starting values of the parameters are taken +// - from file parSCinit (if it is != "") +// - or from the arrays params and/or steps +// +//---------------------------------------------------------------------- +Bool_t MMyFindSuperCuts::FindParams(TString parSCinit, + TArrayD ¶ms, TArrayD &steps ) +{ + + fCalcHadTrain->SetSizeCuts(fSizeCutLow,fSizeCutUp); + fCalcHadTest->SetSizeCuts(fSizeCutLow,fSizeCutUp); + + + // Setup the event loop which will be executed in the + // fcnSupercuts function of MINUIT + // + // parSCinit is the name of the file containing the initial values + // of the parameters; + // if parSCinit = "" 'params' and 'steps' are taken as initial values + + *fLog << "=============================================" << endl; + *fLog << "Setup event loop for fcnSupercuts" << endl; + + + // + // Check that all the necessary containers already exists + // + if (fHadronnessName.IsNull()) + { + *fLog << "MyFindSuperCuts::FindParams; hadronness name is not defined... aborting" + << endl; + return kFALSE; + } + + if (fMatrixTrain == NULL) + { + *fLog << "MMyFindSuperCuts::FindParams; training matrix is not defined... aborting" + << endl; + return kFALSE; + } + + if (fMatrixTrain->GetM().GetNrows() <= 0) + { + *fLog << "MMyFindSuperCuts::FindParams; training matrix has no entries" + << endl; + return kFALSE; + } + + + + // + // ParList + // + MParList parlistfcn; + MTaskList tasklistfcn; + parlistfcn.AddToList(&tasklistfcn); + + + + // + // create container for the supercut parameters + // and set them to their initial values + // + MMySuperCuts super; + + // take initial values from file parSCinit + if (parSCinit != "") + { + TFile inparam(parSCinit); + super.Read("MMySuperCuts"); + inparam.Close(); + *fLog << "MMyFindSuperCuts::FindParams; initial values of parameters are taken from file " << parSCinit << endl; + super.Print(); + } + + // take initial values from 'params' and/or 'steps' + //else if (params.GetSize() != 0 || steps.GetSize() != 0 ) + if (params.GetSize() != 0 || steps.GetSize() != 0 ) + { + if (params.GetSize() != 0) + { + *fLog << "MMyFindSuperCuts::FindParams; initial values of parameters are taken from 'params'" << endl; + super.SetParameters(params); + } + if (steps.GetSize() != 0) + { + *fLog << "MMyFindSuperCuts::FindParams; initial step sizes are taken from 'steps'" + << endl; + super.SetStepsizes(steps); + } + } + else + { + *fLog << "MMyFindSuperCuts::FindParams; initial values and step sizes are taken from the MMySuperCuts constructor" + << endl; + } + + + // Alpha binning + const TString mh3Name = "AlphaFcn"; + MBinning binsalpha("Binning"+mh3Name); + // binsalpha.SetEdges(54, -12.0, 96.0); + binsalpha.SetEdges(45, 0.0, 90.0); + + // Alpha plot |alpha| + MH3 alpha("MatrixTrain[6]"); // |alpha| assumed to be in column 6 + alpha.SetName(mh3Name); + + *fLog << warn << "WARNING----->ALPHA IS ASSUMED TO BE IN COLUMN 6!!!!!!!!" + << endl; + + + // book histograms to be filled during the optimization + // ! not in the event loop ! + MHSupercuts plotsuper; + plotsuper.SetupFill(&parlistfcn); + + + // Add to ParList + parlistfcn.AddToList(&super); + parlistfcn.AddToList(fCam); + parlistfcn.AddToList(fMatrixTrain); + parlistfcn.AddToList(&binsalpha); + parlistfcn.AddToList(&alpha); + parlistfcn.AddToList(&plotsuper); + + + + // + // MTaskList + // + MMatrixLoop loop(fMatrixTrain); // loop over rows of matrix + + // calculate supercuts hadronness + //fCalcHadTrain->SetHadronnessName(fHadronnessName); + + // apply the supercuts + MF scfilter(fHadronnessName+".fHadronness>0.5"); + MContinue supercuts(&scfilter); + + // Fill Alpha Plot + MFillH fillalpha(&alpha); + + + // Add to TaskList + tasklistfcn.AddToList(&loop); + tasklistfcn.AddToList(fCalcHadTrain); + tasklistfcn.AddToList(&supercuts); + tasklistfcn.AddToList(&fillalpha); + + + + // + // EventLoop + // + MEvtLoop evtloopfcn("EvtLoopFCN"); + evtloopfcn.SetParList(&parlistfcn); + *fLog << "Event loop for fcnSupercuts has been setup" << endl; + + + + + //----------------------------------------------------------------------- + // + //---------- Start of minimization part -------------------- + // + // Do the minimization with MINUIT + // + // Be careful: This is not thread safe + // + *fLog << "========================================================"<< endl; + *fLog << "Start minimization for supercuts" << endl; + + + // ------------------------------------------- + // prepare call to MINUIT + // + + // get initial values of parameters + fVinit = super.GetParameters(); + fStep = super.GetStepsizes(); + + TString name[fVinit.GetSize()]; + fStep.Set(fVinit.GetSize()); + fLimlo.Set(fVinit.GetSize()); + fLimup.Set(fVinit.GetSize()); + fFix.Set(fVinit.GetSize()); + + fNpar = fVinit.GetSize(); + + + + for (UInt_t i=0; i= 48) + // { + // fStep[i] = 0.0; + // fFix[i] = 1; + // } + // + + // //} +// // vary only first 23 parameters +// // +// for (UInt_t i=0; i= 24) +// { +// fStep[i] = 0.0; +// fFix[i] = 1; +// } +// } + + + // ------------------------------------------- + // call MINUIT + + TStopwatch clock; + clock.Start(); + + *fLog << "before calling CallMinuit" << endl; + + MMinuitInterface inter; + Bool_t rc = inter.CallMinuit(fcnSuperCuts, name, + fVinit, fStep, fLimlo, fLimup, fFix, + &evtloopfcn, "SIMPLEX", kFALSE); + + *fLog << "after calling CallMinuit" << endl; + + *fLog << "Time spent for the minimization in MINUIT : " << endl;; + clock.Stop(); + clock.Print(); + + plotsuper.DrawClone(); + + if (!rc) + return kFALSE; + + *fLog << "Minimization for supercuts finished" << endl; + *fLog << "========================================================" << endl; + + + + // ----------------------------------------------------------------- + // in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put + // into MMySuperCuts + + // + // Write optimum parameter values onto file filenameParam + // + TFile outparam(fFilenameParam, "RECREATE"); + super.Write(); + outparam.Close(); + + *fLog << "Optimum parameter values for supercuts were written onto file '" + << fFilenameParam << "' :" << endl; + + const TArrayD &check = super.GetParameters(); + for (Int_t i=0; iGetM().GetNrows() <= 0) + { + *fLog << "MMyFindSuperCuts::TestParams; test matrix has no entries" + << endl; + return kFALSE; + } + + // ------------- TEST the supercuts ------------------------------ + // + *fLog << "Test the supercuts on the test sample" << endl; + + // ----------------------------------------------------------------- + // read optimum parameter values from file filenameParam + // into array 'supercutsPar' + + TFile inparam(fFilenameParam); + MMySuperCuts scin; + scin.Read("MMySuperCuts"); + inparam.Close(); + + *fLog << "Optimum parameter values for supercuts were read from file '"; + *fLog << fFilenameParam << "' :" << endl; + + const TArrayD &supercutsPar = scin.GetParameters(); + for (Int_t i=0; iSetHadronnessName(fHadronnessName); + + + //------------------------------------------- + + MMatrixLoop loop(fMatrixTest); + + // plot alpha before applying the supercuts + const TString mh3NameB = "AlphaBefore"; + MBinning binsalphabef("Binning"+mh3NameB); + binsalphabef.SetEdges(54, -12.0, 96.0); + + // |alpha| is assumed to be in column 7 of the matrix + MH3 alphabefore("MatrixTest[6]"); + alphabefore.SetName(mh3NameB); + + TH1 &alphahistb = alphabefore.GetHist(); + alphahistb.SetName("alphaBefore-TestParams"); + + MFillH fillalphabefore(&alphabefore); + fillalphabefore.SetName("FillAlphaBefore"); + + // apply the supercuts + MF scfilter(fHadronnessName+".fHadronness>0.5"); + MContinue applysupercuts(&scfilter); + + // plot alpha after applying the supercuts + const TString mh3NameA = "AlphaAfter"; + MBinning binsalphaaft("Binning"+mh3NameA); + binsalphaaft.SetEdges(54, -12.0, 96.0); + + MH3 alphaafter("MatrixTest[6]"); + alphaafter.SetName(mh3NameA); + + TH1 &alphahista = alphaafter.GetHist(); + alphahista.SetName("alphaAfter-TestParams"); + + MFillH fillalphaafter(&alphaafter); + fillalphaafter.SetName("FillAlphaAfter"); + + //****************************** + // entries in MParList + + parlist2.AddToList(&tasklist2); + + parlist2.AddToList(&supercuts); + + parlist2.AddToList(fCam); + parlist2.AddToList(fMatrixTest); + + parlist2.AddToList(&binsalphabef); + parlist2.AddToList(&alphabefore); + + parlist2.AddToList(&binsalphaaft); + parlist2.AddToList(&alphaafter); + + //****************************** + // entries in MTaskList + + tasklist2.AddToList(&loop); + tasklist2.AddToList(&fillalphabefore); + + tasklist2.AddToList(fCalcHadTest); + tasklist2.AddToList(&applysupercuts); + + tasklist2.AddToList(&fillalphaafter); + + //****************************** + + MProgressBar bar2; + MEvtLoop evtloop2; + evtloop2.SetParList(&parlist2); + evtloop2.SetName("EvtLoopTestParams"); + evtloop2.SetProgressBar(&bar2); + Int_t maxevents2 = -1; + if (!evtloop2.Eventloop(maxevents2)) + return kFALSE; + + tasklist2.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // draw the alpha plots + + MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3"); + TH1 &alphaHist1 = alphaBefore->GetHist(); + alphaHist1.SetXTitle("|\\alpha| [\\circ]"); + + MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3"); + TH1 &alphaHist2 = alphaAfter->GetHist(); + alphaHist2.SetXTitle("|\\alpha| [\\circ]"); + alphaHist2.SetName("alpha-TestParams"); + + TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300); + c->Divide(2,1); + + gROOT->SetSelectedPad(NULL); + + c->cd(1); + alphaHist1.DrawCopy(); + + c->cd(2); + alphaHist2.DrawCopy(); + + + + //------------------------------------------- + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + const Double_t alphasig = 20.0; + const Double_t alphamin = 30.0; + const Double_t alphamax = 90.0; + const Int_t degree = 2; + const Bool_t drawpoly = kTRUE; + const Bool_t fitgauss = kTRUE; + const Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + + const Double_t significance = findsig.GetSignificance(); + const Double_t alphasi = findsig.GetAlphasi(); + + *fLog << "Significance of gamma signal after supercuts in test sample : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + //------------------------------------------- + + + *fLog << "Test of supercuts part finished" << endl; + *fLog << "======================================================" << endl; + + return kTRUE; +} + Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMyFindSuperCuts.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMyFindSuperCuts.h (revision 6310) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMyFindSuperCuts.h (revision 6310) @@ -0,0 +1,149 @@ +#ifndef MARS_MMyFindSuperCuts +#define MARS_MMyFindSuperCuts + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif +#ifndef ROOT_TArrayI +#include +#endif + +class MFilter; +class MEvtLoop; +class MH3; +class MMySuperCutsCalc; +class MGeomCam; +class MHMatrix; +/* +#include "MFilter.h" +#include "MEvtLoop.h" +#include "MH3.h" +#include "MMySuperCutsCalc.h" +#include "MGeomCam.h" +#include "MHMatrix.h" +*/ + +class MMyFindSuperCuts : public MParContainer +{ +private: + + TString fFilenameTrain; + TString fFilenameTest; + + Int_t fHowManyTrain; + Int_t fHowManyTest; + + Bool_t fUseOrigDistribution; + + TString fFilenameParam; + + TString fHadronnessName; + + MMySuperCutsCalc *fCalcHadTrain; + MMySuperCutsCalc *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]; + TArrayD fVinit; + TArrayD fStep; + TArrayD fLimlo; + TArrayD fLimup; + TArrayI fFix; + + UInt_t fNpar; + + TString fMethod; + + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + Int_t fErrMinimize; + //-------------------------------------------- + + Float_t fSizeCutLow; + Float_t fSizeCutUp; + + Int_t fOptimizationMode; + +public: + + MMyFindSuperCuts(const char *name=NULL, const char *title=NULL); + ~MMyFindSuperCuts(); + + void SetSizeCutLow(Float_t size) {fSizeCutLow = size;} + void SetSizeCuts(Float_t sizeLow, Float_t sizeUp) + { + fSizeCutLow = sizeLow; + fSizeCutUp = sizeUp; + } + + + + 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, MH3 &href, + const Int_t howmany, const TString &filetrain); + + Bool_t DefineTestMatrix(const TString &name, MH3 &href, + const Int_t howmany, const TString &filetest); + + Bool_t DefineTrainTestMatrix(const TString &name, MH3 &href, + const Int_t howmanytrain, const Int_t howmanytest, + const TString &filetrain, const TString &filetest); + + MHMatrix *GetMatrixTrain() { return fMatrixTrain; } + MHMatrix *GetMatrixTest() { return fMatrixTest; } + + Bool_t ReadMatrix( const TString &filetrain, const TString &filetest); + + Bool_t FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps); + Bool_t TestParams(); + + + void SetOptimizationMode(Int_t mode ) { fOptimizationMode=mode; } + + ClassDef(MMyFindSuperCuts, 1) // Class for optimization of the Supercuts +}; + +#endif + + + + + + + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCuts.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCuts.cc (revision 6310) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCuts.cc (revision 6310) @@ -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, 08/2003 +! Author(s): Thomas Bretz, 08/2003 +! Author(s): Marcos Lopez, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MMySuperCuts // +// // +// this is the container for the parameters of the supercuts // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MMySuperCuts.h" + +#include +#include + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MMySuperCuts); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// constructor +// +// MMySuperCuts::MMySuperCuts(const char *name, const char *title) +// : fParameters(18), fStepsizes(18), +// fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+3), +// fWidthUp(fParameters.GetArray()+6), fWidthLo(fParameters.GetArray()+9), +// fDistUp(fParameters.GetArray()+12), fDistLo(fParameters.GetArray()+15) +// { +MMySuperCuts::MMySuperCuts(const char *name, const char *title) + : fParameters(48), fStepsizes(48), + fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+4), + fWidthUp(fParameters.GetArray()+4*2), fWidthLo(fParameters.GetArray()+4*3), + fDistUp(fParameters.GetArray()+4*4), fDistLo(fParameters.GetArray()+4*5), + fAsymUp(fParameters.GetArray()+4*6), fAsymLo(fParameters.GetArray()+4*7), + fConcUp(fParameters.GetArray()+4*8), fConcLo(fParameters.GetArray()+4*9), + fLeakage1Up(fParameters.GetArray()+4*10), fLeakage1Lo(fParameters.GetArray()+4*11) + +{ + + fName = name ? name : "MMySuperCuts"; + fTitle = title ? title : "Container for the supercut parameters"; + + // set supercut parameters to their default values + InitParameters(); +} + + +// -------------------------------------------------------------------------- +// +// set default values for the supercut parameters +// +void MMySuperCuts::InitParameters() +{ + //--------------------------------------------------- + // these are the default values + // + fLengthUp[0] = 0.302; + fLengthUp[1] = 0.0;//0.03; + fLengthUp[2] = 0.0; + fLengthUp[3] = 0.0; +// fLengthUp[4] = 0.0; +// fLengthUp[5] = 0.0; +// fLengthUp[6] = 0.0; +// fLengthUp[7] = 0.0; + + fLengthLo[0] = 0.052; + fLengthLo[1] = 0.; + fLengthLo[2] = 0.; + fLengthLo[3] = 0.; +// fLengthLo[4] = 0.; +// fLengthLo[5] = 0.; +// fLengthLo[6] = 0.; +// fLengthLo[7] = 0.; + + fWidthUp[0] = 0.121; + fWidthUp[1] = 0.0; + fWidthUp[2] = 0.0; + fWidthUp[3] = 0.0; +// fWidthUp[4] = 0.0; +// fWidthUp[5] = 0.0; +// fWidthUp[6] = 0.0; +// fWidthUp[7] = 0.0; + + fWidthLo[0] = 0.057; + fWidthLo[1] = 0.; + fWidthLo[2] = 0.; + fWidthLo[3] = 0.; +// fWidthLo[4] = 0.; +// fWidthLo[5] = 0.; +// fWidthLo[6] = 0.; +// fWidthLo[7] = 0.; + + fDistUp[0] = 1.; + fDistUp[1] = 0.0; + fDistUp[2] = 0.0; + fDistUp[3] = 0.0; +// fDistUp[4] = 0.0; +// fDistUp[5] = 0.0; +// fDistUp[6] = 0.0; +// fDistUp[7] = 0.0; + + fDistLo[0] = 0.3; + fDistLo[1] = 0.0; + fDistLo[2] = 0.0; + fDistLo[3] = 0.0; +// fDistLo[4] = 0.0; +// fDistLo[5] = 0.0; +// fDistLo[6] = 0.0; +// fDistLo[7] = 0.0; + + // dummy values + + fAsymUp[0] = 100; + fAsymUp[1] = 0.0; + fAsymUp[2] = 0.0; + fAsymUp[3] = 0.0; + + fAsymLo[0] = -100; + fAsymLo[1] = 0.0; + fAsymLo[2] = 0.0; + fAsymLo[3] = 0.0; + + + fConcUp[0] = 1.; + fConcUp[1] = 0.0; + fConcUp[2] = 0.0; + fConcUp[3] = 0.0; + + fConcLo[0] = 0; + fConcLo[1] = 0.0; + fConcLo[2] = 0.0; + fConcLo[3] = 0.0; + + fLeakage1Up[0] = 1.; + fLeakage1Up[1] = 0.0; + fLeakage1Up[2] = 0.0; + fLeakage1Up[3] = 0.0; + // fLeakage1Up[4] = 0.0; +// fLeakage1Up[5] = 0.0; +// fLeakage1Up[6] = 0.0; +// fLeakage1Up[7] = 0.0; + + fLeakage1Lo[0] = -1.e10; + fLeakage1Lo[1] = 0.0; + fLeakage1Lo[2] = 0.0; + fLeakage1Lo[3] = 0.0; + // fLeakage1Lo[4] = 0.0; +// fLeakage1Lo[5] = 0.0; +// fLeakage1Lo[6] = 0.0; +// fLeakage1Lo[7] = 0.0; + + + //--------------------------------------------------- + // fStepsizes + // if == 0.0 the parameter will be fixed in the minimization + // != 0.0 initial step sizes for the parameters + + // LengthUp + fStepsizes[0] = 0.01; + fStepsizes[1] = 0.01; + fStepsizes[2] = 0.1; + fStepsizes[3] = 0.1; +// fStepsizes[4] = 0.0; +// fStepsizes[5] = 0.0; +// fStepsizes[6] = 0.0; +// fStepsizes[7] = 0.0; + + // LengthLo + fStepsizes[4] = 0.01; + fStepsizes[5] = 0.01; + fStepsizes[6] = 0.1; + fStepsizes[7] = 0.1; +// fStepsizes[12] = 0.0; +// fStepsizes[13] = 0.0; +// fStepsizes[14] = 0.0; +// fStepsizes[15] = 0.0; + + // WidthUp + fStepsizes[8] = 0.01; + fStepsizes[9] = 0.01; + fStepsizes[10] = 0.1; + fStepsizes[11] = 0.1; +// fStepsizes[20] = 0.0; +// fStepsizes[21] = 0.0; +// fStepsizes[22] = 0.0; +// fStepsizes[23] = 0.0; + + // WidthLo + fStepsizes[12] = 0.01; + fStepsizes[13] = 0.01; + fStepsizes[14] = 0.1; + fStepsizes[15] = 0.1; +// fStepsizes[28] = 0.0; +// fStepsizes[29] = 0.0; +// fStepsizes[30] = 0.0; +// fStepsizes[31] = 0.0; + + // DistUp + fStepsizes[16] = 0.1; + fStepsizes[17] = 0.1; + fStepsizes[18] = 0.1; + fStepsizes[19] = 0.1; +// fStepsizes[36] = 0.0; +// fStepsizes[37] = 0.0; +// fStepsizes[38] = 0.0; +// fStepsizes[39] = 0.0; + + // DistLo + fStepsizes[20] = 0.1; + fStepsizes[21] = 0.1; + fStepsizes[22] = 0.1; + fStepsizes[23] = 0.1; +// fStepsizes[44] = 0.0; +// fStepsizes[45] = 0.0; +// fStepsizes[46] = 0.0; +// fStepsizes[47] = 0.0; + + for(int i=24;i<48;i++) + fStepsizes[i] = .1; + +} + + +// -------------------------------------------------------------------------- +// +// Set the parameter values from the array 'd' +// +// +Bool_t MMySuperCuts::SetParameters(const TArrayD &d) +{ + if (d.GetSize() != fParameters.GetSize()) + { + *fLog << err << "Sizes of d and of fParameters are different : " + << d.GetSize() << ", " << fParameters.GetSize() << endl; + return kFALSE; + } + + fParameters = d; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set the step sizes from the array 'd' +// +// +Bool_t MMySuperCuts::SetStepsizes(const TArrayD &d) +{ + if (d.GetSize() != fStepsizes.GetSize()) + { + *fLog << err << "Sizes of d and of fStepsizes are different : " + << d.GetSize() << ", " << fStepsizes.GetSize() << endl; + return kFALSE; + } + + fStepsizes = d; + + return kTRUE; +} + + + +void MMySuperCuts::Print(const Option_t* opt) const +{ + // for(Int_t i=0; i<48; i++) + // cout << i << " " << fParameters[i] << endl; + + printf("Hillas_up/lo = A + B*log(SIZE) + C*log^2(SIZE) + D*DIST^2\n"); + printf(" A B C D\n"); + printf("Length_Up [0-3]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[0],fParameters[1],fParameters[2],fParameters[3]); + printf("Length_Lo [4-7]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[4],fParameters[5],fParameters[6],fParameters[7]); + + printf("Width_Up [8-11]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[8],fParameters[9],fParameters[10],fParameters[11]); + printf("Width_Lo [12-15]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[12],fParameters[13],fParameters[14],fParameters[15]); + + printf("Dist_Up [16-19]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[16],fParameters[17],fParameters[18],fParameters[19]); + printf("Dist_Lo [20-23]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[20],fParameters[21],fParameters[22],fParameters[23]); + + printf("Asym_Up [24-27]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[24],fParameters[25],fParameters[26],fParameters[27]); + printf("Asym_Lo [28-31]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[28],fParameters[29],fParameters[30],fParameters[31]); + + printf("Conc_Up [32-35]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[32],fParameters[33],fParameters[34],fParameters[35]); + printf("Conc_Lo [36-39]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[36],fParameters[37],fParameters[38],fParameters[39]); + + printf("Leak_Up [40-43]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[40],fParameters[41],fParameters[42],fParameters[43]); + printf("Leak_Lo [44-47]: %7.3f %7.3f %7.3f %7.3f\n",fParameters[44],fParameters[45],fParameters[46],fParameters[47]); +} + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCuts.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCuts.h (revision 6310) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCuts.h (revision 6310) @@ -0,0 +1,72 @@ +#ifndef MARS_MMySuperCuts +#define MARS_MMySuperCuts + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MMySuperCuts : public MParContainer +{ + private: + + TArrayD fParameters; // supercut parameters + TArrayD fStepsizes; // step sizes of supercut parameters + + Double_t *fLengthUp; //! + Double_t *fLengthLo; //! + Double_t *fWidthUp; //! + Double_t *fWidthLo; //! + Double_t *fDistUp; //! + Double_t *fDistLo; //! + + Double_t *fAsymUp; //! + Double_t *fAsymLo; //! + Double_t *fConcUp; //! + Double_t *fConcLo; //! + Double_t *fLeakage1Up; //! + Double_t *fLeakage1Lo; //! + + + + + public: + + MMySuperCuts(const char *name=NULL, const char *title=NULL); + + void InitParameters(); + + Bool_t SetParameters(const TArrayD &d); + Bool_t SetStepsizes(const TArrayD &d); + + const TArrayD &GetParameters() const { return fParameters; } + const TArrayD &GetStepsizes() const { return fStepsizes; } + + const Double_t *GetLengthUp() const { return fLengthUp; } + const Double_t *GetLengthLo() const { return fLengthLo; } + const Double_t *GetWidthUp() const { return fWidthUp; } + const Double_t *GetWidthLo() const { return fWidthLo; } + const Double_t *GetDistUp() const { return fDistUp; } + const Double_t *GetDistLo() const { return fDistLo; } + +const Double_t *GetAsymUp() const { return fAsymUp; } + const Double_t *GetAsymLo() const { return fAsymLo; } + + const Double_t *GetConcUp() const { return fConcUp; } + const Double_t *GetConcLo() const { return fConcLo; } + const Double_t *GetLeakage1Up() const { return fLeakage1Up; } + const Double_t *GetLeakage1Lo() const { return fLeakage1Lo; } + + + void Print(const Option_t *opt ="") const; + + ClassDef(MMySuperCuts, 1) // A container for the Supercut parameters +}; + +#endif + + + Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCutsCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCutsCalc.cc (revision 6310) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCutsCalc.cc (revision 6310) @@ -0,0 +1,551 @@ +/* ======================================================================== *\ +! +! * +! * 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, 04/2003 +! Author(s): Marcos Lopez, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MMySuperCutsCalc // +// // +// this class calculates the hadronness for the supercuts // +// the parameters of the supercuts are taken // +// from the container MSupercuts // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MMySuperCutsCalc.h" + +#include +#include + +#include "TFile.h" +#include "TArrayD.h" + +#include "MParList.h" +#include "MHillasExt.h" +#include "MHillasSrc.h" +#include "MNewImagePar.h" +#include "MMcEvt.hxx" +#include "MCerPhotEvt.h" +#include "MGeomCam.h" +#include "MHadronness.h" +#include "MHMatrix.h" +#include "MMySuperCuts.h" +#include "MGeomCamMagic.h" +#include "MPointingPos.h" +#include "MNewImagePar.h" +#include "MHillasExt.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MMySuperCutsCalc); + +using namespace std; + + +// -------------------------------------------------------------------------- +// +// constructor +// +MMySuperCutsCalc::MMySuperCutsCalc(const char *hilname, const char *hilsrcname, const char *name, const char *title) + : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname), + fHilExtName("MHillasExt"), fNewParName("MNewImagePar"), + fSuperName("MMySuperCuts") , + fNoDistCut(kFALSE) ,fSizeCutLow(2000),fSizeCutUp(10000000) + +{ + fName = name ? name : "MMySuperCutsCalc"; + fTitle = title ? title : "Class to evaluate the Supercuts"; + + fMatrix = NULL; +} + + +// -------------------------------------------------------------------------- +// +// +Int_t MMySuperCutsCalc::PreProcess(MParList *pList) +{ + MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!cam) + { + *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = cam->GetConvMm2Deg(); + + fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName); + if (!fHadronness) + { + *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl; + return kFALSE; + } + + fSuper = (MMySuperCuts*)pList->FindObject(fSuperName, "MMySuperCuts"); + if (!fSuper) + { + *fLog << err << fSuperName << " [MMySuperCuts] not found... aborting." << endl; + return kFALSE; + } + + if (fMatrix) + return kTRUE; + + //----------------------------------------------------------- + fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); + if (!fHil) + { + *fLog << err << fHilName << " [MHillas] not found... aborting." << endl; + return kFALSE; + } + + fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt"); + if (!fHilExt) + { + *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl; + return kFALSE; + } + + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); + if (!fHilSrc) + { + *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; + return kFALSE; + } + + fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar"); + if (!fNewPar) + { + *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl; + return kFALSE; + } + + // fPointingPos = (MPointingPos*)pList->FindObject("MPointingPos"); +// if (!fPointingPos) +// { +// *fLog << err << " [MPointingPos] not found... aborting." << endl; +// // return kFALSE; +// } + fNewImagePar = (MNewImagePar*)pList->FindObject("MNewImagePar"); + if (!fNewImagePar) + { + *fLog << err << " [MNewImagePar] not found... aborting." << endl; + return kFALSE; + } + + fHillasExt = (MHillasExt*)pList->FindObject("MHillasExt"); + if (!fHillasExt) + { + *fLog << err << " [MHillasExt] not found... aborting." << endl; + return kFALSE; + } + +// fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); +// if (!fMcEvt) +// { +// *fLog << err << "MMcEvt not found... aborting." << endl; +// return kFALSE; +// } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculation of upper and lower limits +// +Double_t MMySuperCutsCalc::CtsMCut(const Double_t* a, Double_t ls, Double_t ls2, Double_t ct, Double_t dd2) const +{ + // define cut-function + // + // dNOMLOGSIZE = 5.0 (=log(150.0) + // dNOMCOSZA = 1.0 + // + // a: array of cut parameters + // ls: log(SIZE) - dNOMLOGSIZE + // ls2: ls^2 + // ct: Cos(ZA.) - dNOMCOSZA + // dd2: DIST^2 + +// const Double_t limit = +// a[0] + a[1] * dd2 + a[2] * ct + +// ls * (a[3] + a[4] * dd2 + a[5] * ct) + +// ls2 * (a[6] + a[7] * dd2); + +// const Double_t limit = a[0] + ls*a[1] + ls2*a[2] + ct*a[3]; + const Double_t limit = a[0] + ls*a[1] + ls2*a[2] + dd2*a[3]; + + //*fLog << "MMySuperCutsCalc::CtsMCut; *a = " + // << *a << ", " << *(a+1) << ", " << *(a+2) << ", " + // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", " + // << *(a+6) << ", " << *(a+7) << endl; + + //*fLog << "MMySuperCutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls + // << ", " << ls2 << ", " << ct << ", " << dd2 << ", " + // << limit << endl; + + return limit; +} + +// -------------------------------------------------------------------------- +// +// Returns the mapped value from the Matrix +// +Double_t MMySuperCutsCalc::GetVal(Int_t i) const +{ + + Double_t val = (*fMatrix)[fMap[i]]; + + + //*fLog << "MMySuperCutsCalc::GetVal; i, fMatrix, fMap, val = " + // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl; + + + return val; +} + +// -------------------------------------------------------------------------- +// +// You can use this function if you want to use a MHMatrix instead of the +// given containers. This function adds all necessary columns to the +// given matrix. Afterward you should fill the matrix with the corresponding +// data (eg from a file by using MHMatrix::Fill). If you now loop +// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process +// will take the values from the matrix instead of the containers. +// +void MMySuperCutsCalc::InitMapping(MHMatrix *mat) +{ + if (fMatrix) + return; + + fMatrix = mat; + + // //fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta"); +// fMap[0] = fMatrix->AddColumn("MHillas.fWidth"); +// fMap[1] = fMatrix->AddColumn("MHillas.fLength"); +// fMap[2] = fMatrix->AddColumn("MHillas.fSize"); +// fMap[3] = fMatrix->AddColumn("MHillas.fMeanX"); +// fMap[4] = fMatrix->AddColumn("MHillas.fMeanY"); +// fMap[5] = fMatrix->AddColumn("MHillasSrc.fDist"); +// fMap[6] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)"); +// //fMap[7] = fMatrix->AddColumn("MPointingPos.fZd"); +// fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)"); +// fMap[8] = fMatrix->AddColumn("MNewImagePar.fConc"); +// fMap[9] = fMatrix->AddColumn("MHillasExt.fAsym"); +// fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1"); + + fMap[0] = fMatrix->AddColumn("MHillas.fWidth"); + fMap[1] = fMatrix->AddColumn("MHillas.fLength"); + fMap[2] = fMatrix->AddColumn("MHillas.fSize"); + fMap[3] = fMatrix->AddColumn("MHillas.fMeanX"); + fMap[4] = fMatrix->AddColumn("MHillas.fMeanY"); + fMap[5] = fMatrix->AddColumn("MHillasSrc.fDist"); + fMap[6] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)"); + fMap[7] = fMatrix->AddColumn("MNewImagePar.fConc"); + fMap[8] = fMatrix->AddColumn("MHillasExt.fAsym"); + fMap[9] = fMatrix->AddColumn("MNewImagePar.fLeakage1"); + + //fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); + //fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc"); + //fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1"); +} + +// --------------------------------------------------------------------------- +// +// Evaluate dynamical supercuts +// +// set hadronness to 0.25 if cuts are fullfilled +// 0.75 otherwise +// +Int_t MMySuperCutsCalc::Process() +{ + + + //const Double_t kNomLogSize = 7.6;//=log(2000) 4.1; + const Double_t kNomLogSize = log(fSizeCutLow); + + const Double_t kNomCosZA = 1.0; + + //const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta(); + + // const Double_t width0 = fMatrix ? GetVal(0) : fHil->GetWidth(); +// const Double_t length0 = fMatrix ? GetVal(1) : fHil->GetLength(); +// const Double_t size = fMatrix ? GetVal(2) : fHil->GetSize(); +// const Double_t meanx = fMatrix ? GetVal(3) : fHil->GetMeanX(); +// const Double_t meany = fMatrix ? GetVal(4) : fHil->GetMeanY(); +// const Double_t dist0 = fMatrix ? GetVal(5) : fHilSrc->GetDist(); +// // const Double_t theta = fMatrix ? GetVal(7) : fPointingPos->GetZd(); +// const Double_t conc = fMatrix ? GetVal(8) : fNewImagePar->GetConc(); +// Double_t asym = fMatrix ? GetVal(9) : fHillasExt->GetAsym(); +// const Double_t leakage = fMatrix ? GetVal(10): fNewImagePar->GetLeakage1(); + +// const Double_t theta =0; + + const Double_t width0 = fMatrix ? GetVal(0) : fHil->GetWidth(); + const Double_t length0 = fMatrix ? GetVal(1) : fHil->GetLength(); + const Double_t size = fMatrix ? GetVal(2) : fHil->GetSize(); + const Double_t meanx = fMatrix ? GetVal(3) : fHil->GetMeanX(); + const Double_t meany = fMatrix ? GetVal(4) : fHil->GetMeanY(); + const Double_t dist0 = fMatrix ? GetVal(5) : fHilSrc->GetDist(); + const Double_t conc = fMatrix ? GetVal(7) : fNewImagePar->GetConc(); + Double_t asym = fMatrix ? GetVal(8) : fHillasExt->GetAsym(); + const Double_t leakage = fMatrix ? GetVal(9): fNewImagePar->GetLeakage1(); + + + Double_t help; + if (!fMatrix) + help = TMath::Sign(fHilExt->GetM3Long(), fHilSrc->GetCosDeltaAlpha()); + + // const Double_t asym0 = fMatrix ? GetVal(8) : help; +// const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc(); +// const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1(); + + const Double_t newdist = dist0 * fMm2Deg; + + const Double_t dist2 = meanx*meanx + meany*meany; + const Double_t dist = sqrt(dist2) * fMm2Deg; + const Double_t dd2 = dist*dist; + + const Double_t dmls = log(size) - kNomLogSize; + Double_t dmls2 = dmls * dmls; + + //const Double_t dmcza = cos(theta/kRad2Deg) - kNomCosZA; + + const Double_t length = length0 * fMm2Deg; + const Double_t width = width0 * fMm2Deg; + //const Double_t asym = asym0 * fMm2Deg; + + // *fLog << "Using SizeCut = " << fSizeCut << endl; + /* + *fLog << "newdist, length, width, asym, dist, conc, leakage = " + << newdist << ", " << length << ", " << width << ", " + << asym << ", " << dist << ", " << conc << ", " << leakage + << endl; + + *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = " + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + fHadronness = hadronness; + + << CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLeakage1Up(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLeakage1Lo(), dmls, dmcza, dmls2, dd2) << ", " + << endl; + */ + + double dmcza = 0; + + // cout << fSizeCutLow << " " << fSizeCutUp << endl; + + if(fNoDistCut) + { + if ( size > fSizeCutLow && size < fSizeCutUp && + + //newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmls2) && + //newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmls2) && + + length< CtsMCut (fSuper->GetLengthUp(), dmls, dmls2, dmcza, dd2) && + length> CtsMCut (fSuper->GetLengthLo(), dmls, dmls2, dmcza, dd2) && + + width< CtsMCut (fSuper->GetWidthUp(), dmls, dmls2, dmcza, dd2) && + width > CtsMCut (fSuper->GetWidthLo(), dmls, dmls2, dmcza, dd2) && + + conc< CtsMCut (fSuper->GetConcUp(), dmls, dmls2, dmcza, dd2) && + conc > CtsMCut (fSuper->GetConcLo(), dmls, dmls2, dmcza, dd2) && + + asym< CtsMCut (fSuper->GetAsymUp(), dmls, dmls2, dmcza, dd2) && + asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmls2, dmcza, dd2) && + + leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls,dmls2, dmcza, dd2) && + leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmls2, dmcza, dd2) + + + ) + + fHadronness->SetHadronness(0.25); + else + fHadronness->SetHadronness(0.75); + + } + else + { + //cout << asym << endl; + if (size > fSizeCutLow && size < fSizeCutUp && + + newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmls2, dmcza, dd2) && + newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmls2, dmcza, dd2) && + + length < CtsMCut (fSuper->GetLengthUp(),dmls, dmls2, dmcza, dd2) && + length > CtsMCut (fSuper->GetLengthLo(),dmls, dmls2, dmcza, dd2) && + + width < CtsMCut (fSuper->GetWidthUp(),dmls, dmls2, dmcza, dd2) && + width > CtsMCut (fSuper->GetWidthLo(),dmls, dmls2, dmcza, dd2) && + + asym< CtsMCut (fSuper->GetAsymUp(), dmls, dmls2, dmcza, dd2) && + asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmls2, dmcza, dd2) && + + conc< CtsMCut (fSuper->GetConcUp(), dmls, dmls2, dmcza, dd2) && + conc > CtsMCut (fSuper->GetConcLo(), dmls, dmls2, dmcza, dd2) && + + leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls,dmls2, dmcza, dd2) && + leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmls2, dmcza, dd2) + + ) + + fHadronness->SetHadronness(0.25); + else + fHadronness->SetHadronness(0.75); + } + + + + + + + + //*fLog << "SChadroness = " << fHadronness->GetHadronness() << endl; + + fHadronness->SetReadyToSave(); + + return kTRUE; +} + + + + + +// --------------------------------------------------------------------------- +// +// Return kTRUE if the event pass the dist cut, ie., if it is a gamma-like +// event. +// +Bool_t MMySuperCutsCalc::CalcDistCut(MHillasSrc* hillasSrc) +{ + + + + //const Double_t kNomLogSize = 7.6;//=log(2000) 4.1; + const Double_t kNomLogSize = log(fSizeCutLow); + + const Double_t kNomCosZA = 1.0; + + //const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta(); + + // const Double_t width0 = fMatrix ? GetVal(0) : fHil->GetWidth(); +// const Double_t length0 = fMatrix ? GetVal(1) : fHil->GetLength(); +// const Double_t size = fMatrix ? GetVal(2) : fHil->GetSize(); +// const Double_t meanx = fMatrix ? GetVal(3) : fHil->GetMeanX(); +// const Double_t meany = fMatrix ? GetVal(4) : fHil->GetMeanY(); +// const Double_t dist0 = fMatrix ? GetVal(5) : fHilSrc->GetDist(); +// // const Double_t theta = fMatrix ? GetVal(7) : fPointingPos->GetZd(); +// const Double_t conc = fMatrix ? GetVal(8) : fNewImagePar->GetConc(); +// Double_t asym = fMatrix ? GetVal(9) : fHillasExt->GetAsym(); +// const Double_t leakage = fMatrix ? GetVal(10): fNewImagePar->GetLeakage1(); + +// const Double_t theta =0; + + const Double_t width0 = fMatrix ? GetVal(0) : fHil->GetWidth(); + const Double_t length0 = fMatrix ? GetVal(1) : fHil->GetLength(); + const Double_t size = fMatrix ? GetVal(2) : fHil->GetSize(); + const Double_t meanx = fMatrix ? GetVal(3) : fHil->GetMeanX(); + const Double_t meany = fMatrix ? GetVal(4) : fHil->GetMeanY(); + const Double_t dist0 = fMatrix ? GetVal(5) : fHilSrc->GetDist(); + const Double_t conc = fMatrix ? GetVal(7) : fNewImagePar->GetConc(); + Double_t asym = fMatrix ? GetVal(8) : fHillasExt->GetAsym(); + const Double_t leakage = fMatrix ? GetVal(9): fNewImagePar->GetLeakage1(); + + + Double_t help; + if (!fMatrix) + help = TMath::Sign(fHilExt->GetM3Long(), fHilSrc->GetCosDeltaAlpha()); + + // const Double_t asym0 = fMatrix ? GetVal(8) : help; +// const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc(); +// const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1(); + + const Double_t newdist = dist0 * fMm2Deg; + + const Double_t dist2 = meanx*meanx + meany*meany; + const Double_t dist = sqrt(dist2) * fMm2Deg; + const Double_t dd2 = dist*dist; + + const Double_t dmls = log(size) - kNomLogSize; + Double_t dmls2 = dmls * dmls; + + //const Double_t dmcza = cos(theta/kRad2Deg) - kNomCosZA; + + const Double_t length = length0 * fMm2Deg; + const Double_t width = width0 * fMm2Deg; + //const Double_t asym = asym0 * fMm2Deg; + + + double dmcza = 0; + + if ( + newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmls2, dmcza, dd2) && + newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmls2, dmcza, dd2) + ) + return kTRUE; + + + return kFALSE; +} + + +// // --------------------------------------------------------------------------- +// // +// // +// // +// Double_t MMySuperCutsCalc::Calc(MMySuperCuts* super, MHillas* hillas, MHillasSrc* hillasSrc, MHadronness* hadronness) +// { + +// MGeomCamMagic cam; +// fMm2Deg = cam.GetConvMm2Deg(); + +// fSuper = super; // Container with the optimazed cuts +// fHil = hillas; +// fHilSrc = hillasSrc; +// fHadronness = hadronness; + +// Process(); + +// return fHadronness->GetHadronness(); +// } + Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCutsCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCutsCalc.h (revision 6310) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MMySuperCutsCalc.h (revision 6310) @@ -0,0 +1,114 @@ +#ifndef MARS_MMySuperCutsCalc +#define MARS_MMySuperCutsCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MParList; +class MHillas; +class MHillasSrc; +class MHillasExt; +class MNewImagePar; +class MMcEvt; +class MCerPhotEvt; +class MGeomCam; +class MHadronness; +class MHMatrix; +class MMySuperCuts; +class MPointingPos; +class MNewImagePar; +class MHillasExt; + +class MMySuperCutsCalc : public MTask +{ + private: + + MHillas *fHil; + MHillasSrc *fHilSrc; + MHillasExt *fHilExt; + MNewImagePar *fNewPar; + MMcEvt *fMcEvt; + MHadronness *fHadronness; //! output container for hadronness + MMySuperCuts *fSuper; // container for supercut parameters + MPointingPos *fPointingPos; + MNewImagePar *fNewImagePar; + MHillasExt *fHillasExt; + + TString fHadronnessName; // name of container to store hadronness + TString fHilName; + TString fHilSrcName; + TString fHilExtName; + TString fNewParName; + TString fSuperName; // name of container for supercut parameters + + Double_t fMm2Deg; //! + + Int_t fMap[11]; //! + MHMatrix *fMatrix; //! + + Bool_t fNoDistCut; + Float_t fSizeCutLow,fSizeCutUp ; + + + Int_t PreProcess(MParList *pList); + Int_t Process(); + + Double_t GetVal(Int_t i) const; + + Double_t CtsMCut(const Double_t* a, Double_t ls, Double_t ls2, Double_t ct, Double_t dd2) const; + + + public: + + MMySuperCutsCalc(const char *hilname="MHillas", + const char *hilsrcname="MHillasSrc", + const char *name=NULL, const char *title=NULL); + + void SetHadronnessName(const TString name) { fHadronnessName = name; } + TString GetHadronnessName() const { return fHadronnessName; } + + void InitMapping(MHMatrix *mat); + void StopMapping() { InitMapping(NULL); } + + + /* Double_t Calc(MMySuperCuts* super, MHillas* hillas, MHillasSrc* hillasSrc, MHadronness* hadronness); */ + + Bool_t CalcDistCut(MHillasSrc* hillasSrc); + + + void SetNoDistCut(Bool_t flag) { fNoDistCut = flag; } + void SetSizeCuts(Float_t sizeLow, Float_t sizeUp) + { + fSizeCutLow = sizeLow; + fSizeCutUp = sizeUp; + } + + + ClassDef(MMySuperCutsCalc, 0) // A class to evaluate the Supercuts +}; + +#endif + + + + + + + + + + + + + + + + + + +