Index: /trunk/MagicSoft/Mars/mtemp/mmpi/MApplySupercuts.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mmpi/MApplySupercuts.cc (revision 6681) +++ /trunk/MagicSoft/Mars/mtemp/mmpi/MApplySupercuts.cc (revision 6681) @@ -0,0 +1,262 @@ +/* +======================================================================== *\ +! +! * +! * 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): Nicola Galante 02/2005 +! +! Copyright: MAGIC Software Development, 2004 +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MApplySupercuts +// +// Apply supercuts to a file containing Hillas parameters +// (usually a star file). +// +// Input: +// rootfile (starfile) of Hillas parameters +// +// Output: +// MHadronness estimated on event +// +///////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include + +#include "MGeomCamMagic.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MRawEvtHeader.h" +#include "MRawRunHeader.h" +#include "MSupercuts.h" +#include "MHadronness.h" +#include "MHillas.h" +#include "MHillasSrc.h" +#include "MHFindSignificance.h" +#include "MEnergyEst.h" +#include "MApplySupercuts.h" + +ClassImp(MApplySupercuts); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor +// +MApplySupercuts::MApplySupercuts(const char *name, const char *title) + : fRunHeader(0), fEvtHeader(0), fSupercuts(0) +{ + fName = name ? name : "MApplySupercuts"; + fTitle = title ? title : "Task to apply supercuts to star files"; + + fPhe = kTRUE; + fPlot = kFALSE; +} + + +// -------------------------------------------------------------------------- +// +Int_t MApplySupercuts::PreProcess(MParList *pList) +{ + fGeom = (MGeomCamMagic *)pList->FindCreateObj("MGeomCamMagic"); + if (!fGeom) + { + *fLog << err << "MGeomCamMagic not found... abort." << endl; + return kFALSE; + } + + fRunHeader = (MRawRunHeader*)pList->FindCreateObj("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << err << "MRawRunHeader not found... abort." << endl; + return kFALSE; + } + + fEvtHeader = (MRawEvtHeader *)pList->FindObject("MRawEvtHeader"); + if (!fEvtHeader) + { + *fLog << err << "MRawEvtHeader not found... but ok." << endl; + } + + fHillas = (MHillas *)pList->FindObject("MHillas"); + if (!fHillas) + { + *fLog << err << "MHillas not found... abort." << endl; + return kFALSE; + } + + fHillasSrc = (MHillasSrc *)pList->FindObject("MHillasSrc"); + if (!fHillasSrc) + { + *fLog << err << "MHillasSrc not found... abort." << endl; + return kFALSE; + } + + fHadronness = (MHadronness *)pList->FindCreateObj("MHadronness",AddSerialNumber("MHadronness")); + if (!fHadronness) + { + *fLog << err << "MHadronness not found... abort." << endl; + return kFALSE; + } + + fEnergyEst = (MEnergyEst *)pList->FindCreateObj("MEnergyEst"); + if (!fEnergyEst) + { + *fLog << err << "MEnergyEst not found... abort." << endl; + return kFALSE; + } + + if(fPlot){ + fAlpha = new TH1F("fHistAlpha","Alpha Plot",30,0,90); + if (!fAlpha) + { + *fLog << err << "Impossible to create Alpha histogram... abort." << endl; + return kFALSE; + } + } + + fHFindSigma = (MHFindSignificance *)pList->FindCreateObj("MHFindSignificance"); + if (!fHFindSigma) + { + *fLog << err << "MHFindSignificance not found... abort." << endl; + return kFALSE; + } + + // Start reading supercuts + fSupercuts = (MSupercuts *)pList->FindCreateObj("MSupercuts",AddSerialNumber("MSupercuts")); + if (!fSupercuts) + { + *fLog << err << "MSupercuts not found... abort." << endl; + return kFALSE; + } + // read the cuts coefficients + TFile inparam(fSCFilename); + fSupercuts->Read("MSupercuts"); + inparam.Close(); + *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from file " + << fSCFilename << endl; + + TArrayD params = fSupercuts->GetParameters(); + + for (Int_t i=0; i<8; i++) + { + fLengthUp[i] = params[i]; + fLengthLo[i] = params[i+8]; + fWidthUp[i] = params[i+16]; + fWidthLo[i] = params[i+24]; + fDistUp[i] = params[i+32]; + fDistLo[i] = params[i+40]; + } + + // end reading supercuts + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +Int_t MApplySupercuts::Process() +{ + fHadronness->SetHadronness(kHadronHadronness); + + Float_t fMm2Deg = fGeom->GetConvMm2Deg(); + Double_t log3000 = TMath::Log(fSizeOffset); + Float_t fsize = fHillas->GetSize(); + if(fPhe) + fsize /= kPhe2Ph; + Float_t fdist = fHillasSrc->GetDist()*fMm2Deg; + Double_t logsize = TMath::Log((Double_t)fsize); + Double_t lgsize = logsize-log3000; + Double_t lgsize2 = lgsize*lgsize; + Double_t dist2 = (Double_t)fdist*fdist - fDistOffset*fDistOffset; + //Double_t dist2 = (fDist-fDistOffset)*(fDist-fDistOffset); + + // parameters: + Float_t flength = (fHillas->GetLength()) * fMm2Deg; + Float_t fwidth = (fHillas->GetWidth())*fMm2Deg; + //Float_t fsize = fHillas.GetSize()/0.18; + //Float_t fmeanx = (fHillas.GetMeanX())*fMm2Deg; + //Float_t fmeany = (fHillas.GetMeanY())*fMm2Deg; + //Float_t falpha = fHillasSrc.GetAlpha(); + //Float_t fDCAdelta = fHillasSrc.GetDCADelta(); + + + + if ( flength < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) && + flength > CalcLimit(fLengthLo, lgsize, lgsize2, dist2) && + fwidth < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) && + fwidth > CalcLimit(fWidthLo, lgsize, lgsize2, dist2) && + fdist < CalcLimit(fDistUp, lgsize, lgsize2, dist2) && + fdist > CalcLimit(fDistLo, lgsize, lgsize2, dist2) ) + { + // gamma candidates! + fHadronness->SetHadronness(kGammaHadronness); + if( fPlot && fHillas->GetSize()>fSizeLow && fHillas->GetSize()Fill(TMath::Floor(TMath::Abs(fHillasSrc->GetAlpha()))); + } + + return kTRUE; +} + + +Int_t MApplySupercuts::PostProcess() +{ + if(fPlot){ + TFile f("shit_file.root","RECREATE"); + fAlpha->Write(); + f.Close(); + TCanvas c; + if(fHFindSigma!=NULL){ + cout << "W el cogno " << fAlphaSignal <FindSigma((TH1 *)fAlpha,fAlphaMin,fAlphaMax,fDegree,fAlphaSignal,kTRUE,kTRUE,kTRUE); + *fLog << "SIGNIFICANCE = " << fHFindSigma->GetSignificance() << endl; + c.SaveAs("Alpha.gif"); + } + else + *fLog << err << "ERROR: fHFindSigma already deleted" << endl; + } + return kTRUE; +} + + +Double_t MApplySupercuts::CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2) +{ + + Double_t limit = a[0] + a[1] * dd2 + + ls * (a[3] + a[4] * dd2) + + ls2 * (a[6] + a[7] * dd2); + + return limit; +} + + +void MApplySupercuts::SetPrintOutSignificance(Double_t bgalphamin=50.0, Double_t bgalphamax=90.0, + Double_t alphasig=15.0, Int_t degree=4) +{ + fAlphaMin = bgalphamin; + fAlphaMax = bgalphamax; + fDegree = degree; + fAlphaSignal = alphasig; + + fPlot = kTRUE; +} + Index: /trunk/MagicSoft/Mars/mtemp/mmpi/MApplySupercuts.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mmpi/MApplySupercuts.h (revision 6681) +++ /trunk/MagicSoft/Mars/mtemp/mmpi/MApplySupercuts.h (revision 6681) @@ -0,0 +1,91 @@ +#ifndef MARS_MApplySupercuts +#define MARS_MApplySupercuts + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MHadronness +#include "MHadronness.h" +#endif + +class MGeomCamMagic; +class MParList; +class MRawEvtHeader; +class MRawRunHeader; +class MSupercuts; +class MHillas; +class MHillasSrc; +class MHadronness; +class MHFindSignificance; +class MEnergyEst; +class TH1F; + +class MApplySupercuts : public MTask +{ + private: + + // Constant to convert size in photoelectrons to photons + static const Float_t kPhe2Ph = 0.18; + static const Double_t kHadronHadronness = 0.75; + static const Double_t kGammaHadronness = 0.25; + static const Double_t fDistOffset = 0.9; + static const Double_t fSizeOffset = 3000.; + + Bool_t fPhe; + Bool_t fPlot; + Double_t fLengthUp[8]; + Double_t fLengthLo[8]; + Double_t fWidthUp[8]; + Double_t fWidthLo[8]; + Double_t fDistUp[8]; + Double_t fDistLo[8]; + Double_t fAlphaMin; // Minimum alpha of background region (deg) + Double_t fAlphaMax; // Maximum alpha of background region (deg) + Double_t fAlphaSignal; // Maximum alpha of signal region (deg - Min=0) + Int_t fDegree; // Degree of polinomial interpolating function + Double_t fSizeUp; + Double_t fSizeLow; + + TString fSCFilename; + + MGeomCamMagic *fGeom; + MRawRunHeader *fRunHeader; + MRawEvtHeader *fEvtHeader; + MSupercuts *fSupercuts; + MHillas *fHillas; + MHillasSrc *fHillasSrc; + + TH1F *fAlpha; + MHadronness *fHadronness; + MHFindSignificance *fHFindSigma; + MEnergyEst *fEnergyEst; + + Double_t CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + public: + MApplySupercuts(const char *name=0, const char *title=0); + + void SetPhotoelectrons() { fPhe = kTRUE; } + void SetPhotons() { fPhe = kFALSE; } + void SetSCFilename(TString filename) { fSCFilename = filename; } + void SetPrintOutSignificance(Double_t bgalphamin, Double_t bgalphamax, + Double_t alphasig, Int_t kDegree); + void SetPlotON() { fPlot = kTRUE; }; + void SetPlotOFF() { fPlot = kFALSE; }; + void SetBGAlphaMin(Double_t alphamin) { fAlphaMin=alphamin; SetPlotON(); }; + void SetBGAlphaMax(Double_t alphamax) { fAlphaMax=alphamax; SetPlotON(); }; + void SetAlphaSig(Double_t alphasig) { fAlphaSignal=alphasig; SetPlotON(); }; + void SetPolDegree(Int_t deg) { fDegree=deg; SetPlotON(); }; + void SetPlotSizeLow(Double_t size) { fSizeLow = size; }; + void SetPlotSizeUp(Double_t size) { fSizeUp = size; }; + + ClassDef(MApplySupercuts, 1) // Task to apply supercuts +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/mmpi/hit.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mmpi/hit.cc (revision 6681) +++ /trunk/MagicSoft/Mars/mtemp/mmpi/hit.cc (revision 6681) @@ -0,0 +1,800 @@ +#include +#include +#include +#include +#include +#include + +#include "MHMatrix.h" +#include "MParList.h" +#include "MTaskList.h" +#include "MWriteRootFile.h" +#include "MReadMarsFile.h" +#include "MReadReports.h" +#include "MApplySupercuts.h" +#include "MHFindSignificance.h" +#include "MEvtLoop.h" +#include "MProgressBar.h" +#include "MContinue.h" +#include "MGeomApply.h" +#include "MHadronness.h" +#include "MFilterList.h" +#include "MF.h" +#include "MFindSupercutsONOFFThetaLoop.h" +#include "MRanForestFill.h" + +#include "MArgs.h" +#include "MArray.h" +#include "MDirIter.h" + +#include "MStatusDisplay.h" + +#include "MSequence.h" +#include "MJStar.h" +#include "MH3.h" +#include "MRFEnergyEst.h" + +#include "MRanForestCalc.h" +#include "MTaskInteractive.h" +#include "MHillas.h" +//#include "MFHadAlpha.h" + +#include "MLog.h" +#include "MLogManip.h" +#include "MPrint.h" + +using namespace std; + + +Int_t PreProcess(MParList *plist); +Int_t Process(); + +static void StartUpMessage() +{ + gLog << all << endl; + + // 1 2 3 4 5 6 + // 123456789012345678901234567890123456789012345678901234567890 + gLog << "========================================================" << endl; + gLog << " Melibea - MARS V" << MARSVER << endl; + gLog << " MARS -- Supercuts Hillas Image Treatment" << endl; + gLog << " Compiled on <" << __DATE__ << ">" << endl; + gLog << " Using ROOT v" << ROOTVER << endl; + gLog << "========================================================" << endl; + + gLog << endl; +} + +static void Usage() +{ + // 1 2 3 4 5 6 7 8 + // 12345678901234567890123456789012345678901234567890123456789012345678901234567890 + gLog << all << endl; + gLog << "Sorry the usage is:" << endl; + + +} + +int main(int argc, char **argv) +{ + StartUpMessage(); + + // + // Evaluate arguments + // + MArgs arg(argc, argv, kTRUE); + + if (arg.HasOnly("-V") || arg.HasOnly("--version")) + return 0; + + if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) + { + Usage(); + return -1; + } + + gLog.Setup(arg); + + const Bool_t kOptimize = arg.HasOnlyAndRemove("--optimize"); + const Bool_t kSupercuts = arg.HasOnlyAndRemove("--supercuts"); + const Bool_t kOverwrite = arg.HasOnlyAndRemove("-f"); + const Bool_t kPrintSeq = arg.HasOnlyAndRemove("--print-seq"); + const Bool_t kBatch = arg.HasOnlyAndRemove("-b"); + const Bool_t kPlot = arg.HasOnlyAndRemove("-p"); + const Bool_t kMc = arg.HasOnlyAndRemove("--mc"); + const Bool_t kUseSeq = !arg.HasOnlyAndRemove("--nouseseq"); + + // Variable that decides whether the optimized cuts are used + // in the test sample. + const Bool_t TestParams = !arg.HasOnlyAndRemove("--opttest"); + const Bool_t CombineCosThetaBinsForTrainSample = arg.HasOnlyAndRemove("--combcosthtrain"); + const Bool_t CombineCosThetaBinsForTestSample = arg.HasOnlyAndRemove("--combcosthtest"); + const Double_t whichfractiontrain = arg.GetFloatAndRemove("--ontrainfrac=",0.5); + const Double_t whichfractiontest = arg.GetFloatAndRemove("--ontestfrac=",0.5); + const Double_t whichfractiontrainOFF = arg.GetFloatAndRemove("--offtrainfrac=",0.5); + const Double_t whichfractiontestOFF = arg.GetFloatAndRemove("--offtestfrac=",0.5); + const Double_t gammaeff = arg.GetFloatAndRemove("--geff=",0.6); + const Double_t alphasig = arg.GetFloatAndRemove("--alphasig=",6); + const Double_t alphabkgmin = arg.GetFloatAndRemove("--alphabgmin=",20); + const Double_t alphabkgmax = arg.GetFloatAndRemove("--alphabgmax=",80); + const Double_t SizeLow = arg.GetFloatAndRemove("--sizelow=",60); + const Double_t SizeUp = arg.GetFloatAndRemove("--sizeup=",1000000); + const Double_t LeakageMax = arg.GetFloatAndRemove("--leakagemax=",0.25); + const Double_t DistMax = arg.GetFloatAndRemove("--distmax=",1.0); + const Double_t DistMin = arg.GetFloatAndRemove("--distmin=",0.2); + const Double_t fThetaMax = arg.GetFloatAndRemove("--zdmax=",30.); + const Double_t fThetaMin = arg.GetFloatAndRemove("--zdmin=",0.); + const Int_t NAlphaBins = arg.GetIntAndRemove("--nalphabin=",30); + const Double_t AlphaBinLow = arg.GetFloatAndRemove("--alphabinlow=",0); + const Double_t AlphaBinUp = arg.GetFloatAndRemove("--alphabinup=",90); + const Bool_t NormFactorFromAlphaBkg = !arg.HasOnlyAndRemove("--normbgmeth"); + const Bool_t UseFittedQuantities = !arg.HasOnlyAndRemove("--usefittedq"); + const Bool_t NotUseTheta = !arg.HasOnlyAndRemove("--notuseth"); + const Bool_t UseStaticCuts = !arg.HasOnlyAndRemove("--staticcuts"); + const Bool_t ReadInitParamsFromAsciiFile = !arg.HasOnlyAndRemove("--ascii"); + const Bool_t kPrintFiles = !arg.HasOnlyAndRemove("--printseq"); + const Int_t NInitSCPar = arg.GetIntAndRemove("--ninitscpar=",104); + const Int_t degree = arg.GetIntAndRemove("--degree=",2); + const TString InitSCParamAsciiFile = arg.GetStringAndRemove("--inscparascii=","mtemp/mmpi/asciifiles/StartingValuesForSmallSizes1.txt"); + const TString ONDataFilename = arg.GetStringAndRemove("--inondataopt=","./*CrabNebula*.root"); + const TString OFFDataFilename = arg.GetStringAndRemove("--inoffdataopt=","./*OffCrab*.root"); + const TString PathForFiles = arg.GetStringAndRemove("--outopt=","./"); + + const TString kInpath = arg.GetStringAndRemove("--ind=", ""); + const TString kOutpath = arg.GetStringAndRemove("--out=", ""); + const TString kSCfile = arg.GetStringAndRemove("--scf=", ""); + const TString kSource = arg.GetStringAndRemove("--src=", ""); + const TString kDate = arg.GetStringAndRemove("--date=", ""); + const TString kRFTreeFile = arg.GetStringAndRemove("--rftree=", ""); + const Double_t kAlphaMin = arg.GetFloatAndRemove("--alphamin=",-1); + const Double_t kAlphaMax = arg.GetFloatAndRemove("--alphamax=",-1); + const Double_t kAlphaSig = arg.GetFloatAndRemove("--alphasignal=",-1); + const Int_t kDegree = arg.GetIntAndRemove("--poldeg=",-1); + const Int_t kNumTrees = arg.GetIntAndRemove("--ntrees=",100); + + + if (arg.GetNumOptions()>0) + { + gLog << warn << "WARNING - Unknown commandline options..." << endl; + arg.Print("options"); + gLog << endl; + return -1; + } + + if (arg.GetNumArguments()!=1) + { + Usage(); + return -1; + } + // + // Setup sequence file and check for its existance + // + const TString kSequence = arg.GetArgumentStr(0); + + if (gSystem->AccessPathName(kSequence, kFileExists)) + { + gLog << err << "Sorry, sequence file '" << kSequence << "' doesn't exist." << endl; + return -1; + } + + // + // Setup sequence and check its validity + // + MSequence seq(kSequence); + if (kPrintSeq) + { + gLog << all; + gLog.Separator(kSequence); + seq.Print(); + gLog << endl; + } + if (!seq.IsValid()) + { + gLog << err << "Sequence read but not valid!" << endl << endl; + return -1; + } + + // + // Process print options + // + + MDirIter iter; + + if(kUseSeq) { + gLog << inf; + gLog << "Calculate image parameters from sequence "; + gLog << seq.GetName() << endl; + gLog << endl; + +/* + const Int_t n0 = seq.SetupDatRuns(iter, kInpath, "I"); + const Int_t n1 = seq.GetNumDatRuns(); + + if (n0==0) + { + gLog << err << "ERROR - No input files of sequence found!" << endl; + return kFALSE; + } + if (n0!=n1) + { + gLog << err << "ERROR - Number of files found (" << n0 << ") doesn't match number of files in sequence (" << n1 << ")" << endl; + return kFALSE; + } +*/ + } + + + + // + // Initialize root + // + MArray::Class()->IgnoreTObjectStreamer(); + MParContainer::Class()->IgnoreTObjectStreamer(); + + TApplication app("Star", &argc, argv); + if (!gROOT->IsBatch() && !gClient || gROOT->IsBatch() && !kBatch) + { + gLog << err << "Bombing... maybe your DISPLAY variable is not set correctly!" << endl; + return 1; + } + + + // **************************************************** // + // Setup flags for optimization // + // **************************************************** // + + + // Name of the root file where alpha distributions, TTree objects + // with info about the events and cuts applied and info support histograms + // will be stored. + // Write only the name of the file. The Path + // is the one defined previously + + TString RootFilename = kOutpath + "RootFileDynCuts.root"; + TString ONFiles = ONDataFilename + "*_I_*.root"; + TString OFFFiles = OFFDataFilename + "*_I_*.root"; + // Vector containing the theta bins in which data (ON/OFF train/test) + // will be divided. Actually this vector contains the cosinus of + // these theta bins. The dimension of the vector is N+1, where + // N is the number of theta bins intended. The first component of the + // vector is the low bin edge of the first bin, and the last + // vector component the upper bin edge of the last bin. + + TArrayD CosThetaRangeVector(2); + CosThetaRangeVector[1] = 0.866; //30 + CosThetaRangeVector[0] = 0.5; // 60 + + // Object of MCT1FindSupercutsONOFFThetaLoop created, data that was specified + // above is introduced and ... and the party starts. + + MFindSupercutsONOFFThetaLoop FindSupercuts("MFindSupercutsONOFFThetaLoop", + "Optimizer for the supercuts"); + + + // ****************************************************** // + // Calculate and/or apply dynamical cuts // + // ****************************************************** // + + // Hillas data file + //TString datafile = kInpath + "*_I_CrabNebula_E.root"; + + TString datafile; + if (kUseSeq) + datafile = kInpath + kDate + "*_I_" + kSource + "*.root"; + else + datafile = kInpath + "*" + kSource + "*.root"; + + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + if (kUseSeq) + read.AddFiles(iter); + else + read.AddFile(datafile); + + cout << " datafile is : " << datafile << endl; + + MReadReports readreal; + readreal.AddTree("Events", "MTime.", kTRUE); + readreal.AddTree("Drive"); + readreal.AddTree("EffectiveOnTime"); + if (kUseSeq) + readreal.AddFiles(iter); + else + readreal.AddFile(datafile); + +// readreal.AddToBranchList("MReportDrive.*"); +// readreal.AddToBranchList("MRawEvtHeader.*"); +// readreal.AddToBranchList("MEffectiveOnTime.*"); + + + MGeomApply geomapl; + + MApplySupercuts appcuts; + appcuts.SetSCFilename(kSCfile); + + MFilterList flist; + MF hadrf("MHadronness.fHadronness<0.5"); + MF sizef("MHillas.fSize>300"); + MF leakmax("MNewImagePar.fLeakage1<0.25"); + MF distmax("MHillasSrc.fDist<300"); + MF distmin("MHillasSrc.fDist>60"); + MF widthmin("MHillas.fWidth>15"); + MF widthmax("MHillas.fWidth<120"); + + TString ThetaCutMinString ("MPointingPos.fZd>"); // new! // for magic + ThetaCutMinString += fThetaMin; + MContinue ThetaCutMin(ThetaCutMinString); + ThetaCutMin.SetInverted(); + + TString ThetaCutMaxString ("MPointingPos.fZd<"); // new! // for magic + ThetaCutMaxString += fThetaMax; + MContinue ThetaCutMax(ThetaCutMaxString); + ThetaCutMax.SetInverted(); + + MContinue cont("(MHillasSrc.fDist<300) && (MHillasSrc.fDist>60) && (MNewImagePar.fLeakage1<0.25) && (MHillas.fSize>60.)"); + cont.SetInverted(); + + + MRFEnergyEst rfs; + //rfs.SetFile("/datm3/data/RFEnergyEst/EForestsBeforeCutsLZA.root"); + rfs.SetFile("/datm3/data/RFEnergyEst/EForests19990101_10003_I_MCGamTrainHZA_E_10_5.root"); + + + //TString outname = Form("%s{s/_I_/_Q_}", kOutpath.Data()); + MWriteRootFile &write=*new MWriteRootFile(2, Form("%s{s/_I_/_Q_}", kOutpath.Data()), "RECREATE"); + // Write the Events tree + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MImagePar", "Events"); + write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MRawEvtHeader", "Events", kFALSE); + write.AddContainer("MPointingPos", "Events"); + write.AddContainer("MHadronness","Events", kFALSE); + write.AddContainer("MEnergyEst","Events", kFALSE); + + if(kMc) + { + // Monte Carlo + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MMcTrig", "Events"); + // Monte Carlo Headers + write.AddContainer("MMcTrigHeader", "RunHeaders"); + write.AddContainer("MMcConfigRunHeader", "RunHeaders"); + write.AddContainer("MMcCorsikaRunHeader", "RunHeaders"); + write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); + } + else + { + write.AddContainer("MTime", "Events"); + // Run Header + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MBadPixelsCam", "RunHeaders", kFALSE); + write.AddContainer("MGeomCam", "RunHeaders", kFALSE); + //write.AddContainer("MObservatory", "RunHeaders"); + // Drive + //write.AddContainer("MSrcPosCam", "Drive"); + write.AddContainer("MReportDrive", "Drive", kFALSE); + write.AddContainer("MTimeDrive", "Drive", kFALSE); + // Effective On Time + write.AddContainer("MEffectiveOnTime", "EffectiveOnTime", kFALSE); + write.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime", kFALSE); + } + + +/* + write.AddContainer("MHillas","Events"); + write.AddContainer("MHillasExt","Events"); + write.AddContainer("MHillasSrc","Events"); + write.AddContainer("MImagePar","Events"); + write.AddContainer("MNewImagePar","Events"); + write.AddContainer("MPointingPos","Events"); + write.AddContainer("MRawEvtHeader","Events"); + write.AddContainer("MTime","Events", kFALSE); + write.AddContainer("MHadronness","Events"); + write.AddContainer("MEnergyEst","Events"); + // Write the RunHeader tree + write.AddContainer("MRawRunHeader","RunHeaders", kFALSE); + write.AddContainer("MBadPixelsCam","RunHeaders",kFALSE); + write.AddContainer("MGeomCam","RunHeaders",kFALSE); + // Write the Drive tree + write.AddContainer("MReportDrive","Drive",kFALSE); + write.AddContainer("MTimeDrive","Drive",kFALSE); + // Write the EffectiveOnTime tree + write.AddContainer("MEffectiveOnTime","EffectiveOnTime",kFALSE); + write.AddContainer("MTimeEffectiveOnTime","EffectiveOnTime",kFALSE); + + + write.AddContainer("MMcConfigRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcFadcHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcTrigHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcEvt", "Events", kFALSE); + write.AddContainer("MMcTrig", "Events", kFALSE); +*/ + cout << datafile << endl; + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + + MProgressBar bar; + + switch(kSupercuts) + { + case kTRUE: + switch(kOptimize){ + case kTRUE: + //cout << ONFiles << endl; + FindSupercuts.SetPathForFiles(PathForFiles); + FindSupercuts.SetDataONOFFRootFilenames(ONFiles, OFFFiles); + FindSupercuts.SetFractionTrainTestOnOffEvents(whichfractiontrain, + whichfractiontest, + whichfractiontrainOFF, + whichfractiontestOFF); + FindSupercuts.SetGammaEfficiency(gammaeff); + FindSupercuts.SetAlphaSig(alphasig); + // Bkg alpha region is set + FindSupercuts.SetAlphaBkgMin(alphabkgmin); + FindSupercuts.SetAlphaBkgMax(alphabkgmax); + // alpha bkg and signal region set in object FindSupercuts + // are re-checked in order to be sure that make sense + FindSupercuts.CheckAlphaSigBkg(); + // Degree of the polynomials used to fit the ON OFF data is set + FindSupercuts.SetDegreeON(degree); + FindSupercuts.SetDegreeOFF(degree); + // binning for alpha plots is defined + FindSupercuts.SetAlphaPlotBinining(NAlphaBins, AlphaBinLow, AlphaBinUp); + // Size range is defined + FindSupercuts.SetSizeRange(SizeLow, SizeUp); + FindSupercuts.SetFilters(LeakageMax, DistMax, DistMin); + FindSupercuts.SetNormFactorFromAlphaBkg(NormFactorFromAlphaBkg); + FindSupercuts.SetUseFittedQuantities(UseFittedQuantities); + FindSupercuts.SetVariableUseStaticCuts(UseStaticCuts); + FindSupercuts.SetVariableNotUseTheta(NotUseTheta); + FindSupercuts.SetReadMatricesFromFile(kFALSE);// normal case kFALSE + FindSupercuts.SetTrainParameters(kTRUE); + FindSupercuts.SetSkipOptimization(kFALSE); + + FindSupercuts.SetUseInitialSCParams(kTRUE); + FindSupercuts.SetTestParameters(TestParams); + FindSupercuts.SetHadronnessName("MHadSC"); + FindSupercuts.SetHadronnessNameOFF("MHadOFFSC"); + FindSupercuts.SetAlphaDistributionsRootFilename (RootFilename); + FindSupercuts.SetCosThetaRangeVector (CosThetaRangeVector); + +// energy estimation from Thomas Hengstebeck + + + // Names for all root files (matrices, alpha distributions...) are created + FindSupercuts.SetALLNames(); + + if(ReadInitParamsFromAsciiFile) + { + // Initial SC Parameters and steps are retrieved from Ascii file + if(!FindSupercuts.ReadSCParamsFromAsciiFile(InitSCParamAsciiFile,NInitSCPar)) + { + cout << "Initial SC Parameters could not be read from Ascii file " + << InitSCParamAsciiFile << endl + << "Aborting execution of macro... " << endl; + return -1; + } + } + + // Finally loop over all theta bins defined is executed + if (!FindSupercuts.LoopOverThetaRanges()) + { + cout << "Function MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges()" << endl + << "could not be performed" << endl; + } + + // Nex and Significance are computed vs alphasig + if (!FindSupercuts.ComputeNexSignificanceVSAlphaSig()) + { + cout << "Function MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl + << "could not be performed" << endl; + } + + // Several theta bins are combined to produced a single alpha plot (for train and test) + // with single Nex and significances + if (CombineCosThetaBinsForTrainSample || CombineCosThetaBinsForTestSample) + { + if(!FindSupercuts.ComputeOverallSignificance(CombineCosThetaBinsForTrainSample, + CombineCosThetaBinsForTestSample)) + { + cout << "Function MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl + << "could not be performed" << endl; + } + } + break; + case kFALSE: + if(kAlphaMin!=-1) + appcuts.SetBGAlphaMin(kAlphaMin); + if(kAlphaMax!=-1) + appcuts.SetBGAlphaMax(kAlphaMax); + if(kAlphaSig!=-1) + appcuts.SetAlphaSig(kAlphaSig); + if(kDegree!=-1) + appcuts.SetPolDegree(kDegree); + if(kPlot) + appcuts.SetPlotON(); + appcuts.SetPlotSizeLow(SizeLow); + appcuts.SetPlotSizeUp(SizeUp); + cout << "alphasignal=" << kAlphaSig << endl; + //flist.AddToList(&hadrf); + flist.AddToList(&leakmax); + flist.AddToList(&distmax); + flist.AddToList(&distmin); + flist.AddToList(&widthmin); + flist.AddToList(&widthmax); + flist.AddToList(&sizef); + appcuts.SetFilter(&flist); + //write.SetFilter(&flist); + + tlist.AddToList(kMc ? (MTask*)&read : (MTask*)&readreal); + cout << " is montecarlo ? " << kMc << endl; + + tlist.AddToList(&geomapl); + tlist.AddToList(&ThetaCutMin, "Events"); + tlist.AddToList(&ThetaCutMax, "Events"); + tlist.AddToList(&flist, "Events"); + tlist.AddToList(&appcuts, "Events"); + tlist.AddToList(&rfs, "Events"); + tlist.AddToList(&cont, "Events"); + + //tlist.AddToList(&hsigma); +MPrint print1("MTimeDrive"); +MPrint print2("MEffectiveOnTime"); +// tlist.AddToList(&print1, "Drive"); +// tlist.AddToList(&print2, "EffectiveOnTime"); + + tlist.AddToList(&write); + MEvtLoop shitloop; + shitloop.SetParList(&plist); + // + // Start to loop over all events + // + MProgressBar bar; + shitloop.SetProgressBar(&bar); + + if (!shitloop.Eventloop()) + return -1; + + /* + + if (!shitloop.PreProcess()) + return -1; + + while (tlist.Process()) + { + //MHadronness *fH = (MHadronness *)plist.FindObject("MHadronness"); + //cout<< "Hadronness " << fH->GetHadronness() << endl; + } + + if (!shitloop.PostProcess()) + return -1; + */ + tlist.PrintStatistics(); + + + break; + } + break; + case kFALSE: + // PUT THINGS FOR THE RANDOM FOREST + + + gLog << endl; + gLog.Separator("Applying RF to data"); + + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + //MParList plist; + + //MTaskList tlist; + //plist.AddToList(&tlist); + // + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // first event loop: the trees of the random forest are read in + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // + MReadTree readrf("Tree",kRFTreeFile); + readrf.DisableAutoScheme(); + + MRanForestFill rffill; + rffill.SetNumTrees(kNumTrees); + + tlist.AddToList(&readrf); + tlist.AddToList(&rffill); + + // + // Eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetProgressBar(&bar); + + if (!evtloop.Eventloop()) + return 1; + + + // --------------------------------------------------------------- + // + // Apply RF to data + // --------------------------------------------------------------- + // --------------------------------------------------------------- + MTaskList tlist3; + plist.Replace(&tlist3); + + // Calc RF + MRanForestCalc calc; + + // Fill Hist (Only for MC file with Gammas & Hadrons mixed) + //MFillH fillh3a("MHHadronness"); + //MFillH fillh3b("MHRanForest"); + + + // My intercative task to apply cuts + //MTaskInteractive mytask; + //mytask.SetPreProcess(PreProcess); + //mytask.SetProcess(Process); + + // + // Hadronness cut + // +// MFHadAlpha cut; +// cut.SetCutType(MFHadAlpha::kHadCut ); +// cut.SetInverted(); +// MContinue contHad(&cut); + + + + + // Add to TaskList + tlist3.AddToList(kMc ? (MTask*)&read : (MTask*)&readreal); + //tlist3.AddToList(&SizeCut); + tlist3.AddToList(&calc); + tlist3.AddToList(&rfs); + // tlist3.AddToList(&fillh3a); +// tlist3.AddToList(&fillh3b); + + //tlist.AddToList(&mytask); + //tlist.AddToList(&contHad); + + +// tlist3.AddToList(&flist); + tlist3.AddToList(&ThetaCutMax); + tlist3.AddToList(&ThetaCutMin); + tlist3.AddToList(&write); + + + + + // + // Execute your analysis + // + // MProgressBar bar; + evtloop.SetProgressBar(&bar); + if (!evtloop.Eventloop()) + return 1; + + tlist3.PrintStatistics(); + + // plist.FindObject("MHRanForest")->DrawClone(); +// plist.FindObject("MHHadronness")->DrawClone(); + + + + +// END OF PUTTING THINGS FOR THE RANDOM FOREST + break; + } + +//abelardo + //////////////////////////////////////////////////////////////////// + // + // Second event loop: simply copy the tree MMcEvtBasic in the tree + // "OriginalMC" from the input file to the output file: + delete &write; + if(kMc) + //if(0==1) + { + MParList plist2; + MTaskList tlist2; + + plist2.AddToList(&tlist2); + + MReadMarsFile read2("OriginalMC",datafile); +// read2.AddFile(AnalysisFilename); + read2.DisableAutoScheme(); + + tlist2.AddToList(&read2); + + MWriteRootFile writeOrig(2, Form("%s{s/_I_/_Q_}", kOutpath.Data()),"UPDATE"); + writeOrig.AddContainer("MMcEvtBasic", "OriginalMC", kFALSE); + + tlist2.AddToList(&writeOrig); + + MEvtLoop evtloop2; + bar.SetWindowName("Copying OriginalMC tree..."); + evtloop2.SetProgressBar(&bar); + evtloop2.SetParList(&plist2); + + if (!evtloop2.Eventloop()) + return -1; + + tlist2.PrintStatistics(); + } + // end ab + return 0; +} + +Int_t OptimizeCuts() +{ + return 0; +} + + + + + + +// --------------------------------------------------------------------- +// +// Interactive Task +// +// +// Data members of the InteractiveTask +MParList *fParList; +MTaskList *fTaskList; +MHadronness *fHadronness; +MHillas *fHillas; + + +Int_t PreProcess(MParList *plist) +{ + // Get parameter and tasklist, see Process + fParList = plist; + fTaskList = (MTaskList*)plist->FindObject("MTaskList"); + + // Search in the PreProcess for the objects + fHadronness = (MHadronness*)fParList->FindObject("MHadronness"); + fHillas = (MHillas*)fParList->FindObject("MHillas"); + + return kTRUE; +} + +// -------------------------------------------------------------- +// +Int_t Process() +{ + Float_t size = fHillas->GetSize()/0.18; //photons! + Float_t had = fHadronness->GetHadronness(); + + if(size<500) + return kCONTINUE; + + if(size>500 && size <1000 && had>0.32) + return kCONTINUE; + + if(size>1000 && size <2000 && had>0.28) + return kCONTINUE; + + if(size>2000 && size <4000 && had>0.12) + return kCONTINUE; + + if(size>4000 && size <8000 && had>0.21) + return kCONTINUE; + + if(size>8000 && had>0.49) + return kCONTINUE; + + return kTRUE; +} + + +