#include "CT1EgyEst.C" void InitBinnings(MParList *plist) { gLog << "InitBinnings" << endl; //-------------------------------------------- MBinning *binse = new MBinning("BinningE"); //binse->SetEdgesLog(30, 1.0e2, 1.0e5); //This is Daniel's binning in energy: binse->SetEdgesLog(14, 296.296, 86497.6); plist->AddToList(binse); //-------------------------------------------- MBinning *binssize = new MBinning("BinningSize"); binssize->SetEdgesLog(50, 10, 1.0e5); plist->AddToList(binssize); MBinning *binsdistc = new MBinning("BinningDist"); binsdistc->SetEdges(50, 0, 1.4); plist->AddToList(binsdistc); MBinning *binswidth = new MBinning("BinningWidth"); binswidth->SetEdges(50, 0, 1.0); plist->AddToList(binswidth); MBinning *binslength = new MBinning("BinningLength"); binslength->SetEdges(50, 0, 1.0); plist->AddToList(binslength); MBinning *binsalpha = new MBinning("BinningAlpha"); binsalpha->SetEdges(100, -100, 100); plist->AddToList(binsalpha); MBinning *binsasym = new MBinning("BinningAsym"); binsasym->SetEdges(50, -1.5, 1.5); plist->AddToList(binsasym); MBinning *binsm3l = new MBinning("BinningM3Long"); binsm3l->SetEdges(50, -1.5, 1.5); plist->AddToList(binsm3l); MBinning *binsm3t = new MBinning("BinningM3Trans"); binsm3t->SetEdges(50, -1.5, 1.5); plist->AddToList(binsm3t); //..... MBinning *binsb = new MBinning("BinningSigmabar"); binsb->SetEdges( 100, 0.0, 5.0); plist->AddToList(binsb); MBinning *binth = new MBinning("BinningTheta"); // this is Daniel's binning in theta //Double_t yedge[8] = // {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99}; // this is our binning Double_t yedge[9] = {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; TArrayD yed; yed.Set(9,yedge); binth->SetEdges(yed); plist->AddToList(binth); MBinning *bincosth = new MBinning("BinningCosTheta"); Double_t zedge[9]; for (Int_t i=0; i<9; i++) { zedge[8-i] = cos(yedge[i]/kRad2Deg); } TArrayD zed; zed.Set(9,zedge); bincosth->SetEdges(zed); plist->AddToList(bincosth); MBinning *binsdiff = new MBinning("BinningDiffsigma2"); binsdiff->SetEdges(100, -5.0, 20.0); plist->AddToList(binsdiff); // robert ---------------------------------------------- MBinning *binsalphaf = new MBinning("BinningAlphaFlux"); binsalphaf->SetEdges(100, -100, 100); plist->AddToList(binsalphaf); MBinning *binsdifftime = new MBinning("BinningTimeDiff"); binsdifftime->SetEdges(50, 0, 10); plist->AddToList(binsdifftime); MBinning *binstime = new MBinning("BinningTime"); binstime->SetEdges(50, 44500, 61000); plist->AddToList(binstime); } void DeleteBinnings(MParList *plist) { gLog << "DeleteBinnings" << endl; TObject *bin; //-------------------------------------------- bin = plist->FindObject("BinningE"); if (bin) delete bin; //-------------------------------------------- bin = plist->FindObject("BinningSize"); if (bin) delete bin; bin = plist->FindObject("BinningDist"); if (bin) delete bin; bin = plist->FindObject("BinningWidth"); if (bin) delete bin; bin = plist->FindObject("BinningLength"); if (bin) delete bin; bin = plist->FindObject("BinningAlpha"); if (bin) delete bin; bin = plist->FindObject("BinningAsym"); if (bin) delete bin; bin = plist->FindObject("BinningM3Long"); if (bin) delete bin; bin = plist->FindObject("BinningM3Trans"); if (bin) delete bin; //..... bin = plist->FindObject("BinningSigmabar"); if (bin) delete bin; bin = plist->FindObject("BinningTheta"); if (bin) delete bin; bin = plist->FindObject("BinningCosTheta"); if (bin) delete bin; bin = plist->FindObject("BinningDiffsigma2"); if (bin) delete bin; // robert ---------------------------------------------- bin = plist->FindObject("BinningAlphaFlux"); if (bin) delete bin; bin = plist->FindObject("BinningTimeDiff"); if (bin) delete bin; bin = plist->FindObject("BinningTime"); if (bin) delete bin; } //************************************************************************ void ONOFFCT1Analysis() { gLog.SetNoColors(); if (gRandom) delete gRandom; gRandom = new TRandom3(0); //----------------------------------------------- const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc"; const char *onfile = "~magican/ct1test/wittek/mkn421_00-01"; const char *mcfile = "~magican/ct1test/wittek/mkn421_mc_pedrms_0.001.preproc"; //----------------------------------------------- // path for input for Mars TString inPath = "~magican/ct1test/wittek/marsoutput/optionD/"; // path for output from Mars TString outPath = "~magican/ct1test/wittek/marsoutput/optionD/"; //----------------------------------------------- //TEnv env("macros/CT1env.rc"); //Bool_t printEnv = kFALSE; //************************************************************************ // Job A : // - produce MHSigmaTheta plots for ON and OFF data // - write out (or read in) these MHSigmaTheta plots // - read ON (or OFF or MC) data // - pad the events; // - write root file for ON (or OFF or MC) data (ON1.root, ...); Bool_t JobA = kFALSE; Bool_t WPad = kFALSE; // write out padding histograms ? Bool_t RPad = kFALSE; // read in padding histograms ? Bool_t Wout = kFALSE; // write out root file ON1.root // (or OFF1.root or MC1.root)? // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas // - if CTrainRF = TRUE : create matrices of training events // - if RTree = TRUE : read in trees, otherwise create trees // - calculate hadroness for method of RANDOM FOREST // - update the input files with the hadronesses (ON2.root, OFF2.root // or MC2.root) Bool_t JobB_RF_UP = kFALSE; Bool_t CTrainRF = kFALSE; // create matrices of training events Bool_t RTrainRF = kFALSE; // read in matrices of training events Bool_t RTree = kFALSE; // read in trees (otherwise grow them) Bool_t WRF = kFALSE; // update input root file ? // Job B_SC_UP : read ON2.root (or MC2.root) file // - depending on WParSC : create (or read in) supercuts parameter values // - calculate hadroness for the SUPERCUTS // - update the input files with the hadroness (==>ON3.root or MC3.root) Bool_t JobB_SC_UP = kFALSE; Bool_t CMatrix = kFALSE; // create training and test matrices Bool_t RMatrix = kFALSE; // read training and test matrices from file Bool_t WOptimize = kFALSE; // do optimization using the training sample // and write supercuts parameter values // onto the file parSCfile Bool_t RTest = kFALSE; // test the supercuts using the test matrix Bool_t WSC = kFALSE; // update input root file ? // Job C: // - read ON3.root and MC3.root files // which should have been updated to contain the hadronnesses // for the method of // RF // SUPERCUTS and // - produce Neyman-Pearson plots Bool_t JobC = kFALSE; // Job D : // - select g/h separation method XX // - read ON3 (or MC3) root file // - apply cuts in hadronness // - make plots Bool_t JobD = kFALSE; // Job E_XX : extended version of E_XX (including flux plots) // - select g/h separation method XX // - read MC root file   // - calculate eff. collection area // - optimize energy estimator // - read ON root file // - apply final cuts // - calculate flux // - write root file for ON data after final cuts Bool_t JobE_XX = kTRUE; Bool_t OEEst = kFALSE; // optimize energy estimator Bool_t WEX = kTRUE; // update root file ? Bool_t WRobert = kTRUE; // write out Robert's file ? //************************************************************************ //--------------------------------------------------------------------- // Job A //========= // read ON data file // - produce the 2D-histogram "sigmabar versus Theta" // (SigmaTheta_ON.root) for ON data // (to be used for the padding of the MC gamma data) // - write a file of ON events (ON1.root) // (after the standard cuts, before the g/h separation) // (to be used together with the corresponding MC gamma file (MC1.root) // for the optimization of the g/h separation) if (JobA) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job A" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobA, WPad, RPad, Wout = " << (JobA ? "kTRUE" : "kFALSE") << ", " << (WPad ? "kTRUE" : "kFALSE") << ", " << (RPad ? "kTRUE" : "kFALSE") << ", " << (Wout ? "kTRUE" : "kFALSE") << endl; //-------------------------------------------------- // names of ON and OFF files to be read // for generating the histograms to be used in the padding TString fileON = onfile; TString fileOFF = offfile; gLog << "fileON, fileOFF = " << fileON << ", " << fileOFF << endl; // name of file to conatin the histograms for the padding TString outNameSigTh = outPath; outNameSigTh += "SigmaTheta"; outNameSigTh += ".root"; //-------------------------------------------------- // type of data to be padded //TString typeInput = "ON"; //TString typeInput = "OFF"; TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; // name of input root file if (typeInput == "ON") TString filenamein(onfile); else if (typeInput == "OFF") TString filenamein(offfile); else if (typeInput == "MC") TString filenamein(mcfile); gLog << "data to be padded : " << filenamein << endl; // name of output root file TString outNameImage = outPath; outNameImage += typeInput; outNameImage += "1.root"; gLog << "padded data to be written onto : " << outNameImage << endl; //-------------------------------------------------- //************************************************************ // generate histograms to be used in the padding // // read ON and OFF data files // generate (or read in) the padding histograms for ON and OFF data // and merge these histograms MCT1PadONOFF pad; pad.SetName("MCT1PadONOFF"); pad.SetPadFlag(1); pad.SetDataType(typeInput); // generate the padding histograms if (!RPad) { gLog << "=====================================================" << endl; gLog << "Start generating the padding histograms" << endl; //----------------------------------------- // ON events gLog << "-----------" << endl; gLog << "ON events :" << endl; gLog << "-----------" << endl; MTaskList tliston; MParList pliston; MCT1ReadPreProc readON(fileON); //MFCT1SelBasic selthetaon; //selthetaon.SetCuts(-100.0, 29.5, 35.5); //MContinue contthetaon(&selthetaon); MBlindPixelCalc blindon; blindon.SetUseBlindPixels(); MFCT1SelBasic selbasicon; MContinue contbasicon(&selbasicon); MHBlindPixels blindON("BlindPixelsON"); MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels"); fillblindON.SetName("FillBlind"); MSigmabarCalc sigbarcalcon; MHSigmaTheta sigthON("SigmaThetaON"); MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt"); fillsigthetaON.SetName("FillSigTheta"); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&blindON); pliston.AddToList(&sigthON); //***************************** // entries in MTaskList tliston.AddToList(&readON); //tliston.AddToList(&contthetaon); tliston.AddToList(&blindon); tliston.AddToList(&contbasicon); tliston.AddToList(&fillblindON); tliston.AddToList(&sigbarcalcon); tliston.AddToList(&fillsigthetaON); MProgressBar baron; MEvtLoop evtloopon; evtloopon.SetParList(&pliston); evtloopon.SetProgressBar(&baron); Int_t maxeventson = -1; //Int_t maxeventson = 10000; if ( !evtloopon.Eventloop(maxeventson) ) return; tliston.PrintStatistics(0, kTRUE); blindON.DrawClone(); sigthON.DrawClone(); // save the histograms for the padding TH2D *sigthon = sigthON.GetSigmaTheta(); TH3D *sigpixthon = sigthON.GetSigmaPixTheta(); TH3D *diffpixthon = sigthON.GetDiffPixTheta(); TH2D *blindidthon = blindON.GetBlindId(); TH2D *blindnthon = blindON.GetBlindN(); //----------------------------------------- // OFF events gLog << "------------" << endl; gLog << "OFF events :" << endl; gLog << "------------" << endl; MTaskList tlistoff; MParList plistoff; MCT1ReadPreProc readOFF(fileOFF); MFCT1SelBasic selthetaoff; selthetaoff.SetCuts(-100.0, 29.5, 35.5); MContinue contthetaoff(&selthetaoff); MBlindPixelCalc blindoff; blindoff.SetUseBlindPixels(); MFCT1SelBasic selbasicoff; MContinue contbasicoff(&selbasicoff); MHBlindPixels blindOFF("BlindPixelsOFF"); MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels"); fillblindOFF.SetName("FillBlindOFF"); MSigmabarCalc sigbarcalcoff; MHSigmaTheta sigthOFF("SigmaThetaOFF"); MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt"); fillsigthetaOFF.SetName("FillSigThetaOFF"); //***************************** // entries in MParList plistoff.AddToList(&tlistoff); InitBinnings(&plistoff); plistoff.AddToList(&blindOFF); plistoff.AddToList(&sigthOFF); //***************************** // entries in MTaskList tlistoff.AddToList(&readOFF); //tlistoff.AddToList(&contthetaoff); tlistoff.AddToList(&blindoff); tlistoff.AddToList(&contbasicoff); tlistoff.AddToList(&fillblindOFF); tlistoff.AddToList(&sigbarcalcoff); tlistoff.AddToList(&fillsigthetaOFF); MProgressBar baroff; MEvtLoop evtloopoff; evtloopoff.SetParList(&plistoff); evtloopoff.SetProgressBar(&baroff); Int_t maxeventsoff = -1; //Int_t maxeventsoff = 20000; if ( !evtloopoff.Eventloop(maxeventsoff) ) return; tlistoff.PrintStatistics(0, kTRUE); blindOFF.DrawClone(); sigthOFF.DrawClone(); // save the histograms for the padding TH2D *sigthoff = sigthOFF.GetSigmaTheta(); TH3D *sigpixthoff = sigthOFF.GetSigmaPixTheta(); TH3D *diffpixthoff = sigthOFF.GetDiffPixTheta(); TH2D *blindidthoff = blindOFF.GetBlindId(); TH2D *blindnthoff = blindOFF.GetBlindN(); //----------------------------------------- gLog << "End of generating the padding histograms" << endl; gLog << "=====================================================" << endl; pad.MergeHistograms(sigthon, sigthoff, sigpixthon, sigpixthoff, diffpixthon, diffpixthoff, blindidthon, blindidthoff, blindnthon, blindnthoff); if (WPad) { // write the padding histograms onto a file --------- pad.WriteTargetDist(outNameSigTh); } } // read the padding histograms --------------------------- if (RPad) { pad.ReadTargetDist(outNameSigTh); } //************************************************************ if (Wout) { gLog << "=====================================================" << endl; gLog << "Start the padding" << endl; //----------------------------------------------------------- MTaskList tliston; MParList pliston; char *sourceName = "MSrcPosCam"; MSrcPosCam source(sourceName); // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MCT1ReadPreProc read(filenamein); MFCT1SelBasic seltheta; seltheta.SetCuts(-100.0, 29.5, 35.5); MContinue conttheta(&seltheta); if (typeInput == "ON") { MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", "MCT1PointingCorrCalc"); } MBlindPixelCalc blindbeforepad; blindbeforepad.SetUseBlindPixels(); blindbeforepad.SetName("BlindBeforePadding"); MBlindPixelCalc blind; blind.SetUseBlindPixels(); blind.SetName("BlindAfterPadding"); MFCT1SelBasic selbasic; MContinue contbasic(&selbasic); contbasic.SetName("SelBasic"); MFillH fillblind("BlindPixels[MHBlindPixels]", "MBlindPixels"); fillblind.SetName("HBlind"); MSigmabarCalc sigbarcalc; MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); fillsigtheta.SetName("HSigmaTheta"); MImgCleanStd clean; // calculation of image parameters --------------------- TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; MHillasCalc hcalc; hcalc.SetNameHillas(fHilName); hcalc.SetNameHillasExt(fHilNameExt); hcalc.SetNameNewImgPar(fImgParName); MHillasSrcCalc hsrccalc(sourceName, fHilNameSrc); hsrccalc.SetInput(fHilName); MFillH hfill1("MHHillas", fHilName); hfill1.SetName("HHillas"); MFillH hfill2("MHStarMap", fHilName); hfill2.SetName("HStarMap"); MFillH hfill3("MHHillasExt", fHilNameSrc); hfill3.SetName("HHillasExt"); MFillH hfill4("MHHillasSrc", fHilNameSrc); hfill4.SetName("HHillasSrc"); MFillH hfill5("MHNewImagePar", fImgParName); hfill5.SetName("HNewImagePar"); // -------------------------------------------------- MFCT1SelStandard selstandard(fHilNameSrc); selstandard.SetHillasName(fHilName); selstandard.SetImgParName(fImgParName); selstandard.SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0); MContinue contstandard(&selstandard); contstandard.SetName("SelStandard"); MWriteRootFile write(outNameImage); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); write.AddContainer("MMcEvt", "Events"); write.AddContainer("ThetaOrig", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasExt", "Events"); write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&source); //***************************** // entries in MTaskList tliston.AddToList(&read); //tliston.AddToList(&conttheta); tliston.AddToList(&blindbeforepad); tliston.AddToList(&pad); if (typeInput == "ON") tliston.AddToList(&pointcorr); tliston.AddToList(&blind); tliston.AddToList(&contbasic); tliston.AddToList(&fillblind); tliston.AddToList(&sigbarcalc); tliston.AddToList(&fillsigtheta); tliston.AddToList(&clean); tliston.AddToList(&hcalc); tliston.AddToList(&hsrccalc); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); tliston.AddToList(&contstandard); tliston.AddToList(&write); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&pliston); //evtloop.ReadEnv(env, "", printEnv); evtloop.SetProgressBar(&bar); // evtloop.Write(); Int_t maxevents = -1; //Int_t maxevents = 1000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); //------------------------------------------- // Display the histograms pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone(); pliston.FindObject("BlindPixels", "MHBlindPixels")->DrawClone(); pliston.FindObject("MHHillas")->DrawClone(); pliston.FindObject("MHHillasExt")->DrawClone(); pliston.FindObject("MHHillasSrc")->DrawClone(); pliston.FindObject("MHNewImagePar")->DrawClone(); pliston.FindObject("MHStarMap")->DrawClone(); DeleteBinnings(&pliston); gLog << "End of padding" << endl; gLog << "=====================================================" << endl; } gLog << "Macro CT1Analysis : End of Job A" << endl; gLog << "===================================================" << endl; } //--------------------------------------------------------------------- // Job B_RF_UP //============ // - create (or read in) the matrices of training events for gammas // and hadrons // - create (or read in) the trees // - then read ON1.root (or MC1.root) file // - calculate the hadroness for the method of RANDOM FOREST // - update input root file with the hadroness if (JobB_RF_UP) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = " << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", " << (RTrainRF ? "kTRUE" : "kFALSE") << ", " << (CTrainRF ? "kTRUE" : "kFALSE") << ", " << (RTree ? "kTRUE" : "kFALSE") << ", " << (WRF ? "kTRUE" : "kFALSE") << endl; //-------------------------------------------- // parameters for the random forest Int_t NumTrees = 100; Int_t NumTry = 3; Int_t NdSize = 1; TString hadRFName = "HadRF"; Float_t maxhadronness = 0.23; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; //-------------------------------------------- // file to be updated (ON, OFF or MC) //TString typeInput = "ON"; //TString typeInput = "OFF"; TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; // name of input root file TString filenameData = outPath; filenameData += typeInput; filenameData += "1.root"; gLog << "filenameData = " << filenameData << endl; // name of output root file TString outNameImage = outPath; outNameImage += typeInput; outNameImage += "2.root"; //TString outNameImage = filenameData; gLog << "outNameImage = " << outNameImage << endl; //-------------------------------------------- // files to be read for generating the matrices of training events // "hadrons" : TString filenameHad = outPath; filenameHad += "OFF"; filenameHad += "1.root"; Int_t howManyHadrons = 1000000; gLog << "filenameHad = " << filenameHad << ", howManyHadrons = " << howManyHadrons << endl; // "gammas" : TString filenameMC = outPath; filenameMC += "MC"; filenameMC += "1.root"; Int_t howManyGammas = 50000; gLog << "filenameMC = " << filenameMC << ", howManyGammas = " << howManyGammas << endl; //-------------------------------------------- // files of training events TString outNameGammas = outPath; outNameGammas += "RFmatrix_gammas_"; outNameGammas += "MC"; outNameGammas += ".root"; TString typeMatrixHadrons = "OFF"; gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; TString outNameHadrons = outPath; outNameHadrons += "RFmatrix_hadrons_"; outNameHadrons += typeMatrixHadrons; outNameHadrons += ".root"; MHMatrix matrixg("MatrixGammas"); matrixg.EnableGraphicalOutput(); matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)"); matrixg.AddColumn("MSigmabar.fSigmabar"); matrixg.AddColumn("log10(MHillas.fSize)"); matrixg.AddColumn("MHillasSrc.fDist"); matrixg.AddColumn("MHillas.fWidth"); matrixg.AddColumn("MHillas.fLength"); matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); matrixg.AddColumn("MNewImagePar.fConc"); matrixg.AddColumn("MNewImagePar.fLeakage1"); MHMatrix matrixh("MatrixHadrons"); matrixh.EnableGraphicalOutput(); matrixh.AddColumns(matrixg.GetColumns()); //-------------------------------------------- // file of trees of the random forest TString outRF = outPath; outRF += "RF.root"; //************************************************************************* // read in matrices of training events if (RTrainRF) { const char* mtxName = "MatrixGammas"; gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Get matrix for (gammas)" << endl; gLog << "matrix name = " << mtxName << endl; gLog << "name of root file = " << outNameGammas << endl; gLog << "" << endl; // read in the object with the name 'mtxName' from file 'outNameGammas' // TFile fileg(outNameGammas); matrixg.Read(mtxName); matrixg.Print("SizeCols"); //***************************************************************** const char* mtxName = "MatrixHadrons"; gLog << "" << endl; gLog << "========================================================" << endl; gLog << " Get matrix for (hadrons)" << endl; gLog << "matrix name = " << mtxName << endl; gLog << "name of root file = " << outNameHadrons << endl; gLog << "" << endl; // read in the object with the name 'mtxName' from file 'outNameHadrons' // TFile fileh(outNameHadrons); matrixh.Read(mtxName); matrixh.Print("SizeCols"); } //************************************************************************* // create matrices of training events if (CTrainRF) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << " Create matrices of training events" << endl; gLog << " Gammas :" << endl; MParList plistg; MTaskList tlistg; MFilterList flistg; MParList plisth; MTaskList tlisth; MFilterList flisth; MReadMarsFile readg("Events", filenameMC); readg.DisableAutoScheme(); MReadMarsFile readh("Events", filenameHad); readh.DisableAutoScheme(); MFParticleId fgamma("MMcEvt", '=', kGAMMA); fgamma.SetName("gammaID"); MFParticleId fhadrons("MMcEvt", '!', kGAMMA); fhadrons.SetName("hadronID)"); TString mgname("costhg"); MBinning bing("Binning"+mgname); bing.SetEdges(10, 0., 1.0); MH3 gref("cos(MMcEvt.fTelescopeTheta)"); gref.SetName(mgname); MH::SetBinning(&gref.GetHist(), &bing); for (Int_t i=1; i<=gref.GetNbins(); i++) gref.GetHist().SetBinContent(i, 1.0); MFEventSelector2 selectorg(gref); selectorg.SetNumMax(howManyGammas); selectorg.SetName("selectGammas"); MFillH fillmatg("MatrixGammas"); fillmatg.SetFilter(&flistg); fillmatg.SetName("fillGammas"); //***************************** fill gammas *** // entries in MFilterList flistg.AddToList(&fgamma); flistg.AddToList(&selectorg); //***************************** // entries in MParList plistg.AddToList(&tlistg); InitBinnings(&plistg); plistg.AddToList(&matrixg); //***************************** // entries in MTaskList tlistg.AddToList(&readg); tlistg.AddToList(&flistg); tlistg.AddToList(&fillmatg); //***************************** MProgressBar matrixbar; MEvtLoop evtloopg; evtloopg.SetName("FillGammaMatrix"); evtloopg.SetParList(&plistg); //evtloopg.ReadEnv(env, "", printEnv); evtloopg.SetProgressBar(&matrixbar); Int_t maxevents = -1; if (!evtloopg.Eventloop(maxevents)) return; tlistg.PrintStatistics(0, kTRUE); //***************************** fill hadrons *** gLog << " Hadrons :" << endl; TString mhname("costhh"); MBinning binh("Binning"+mhname); binh.SetEdges(10, 0., 1.0); MH3 href("cos(MMcEvt.fTelescopeTheta)"); href.SetName(mhname); MH::SetBinning(&href.GetHist(), &binh); for (Int_t i=1; i<=href.GetNbins(); i++) href.GetHist().SetBinContent(i, 1.0); MFEventSelector2 selectorh(href); //selectorh.SetNumMax(howManyHadrons); // select as many hadrons as gammas selectorh.SetNumMax(matrixg.GetM().GetNrows()); selectorh.SetName("selectHadrons"); MFillH fillmath("MatrixHadrons"); fillmath.SetFilter(&flisth); fillmath.SetName("fillHadrons"); // entries in MFilterList flisth.AddToList(&fhadrons); flisth.AddToList(&selectorh); //***************************** // entries in MParList plisth.AddToList(&tlisth); InitBinnings(&plisth); plisth.AddToList(&matrixh); //***************************** // entries in MTaskList tlisth.AddToList(&readh); tlisth.AddToList(&flisth); tlisth.AddToList(&fillmath); //***************************** MProgressBar matrixbar; MEvtLoop evtlooph; evtlooph.SetName("FillHadronMatrix"); evtlooph.SetParList(&plisth); //evtlooph.ReadEnv(env, "", printEnv); evtlooph.SetProgressBar(&matrixbar); Int_t maxevents = -1; if (!evtlooph.Eventloop(maxevents)) return; tlisth.PrintStatistics(0, kTRUE); //***************************************************** // write out matrices of training events gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Write out matrices of training events" << endl; //------------------------------------------- // "gammas" gLog << "Gammas :" << endl; matrixg.Print("SizeCols"); TFile writeg(outNameGammas, "RECREATE", ""); matrixg.Write(); gLog << "" << endl; gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file " << outNameGammas << endl; //------------------------------------------- // "hadrons" gLog << "Hadrons :" << endl; matrixh.Print("SizeCols"); TFile writeh(outNameHadrons, "RECREATE", ""); matrixh.Write(); gLog << "" << endl; gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file " << outNameHadrons << endl; } //********** end of creating matrices of training events *********** MRanForest *fRanForest; MRanTree *fRanTree; //----------------------------------------------------------------- // read in the trees of the random forest if (RTree) { MParList plisttr; MTaskList tlisttr; plisttr.AddToList(&tlisttr); MReadTree readtr("TREE", outRF); readtr.DisableAutoScheme(); MRanForestFill rffill; rffill.SetNumTrees(NumTrees); // list of tasks for the loop over the trees tlisttr.AddToList(&readtr); tlisttr.AddToList(&rffill); //------------------- // Execute tree loop // MEvtLoop evtlooptr; evtlooptr.SetName("ReadRFTrees"); evtlooptr.SetParList(&plisttr); if (!evtlooptr.Eventloop()) return; tlisttr.PrintStatistics(0, kTRUE); // get adresses of objects which are used in the next eventloop fRanForest = (MRanForest*)plisttr->FindObject("MRanForest"); if (!fRanForest) { gLog << err << dbginf << "MRanForest not found... aborting." << endl; return kFALSE; } fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); if (!fRanTree) { gLog << err << dbginf << "MRanTree not found... aborting." << endl; return kFALSE; } } //----------------------------------------------------------------- // grow the trees of the random forest (event loop = tree loop) if (!RTree) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Macro CT1Analysis : start growing trees" << endl; MTaskList tlist2; MParList plist2; plist2.AddToList(&tlist2); plist2.AddToList(&matrixg); plist2.AddToList(&matrixh); MRanForestGrow rfgrow2; rfgrow2.SetNumTrees(NumTrees); rfgrow2.SetNumTry(NumTry); rfgrow2.SetNdSize(NdSize); MWriteRootFile rfwrite2(outRF); rfwrite2.AddContainer("MRanTree", "TREE"); MFillH fillh2("MHRanForestGini"); // list of tasks for the loop over the trees tlist2.AddToList(&rfgrow2); tlist2.AddToList(&rfwrite2); tlist2.AddToList(&fillh2); //------------------- // Execute tree loop // MEvtLoop treeloop; treeloop.SetName("GrowRFTrees"); treeloop.SetParList(&plist2); if ( !treeloop.Eventloop() ) return; tlist2.PrintStatistics(0, kTRUE); plist2.FindObject("MHRanForestGini")->DrawClone(); // get adresses of objects which are used in the next eventloop fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); if (!fRanForest) { gLog << err << dbginf << "MRanForest not found... aborting." << endl; return kFALSE; } fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); if (!fRanTree) { gLog << err << dbginf << "MRanTree not found... aborting." << endl; return kFALSE; } } // end of growing the trees of the random forest //----------------------------------------------------------------- //----------------------------------------------------------------- // Update the input files with the RF hadronness // if (WRF) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Update input file '" << filenameData << "' with the RF hadronness" << endl; MTaskList tliston; MParList pliston; // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); //....................................................................... // calculate hadronnes for method of RANDOM FOREST MRanForestCalc rfcalc; rfcalc.SetHadronnessName(hadRFName); //....................................................................... //MWriteRootFile write(outNameImage, "UPDATE"); MWriteRootFile write(outNameImage, "RECREATE"); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); write.AddContainer("MMcEvt", "Events"); write.AddContainer("ThetaOrig", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasExt", "Events"); write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); write.AddContainer(hadRFName, "Events"); //----------------------------------------------------------------- MFCT1SelFinal selfinalgh(fHilNameSrc); selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); selfinalgh.SetHadronnessName(hadRFName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); MFillH fillranfor("MHRanForest"); fillranfor.SetName("HRanForest"); MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); fillhadrf.SetName("HhadRF"); MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(hadRFName); selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); contfinal.SetName("ContSelFinal"); TString mh3name = "abs(Alpha)"; MBinning binsalphaabs("Binning"+mh3name); binsalphaabs.SetEdges(50, -2.0, 98.0); MH3 alphaabs("abs(MHillasSrc.fAlpha)"); alphaabs.SetName(mh3name); MFillH alpha(&alphaabs); alpha.SetName("FillAlphaAbs"); MFillH hfill1("MHHillas", fHilName); hfill1.SetName("HHillas"); MFillH hfill2("MHStarMap", fHilName); hfill2.SetName("HStarMap"); MFillH hfill3("MHHillasExt", fHilNameSrc); hfill3.SetName("HHillasExt"); MFillH hfill4("MHHillasSrc", fHilNameSrc); hfill4.SetName("HHillasSrc"); MFillH hfill5("MHNewImagePar", fImgParName); hfill5.SetName("HNewImagePar"); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(fRanForest); pliston.AddToList(fRanTree); pliston.AddToList(&binsalphaabs); pliston.AddToList(&alphaabs); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&rfcalc); tliston.AddToList(&fillranfor); tliston.AddToList(&fillhadrf); tliston.AddToList(&write); tliston.AddToList(&contfinalgh); tliston.AddToList(&alpha); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); tliston.AddToList(&contfinal); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetName("UpdateRootFile"); evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); Int_t maxevents = -1; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); //------------------------------------------- // Display the histograms // pliston.FindObject("MHRanForest")->DrawClone(); pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); pliston.FindObject("MHHillas")->DrawClone(); pliston.FindObject("MHHillasExt")->DrawClone(); pliston.FindObject("MHHillasSrc")->DrawClone(); pliston.FindObject("MHNewImagePar")->DrawClone(); pliston.FindObject("MHStarMap")->DrawClone(); //------------------------------------------- // fit alpha distribution to get the number of excess events and // calculate significance of gamma signal in the alpha plot MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); TH1 &alphaHist = absalpha->GetHist(); alphaHist.SetXTitle("|alpha| [\\circ]"); alphaHist.SetName("alpha-macro"); Double_t alphasig = 13.1; Double_t alphamin = 30.0; Double_t alphamax = 90.0; Int_t degree = 2; Double_t significance = -99.0; Bool_t drawpoly = kTRUE; Bool_t fitgauss = kTRUE; Bool_t print = kTRUE; MHFindSignificance findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, alphasig, drawpoly, fitgauss, print); significance = findsig.GetSignificance(); Float_t alphasi = findsig.GetAlphasi(); gLog << "For file '" << filenameData << "' : " << endl; gLog << "Significance of gamma signal after supercuts : " << significance << " (for |alpha| < " << alphasi << " degrees)" << endl; findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); //------------------------------------------- DeleteBinnings(&pliston); } gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- // Job B_SC_UP //============ // - create (or read in) optimum supercuts parameter values // // - calculate the hadroness for the supercuts // // - update input root file, including the hadroness if (JobB_SC_UP) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job B_SC_UP" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = " << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", " << (CMatrix ? "kTRUE" : "kFALSE") << ", " << (RMatrix ? "kTRUE" : "kFALSE") << ", " << (WOptimize ? "kTRUE" : "kFALSE") << ", " << (RTest ? "kTRUE" : "kFALSE") << ", " << (WSC ? "kTRUE" : "kFALSE") << endl; //-------------------------------------------- // file which contains the initial parameter values for the supercuts // if parSCinit ="" the initial values are taken from the constructor of // MCT1Supercuts TString parSCinit = outPath; //parSCinit += "parSC_1709d"; parSCinit = ""; gLog << "parSCinit = " << parSCinit << endl; //--------------- // file onto which the optimal parameter values for the supercuts // are written TString parSCfile = outPath; parSCfile += "parSC_2310a"; gLog << "parSCfile = " << parSCfile << endl; //-------------------------------------------- // file to be updated (either ON or MC) //TString typeInput = "ON"; //TString typeInput = "OFF"; TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; // name of input root file TString filenameData = outPath; filenameData += typeInput; filenameData += "2.root"; gLog << "filenameData = " << filenameData << endl; // name of output root file TString outNameImage = outPath; outNameImage += typeInput; outNameImage += "3.root"; //TString outNameImage = filenameData; gLog << "outNameImage = " << outNameImage << endl; //-------------------------------------------- // files to be read for optimizing the supercuts // // for the training TString filenameTrain = outPath; filenameTrain += "ON"; filenameTrain += "1.root"; Int_t howManyTrain = 800000; gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = " << howManyTrain << endl; // for testing TString filenameTest = outPath; filenameTest += "ON"; filenameTest += "1.root"; Int_t howManyTest = 800000; gLog << "filenameTest = " << filenameTest << ", howManyTest = " << howManyTest << endl; //-------------------------------------------- // files to contain the matrices (generated from filenameTrain and // filenameTest) // // for the training TString fileMatrixTrain = outPath; fileMatrixTrain += "MatrixTrainSC"; fileMatrixTrain += ".root"; gLog << "fileMatrixTrain = " << fileMatrixTrain << endl; // for testing TString fileMatrixTest = outPath; fileMatrixTest += "MatrixTestSC"; fileMatrixTest += ".root"; gLog << "fileMatrixTest = " << fileMatrixTest << endl; //--------------------------------------------------------------------- // Training and test matrices : // - either create them and write them onto a file // - or read them from a file MCT1FindSupercuts findsuper; findsuper.SetFilenameParam(parSCfile); findsuper.SetHadronnessName("HadSC"); //-------------------------- // create matrices and write them onto files if (CMatrix) { MH3 &mh3 = *(new MH3("MHillas.fSize")); mh3.SetName("Target distribution for SIZE"); if (filenameTrain == filenameTest) { if ( !findsuper.DefineTrainTestMatrix(filenameTrain, howManyTrain, mh3, howManyTest, mh3, fileMatrixTrain, fileMatrixTest) ) { *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl; return; } } else { if ( !findsuper.DefineTrainMatrix(filenameTrain, howManyTrain, mh3, fileMatrixTrain) ) { *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl; return; } if ( !findsuper.DefineTestMatrix( filenameTest, howManyTest, mh3, fileMatrixTest) ) { *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl; return; } } } //-------------------------- // read matrices from files // if (RMatrix) findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest); //-------------------------- //--------------------------------------------------------------------- // optimize supercuts using the training sample // // the initial values are taken // - from the file parSCinit (if != "") // - or from the arrays params and steps (if their sizes are != 0) // - or from the MCT1Supercuts constructor if (WOptimize) { gLog << "CT1Analysis.C : optimize the supercuts using the training matrix" << endl; TArrayD params(0); TArrayD steps(0); if (parSCinit == "") { Double_t vparams[104] = { // LengthUp 0.315585, 0.001455, 0.203198, 0.005532, -0.001670, -0.020362, 0.007388, -0.013463, // LengthLo 0.151530, 0.028323, 0.510707, 0.053089, 0.013708, 2.357993, 0.000080, -0.007157, // WidthUp 0.145412, -0.001771, 0.054462, 0.022280, -0.009893, 0.056353, 0.020711, -0.016703, // WidthLo 0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101, 0.006126, -0.002849, // DistUp 1.787943, 0.0, 2.942310, 0.199815, 0.0, 0.249909, 0.189697, 0.0, // DistLo 0.589406, 0.0, -0.083964,-0.007975, 0.0, 0.045374, -0.001750, 0.0, // AsymUp 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // AsymLo -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // ConcUp 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // ConcLo -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // Leakage1Up 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // Leakage1Lo -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // AlphaUp 13.12344, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; Double_t vsteps[104] = { // LengthUp 0.03, 0.0002, 0.02, 0.0006, 0.0002, 0.002, 0.0008, 0.002, // LengthLo 0.02, 0.003, 0.05, 0.006, 0.002, 0.3, 0.0001, 0.0008, // WidthUp 0.02, 0.0002, 0.006, 0.003, 0.002, 0.006, 0.002, 0.002, // WidthLo 0.009, 0.0007, 0.008, 0.0004, 0.0005, 0.002, 0.0007, 0.003, // DistUp 0.2, 0.0, 0.3, 0.02, 0.0, 0.03, 0.02, 0.0 // DistLo 0.06, 0.0, 0.009, 0.0008, 0.0, 0.005, 0.0002, 0.0 // AsymUp 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // AsymLo 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // ConcUp 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // ConcLo 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // Leakage1Up 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // Leakage1Lo 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // AlphaUp 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; params.Set(104, vparams); steps.Set (104, vsteps ); } Bool_t rf; rf = findsuper.FindParams(parSCinit, params, steps); if (!rf) { gLog << "CT1Analysis.C : optimization of supercuts failed" << endl; return; } } //-------------------------------------- // test the supercuts on the test sample // if (RTest) { gLog << "CT1Analysis.C : test the supercuts on the test matrix" << endl; Bool_t rt = findsuper.TestParams(); if (!rt) { gLog << "CT1Analysis.C : test of supercuts on the test matrix failed" << endl; } } //----------------------------------------------------------------- // Update the input files with the SC hadronness // if (WSC) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Update input file '" << filenameData << "' with the SC hadronness" << endl; //---------------------------------------------------- // read in optimum parameter values for the supercuts TFile inparam(parSCfile); MCT1Supercuts scin; scin.Read("MCT1Supercuts"); inparam.Close(); gLog << "Parameter values for supercuts were read in from file '" << parSCfile << "'" << endl; TArrayD supercutsPar; supercutsPar = scin.GetParameters(); TArrayD supercutsStep; supercutsStep = scin.GetStepsizes(); gLog << "Parameter values for supercuts : " << endl; for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; //....................................................................... // calculation of hadroness for the supercuts // (=0.25 if fullfilled, =0.75 otherwise) TString hadSCName = "HadSC"; MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc); sccalc.SetHadronnessName(hadSCName); //....................................................................... //MWriteRootFile write(outNameImage, "UPDATE"); //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE"); MWriteRootFile write(outNameImage, "RECREATE"); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); write.AddContainer("MMcEvt", "Events"); write.AddContainer("ThetaOrig", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasExt", "Events"); write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); write.AddContainer("HadRF", "Events"); write.AddContainer(hadSCName, "Events"); //----------------------------------------------------------------- // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); Float_t maxhadronness = 0.40; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; MFCT1SelFinal selfinalgh(fHilNameSrc); selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); selfinalgh.SetHadronnessName(hadSCName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); fillhadsc.SetName("HhadSC"); MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(hadSCName); selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); contfinal.SetName("ContSelFinal"); TString mh3name = "abs(Alpha)"; MBinning binsalphaabs("Binning"+mh3name); binsalphaabs.SetEdges(50, -2.0, 98.0); MH3 alphaabs("abs(MHillasSrc.fAlpha)"); alphaabs.SetName(mh3name); TH1 &alphahist = alphaabs->GetHist(); MFillH alpha(&alphaabs); alpha.SetName("FillAlphaAbs"); MFillH hfill1("MHHillas", fHilName); hfill1.SetName("HHillas"); MFillH hfill2("MHStarMap", fHilName); hfill2.SetName("HStarMap"); MFillH hfill3("MHHillasExt", fHilNameSrc); hfill3.SetName("HHillasExt"); MFillH hfill4("MHHillasSrc", fHilNameSrc); hfill4.SetName("HHillasSrc"); MFillH hfill5("MHNewImagePar", fImgParName); hfill5.SetName("HNewImagePar"); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&supercuts); pliston.AddToList(&binsalphaabs); pliston.AddToList(&alphaabs); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&sccalc); tliston.AddToList(&fillhadsc); tliston.AddToList(&write); tliston.AddToList(&contfinalgh); tliston.AddToList(&alpha); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); tliston.AddToList(&contfinal); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); Int_t maxevents = -1; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); //------------------------------------------- // Display the histograms // pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); pliston.FindObject("MHHillas")->DrawClone(); pliston.FindObject("MHHillasExt")->DrawClone(); pliston.FindObject("MHHillasSrc")->DrawClone(); pliston.FindObject("MHNewImagePar")->DrawClone(); pliston.FindObject("MHStarMap")->DrawClone(); //------------------------------------------- // fit alpha distribution to get the number of excess events and // calculate significance of gamma signal in the alpha plot MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); TH1 &alphaHist = absalpha->GetHist(); alphaHist.SetXTitle("|alpha| [\\circ]"); alphaHist.SetName("alpha-macro"); Double_t alphasig = 13.1; Double_t alphamin = 30.0; Double_t alphamax = 90.0; Int_t degree = 2; Double_t significance = -99.0; Bool_t drawpoly = kTRUE; Bool_t fitgauss = kTRUE; Bool_t print = kTRUE; MHFindSignificance findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, alphasig, drawpoly, fitgauss, print); significance = findsig.GetSignificance(); Float_t alphasi = findsig.GetAlphasi(); gLog << "For file '" << filenameData << "' : " << endl; gLog << "Significance of gamma signal after supercuts : " << significance << " (for |alpha| < " << alphasi << " degrees)" << endl; findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); //------------------------------------------- DeleteBinnings(&pliston); } gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- // Job C //====== // - read ON1 and MC1 data files // which should have been updated to contain the hadronnesses // for the method of Random Forest and for the SUPERCUTS // - produce Neyman-Pearson plots if (JobC) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job C" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobC = " << (JobC ? "kTRUE" : "kFALSE") << endl; // name of input data file TString filenameData = outPath; filenameData += "OFF"; filenameData += "3.root"; gLog << "filenameData = " << filenameData << endl; // name of input MC file TString filenameMC = outPath; filenameMC += "MC"; filenameMC += "3.root"; gLog << "filenameMC = " << filenameMC << endl; //----------------------------------------------------------------- MTaskList tliston; MParList pliston; // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameMC); read.AddFile(filenameData); read.DisableAutoScheme(); //....................................................................... // names of hadronness containers TString hadSCName = "HadSC"; TString hadRFName = "HadRF"; //....................................................................... TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; Float_t maxhadronness = 0.40; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; MFCT1SelFinal selfinalgh(fHilNameSrc); selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); selfinalgh.SetHadronnessName(hadRFName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); fillhadsc.SetName("HhadSC"); MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); fillhadrf.SetName("HhadRF"); MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(hadRFName); selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); contfinal.SetName("ContSelFinal"); MFillH hfill1("MHHillas", fHilName); hfill1.SetName("HHillas"); MFillH hfill2("MHStarMap", fHilName); hfill2.SetName("HStarMap"); MFillH hfill3("MHHillasExt", fHilNameSrc); hfill3.SetName("HHillasExt"); MFillH hfill4("MHHillasSrc", fHilNameSrc); hfill4.SetName("HHillasSrc"); MFillH hfill5("MHNewImagePar", fImgParName); hfill5.SetName("HNewImagePar"); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&fillhadsc); tliston.AddToList(&fillhadrf); tliston.AddToList(&contfinalgh); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); tliston.AddToList(&contfinal); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); Int_t maxevents = -1; //Int_t maxevents = 35000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); //------------------------------------------- // Display the histograms // pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); pliston.FindObject("MHHillas")->DrawClone(); pliston.FindObject("MHHillasExt")->DrawClone(); pliston.FindObject("MHHillasSrc")->DrawClone(); pliston.FindObject("MHNewImagePar")->DrawClone(); pliston.FindObject("MHStarMap")->DrawClone(); DeleteBinnings(&pliston); gLog << "Macro CT1Analysis : End of Job C" << endl; gLog << "===================================================" << endl; } //--------------------------------------------------------------------- // Job D //====== // - select g/h separation method XX // - read ON2 (or MC2) root file // - apply cuts in hadronness // - make plots if (JobD) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job D" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobD = " << (JobD ? "kTRUE" : "kFALSE") << endl; // type of data to be analysed TString typeData = "ON"; //TString typeData = "OFF"; //TString typeData = "MC"; gLog << "typeData = " << typeData << endl; TString ext = "3.root"; //------------------------------ // selection of g/h separation method // and definition of final selections //TString XX("SC"); TString XX("RF"); TString fhadronnessName("Had"); fhadronnessName += XX; gLog << "fhadronnessName = " << fhadronnessName << endl; // maximum values of the hadronness, |ALPHA| and DIST Float_t maxhadronness = 0.233; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; gLog << "Maximum values of hadronness, |ALPHA| and DIST = " << maxhadronness << ", " << maxalpha << ", " << maxdist << endl; //------------------------------ // name of data file to be analysed TString filenameData(outPath); filenameData += typeData; filenameData += ext; gLog << "filenameData = " << filenameData << endl; //************************************************************************* // // Analyse the data // MTaskList tliston; MParList pliston; // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); //----------------------------------------------------------------- // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); MFCT1SelFinal selfinalgh(fHilNameSrc); selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); selfinalgh.SetHadronnessName(fhadronnessName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); fillhadsc.SetName("HhadSC"); MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); fillhadrf.SetName("HhadRF"); TString mh3name = "abs(Alpha)"; MBinning binsalphaabs("Binning"+mh3name); binsalphaabs.SetEdges(50, -2.0, 98.0); MH3 alphaabs("abs(MHillasSrc.fAlpha)"); alphaabs.SetName(mh3name); TH1 &alphahist = alphaabs->GetHist(); MFillH alpha(&alphaabs); alpha.SetName("FillAlphaAbs"); MFillH hfill1("MHHillas", fHilName); hfill1.SetName("HHillas"); MFillH hfill2("MHStarMap", fHilName); hfill2.SetName("HStarMap"); MFillH hfill3("MHHillasExt", fHilNameSrc); hfill3.SetName("HHillasExt"); MFillH hfill4("MHHillasSrc", fHilNameSrc); hfill4.SetName("HHillasSrc"); MFillH hfill5("MHNewImagePar", fImgParName); hfill5.SetName("HNewImagePar"); MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(fhadronnessName); selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); contfinal.SetName("ContSelFinal"); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&binsalphaabs); pliston.AddToList(&alphaabs); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&contfinalgh); tliston.AddToList(&fillhadsc); tliston.AddToList(&fillhadrf); tliston.AddToList(&alpha); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); tliston.AddToList(&contfinal); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); Int_t maxevents = -1; //Int_t maxevents = 10000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); //------------------------------------------- // Display the histograms // pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); pliston.FindObject("MHHillas")->DrawClone(); pliston.FindObject("MHHillasExt")->DrawClone(); pliston.FindObject("MHHillasSrc")->DrawClone(); pliston.FindObject("MHNewImagePar")->DrawClone(); pliston.FindObject("MHStarMap")->DrawClone(); //------------------------------------------- // fit alpha distribution to get the number of excess events and // calculate significance of gamma signal in the alpha plot MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); TH1 &alphaHist = absalpha->GetHist(); alphaHist.SetXTitle("|alpha| [\\circ]"); alphaHist.SetName("alpha-JobD"); Double_t alphasig = 13.1; Double_t alphamin = 30.0; Double_t alphamax = 90.0; Int_t degree = 2; Double_t significance = -99.0; Bool_t drawpoly = kTRUE; Bool_t fitgauss = kTRUE; Bool_t print = kTRUE; MHFindSignificance findsig; findsig.SetRebin(kTRUE); findsig.SetReduceDegree(kFALSE); findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, alphasig, drawpoly, fitgauss, print); significance = findsig.GetSignificance(); Float_t alphasi = findsig.GetAlphasi(); gLog << "For file '" << filenameData << "' : " << endl; gLog << "Significance of gamma signal after supercuts : " << significance << " (for |alpha| < " << alphasi << " degrees)" << endl; findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); //------------------------------------------- DeleteBinnings(&pliston); gLog << "Macro CT1Analysis : End of Job D" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- // Job E_XX //========= // - select g/h separation method XX // - read MC_XX2.root file   // - calculate eff. collection area // - read ON_XX2.root file // - apply final cuts // - calculate flux // - write root file for ON data after final cuts (ON_XX3.root)) if (JobE_XX) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job E_XX" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = " << (JobE_XX ? "kTRUE" : "kFALSE") << ", " << (OEEst ? "kTRUE" : "kFALSE") << ", " << (WEX ? "kTRUE" : "kFALSE") << endl; // type of data to be analysed //TString typeData = "ON"; //TString typeData = "OFF"; TString typeData = "MC"; gLog << "typeData = " << typeData << endl; TString typeMC = "MC"; TString ext = "3.root"; TString extout = "4.root"; //------------------------------ // selection of g/h separation method // and definition of final selections //TString XX("SC"); TString XX("RF"); TString fhadronnessName("Had"); fhadronnessName += XX; gLog << "fhadronnessName = " << fhadronnessName << endl; // maximum values of the hadronness, |ALPHA| and DIST Float_t maxhadronness = 0.23; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; gLog << "Maximum values of hadronness, |ALPHA| and DIST = " << maxhadronness << ", " << maxalpha << ", " << maxdist << endl; //------------------------------ // name of MC file to be used for optimizing the energy estimator TString filenameOpt(outPath); filenameOpt += typeMC; filenameOpt += ext; gLog << "filenameOpt = " << filenameOpt << endl; //------------------------------ // name of file containing the parameters of the energy estimator TString energyParName(outPath); energyParName += "energyest_"; energyParName += XX; energyParName += ".root"; gLog << "energyParName = " << energyParName << endl; //------------------------------ // name of MC file to be used for calculating the eff. collection areas TString filenameArea(outPath); filenameArea += typeMC; filenameArea += ext; gLog << "filenameArea = " << filenameArea << endl; //------------------------------ // name of file containing the eff. collection areas TString collareaName(outPath); collareaName += "area_"; collareaName += XX; collareaName += ".root"; gLog << "collareaName = " << collareaName << endl; //------------------------------ // name of data file to be analysed TString filenameData(outPath); filenameData += typeData; filenameData += ext; gLog << "filenameData = " << filenameData << endl; //------------------------------ // name of output data file (after the final cuts) TString filenameDataout(outPath); filenameDataout += typeData; filenameDataout += "_"; filenameDataout += XX; filenameDataout += extout; gLog << "filenameDataout = " << filenameDataout << endl; //------------------------------ // name of file containing histograms for flux calculastion TString filenameResults(outPath); filenameResults += typeData; filenameResults += "Results_"; filenameResults += XX; filenameResults += extout; gLog << "filenameResults = " << filenameResults << endl; //==================================================================== gLog << "-----------------------------------------------" << endl; gLog << "Start calculation of effective collection areas" << endl; MParList parlist; MTaskList tasklist; //--------------------------------------- // Setup the tasks to be executed // MReadMarsFile reader("Events", filenameArea); reader.DisableAutoScheme(); MFCT1SelFinal cuthadrons; cuthadrons.SetHadronnessName(fhadronnessName); cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); MContinue conthadrons(&cuthadrons); MHMcCT1CollectionArea collarea; collarea.SetEaxis(MHMcCT1CollectionArea::kLinear); MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); filler.SetName("CollectionArea"); //******************************** // entries in MParList parlist.AddToList(&tasklist); InitBinnings(&parlist); parlist.AddToList(&collarea); //******************************** // entries in MTaskList tasklist.AddToList(&reader); tasklist.AddToList(&conthadrons); tasklist.AddToList(&filler); //******************************** //----------------------------------------- // Execute event loop // MEvtLoop magic; magic.SetParList(&parlist); MProgressBar bar; magic.SetProgressBar(&bar); if (!magic.Eventloop()) return; tasklist.PrintStatistics(0, kTRUE); // Calculate effective collection areas // and display the histograms // //MHMcCT1CollectionArea *collarea = // (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea"); collarea.CalcEfficiency(); collarea.DrawClone(); // save binnings for call to CT1EEst MBinning *binsE = (MBinning*)parlist.FindObject("BinningE"); if (!binsE) { gLog << "Object 'BinningE' not found in MParList" << endl; return; } MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta"); if (!binsTheta) { gLog << "Object 'BinningTheta' not found in MParList" << endl; return; } //--------------------------------------------- // Write histograms to a file // TFile f(collareaName, "RECREATE"); collarea.GetHist()->Write(); collarea.GetHAll()->Write(); collarea.GetHSel()->Write(); f.Close(); gLog << "Collection area plots written onto file " << collareaName << endl; gLog << "Calculation of effective collection areas done" << endl; gLog << "-----------------------------------------------" << endl; //------------------------------------------------------------------ TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; if (OEEst) { //=========================================================== // // Optimization of energy estimator // gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl; TString inpath(""); TString outpath(""); Int_t howMany = 2000; CT1EEst(inpath, filenameOpt, outpath, energyParName, fHilName, fHilNameSrc, fhadronnessName, howMany, maxhadronness, maxalpha, maxdist, binsE, binsTheta); gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl; } if (WEX) { //----------------------------------------------------------- // // Read in parameters of energy estimator ("MMcEnergyEst") // and migration matrix ("MHMcEnergyMigration") // gLog << "================================================================" << endl; gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '" << energyParName << "'" << endl; TFile enparam(energyParName); enparam.ls(); MMcEnergyEst mcest("MMcEnergyEst"); mcest.Read("MMcEnergyEst"); //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst")); gLog << "Parameters of energy estimator were read in" << endl; gLog << "Read in Migration matrix" << endl; MHMcEnergyMigration mighiston("MHMcEnergyMigration"); mighiston.Read("MHMcEnergyMigration"); //MHMcEnergyMigration &mighiston = // *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration")); gLog << "Migration matrix was read in" << endl; TArrayD parA(mcest.GetNumCoeffA()); TArrayD parB(mcest.GetNumCoeffB()); for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); //....................................................................... gLog << "CT1Analysis.C : write root file '" << filenameDataout << "'" << endl; //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); MWriteRootFile write(filenameDataout, "RECREATE"); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); write.AddContainer("MMcEvt", "Events"); write.AddContainer("ThetaOrig", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasExt", "Events"); write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); //write.AddContainer("HadNN", "Events"); write.AddContainer("HadSC", "Events"); write.AddContainer("HadRF", "Events"); write.AddContainer("MEnergyEst", "Events"); //----------------------------------------------------------------- // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); MFCT1SelFinal selfinalgh(fHilNameSrc); selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); selfinalgh.SetHadronnessName(fhadronnessName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); //fillhadnn.SetName("HhadNN"); MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); fillhadsc.SetName("HhadSC"); MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); fillhadrf.SetName("HhadRF"); //--------------------------- // calculate estimated energy MEnergyEstParam eeston(fHilName); eeston.Add(fHilNameSrc); eeston.SetCoeffA(parA); eeston.SetCoeffB(parB); //--------------------------- // calculate estimated energy using Daniel's parameters //MEnergyEstParamDanielMkn421 eeston(fHilName); //eeston.Add(fHilNameSrc); //eeston.SetCoeffA(parA); //eeston.SetCoeffB(parB); //--------------------------- MFillH hfill1("MHHillas", fHilName); hfill1.SetName("HHillas"); MFillH hfill2("MHStarMap", fHilName); hfill2.SetName("HStarMap"); MFillH hfill3("MHHillasExt", fHilNameSrc); hfill3.SetName("HHillasExt"); MFillH hfill4("MHHillasSrc", fHilNameSrc); hfill4.SetName("HHillasSrc"); MFillH hfill5("MHNewImagePar", fImgParName); hfill5.SetName("HNewImagePar"); //--------------------------- // new from Robert MFillH hfill6("MHTimeDiffTheta", "MMcEvt"); hfill6.SetName("HTimeDiffTheta"); MFillH hfill6a("MHTimeDiffTime", "MMcEvt"); hfill6a.SetName("HTimeDiffTime"); MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc); hfill7.SetName("HAlphaEnergyTheta"); MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc); hfill7a.SetName("HAlphaEnergyTime"); MFillH hfill7b("MHThetabarTime", fHilNameSrc); hfill7b.SetName("HThetabarTime"); MFillH hfill7c("MHEnergyTime", "MMcEvt"); hfill7c.SetName("HEnergyTime"); //--------------------------- MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(fhadronnessName); selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); contfinal.SetName("ContSelFinal"); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); //***************************** // entries in MTaskList tliston.AddToList(&read); // robert tliston.AddToList(&hfill6); //timediff tliston.AddToList(&hfill6a); //timediff tliston.AddToList(&contfinalgh); tliston.AddToList(&eeston); tliston.AddToList(&write); //tliston.AddToList(&fillhadnn); tliston.AddToList(&fillhadsc); tliston.AddToList(&fillhadrf); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); //robert tliston.AddToList(&hfill7); tliston.AddToList(&hfill7a); tliston.AddToList(&hfill7b); tliston.AddToList(&hfill7c); tliston.AddToList(&contfinal); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); Int_t maxevents = -1; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); //------------------------------------------- // Display the histograms // //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); gLog << "before hadRF" << endl; pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); gLog << "before hadSC" << endl; pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); gLog << "before MHHillas" << endl; pliston.FindObject("MHHillas")->DrawClone(); gLog << "before MHHillasExt" << endl; pliston.FindObject("MHHillasExt")->DrawClone(); gLog << "before MHHillasSrc" << endl; pliston.FindObject("MHHillasSrc")->DrawClone(); gLog << "before MHNewImagePar" << endl; pliston.FindObject("MHNewImagePar")->DrawClone(); gLog << "before MHStarMap" << endl; pliston.FindObject("MHStarMap")->DrawClone(); gLog << "before DeleteBinnings" << endl; DeleteBinnings(&pliston); gLog << "before Robert's code" << endl; //rwagner write all relevant histograms onto a file if (WRobert) { gLog << "=======================================================" << endl; gLog << "Write results onto file '" << filenameResults << "'" << endl; TFile outfile(filenameResults,"recreate"); MHHillasSrc* hillasSrc = (MHHillasSrc*)(pliston->FindObject("MHHillasSrc")); TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); alphaHist->Write(); gLog << "Alpha plot has been written out" << endl; MHAlphaEnergyTheta* aetH = (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta")); TH3D* aetHist = (TH3D*)(aetH->GetHist()); aetHist->SetName("aetHist"); aetHist->Write(); gLog << "AlphaEnergyTheta plot has been written out" << endl; MHAlphaEnergyTime* aetH2 = (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime")); TH3D* aetHist2 = (TH3D*)(aetH2->GetHist()); aetHist2->SetName("aetimeHist"); // aetHist2->DrawClone(); aetHist2->Write(); gLog << "AlphaEnergyTime plot has been written out" << endl; MHThetabarTime* tbt = (MHThetabarTime*)(pliston->FindObject("MHThetabarTime")); TProfile* tbtHist = (TProfile*)(tbt->GetHist()); tbtHist->SetName("tbtHist"); tbtHist->Write(); gLog << "ThetabarTime plot has been written out" << endl; MHEnergyTime* ent = (MHEnergyTime*)(pliston->FindObject("MHEnergyTime")); TH2D* entHist = (TH2D*)(ent->GetHist()); entHist->SetName("entHist"); entHist->Write(); gLog << "EnergyTime plot has been written out" << endl; MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta"); TH2D* timeHist = (TH2D*)(time->GetHist()); timeHist->SetName("MHTimeDiffTheta"); timeHist->SetTitle("Time diffs"); timeHist->Write(); gLog << "TimeDiffTheta plot has been written out" << endl; MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime"); TH2D* timeHist2 = (TH2D*)(time2->GetHist()); timeHist2->SetName("MHTimeDiffTime"); timeHist2->SetTitle("Time diffs"); timeHist2->Write(); gLog << "TimeDiffTime plot has been written out" << endl; //rwagner write also collareas to same file collarea->GetHist()->Write(); collarea->GetHAll()->Write(); collarea->GetHSel()->Write(); gLog << "Effective collection areas have been written out" << endl; //rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data //rwagner (ON/OFF/MC), MJDmin, MJDmax to this file gLog << "before closing outfile" << endl; //outfile.Close(); gLog << "Results were written onto file '" << filenameResults << "'" << endl; gLog << "=======================================================" << endl; } } gLog << "Macro CT1Analysis : End of Job E_XX" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- }