#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"); rfs.SetFile("/home/pcmagic16/mazin/mars/files/EForestsHZA.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; }