//#include "MagicEgyEst.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 ONOFFAnalysis() { gLog.SetNoColors(); if (gRandom) delete gRandom; gRandom = new TRandom3(0); //----------------------------------------------- //const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; const char *offfile = "/data/MAGIC/rootdata/2003_11_29/20031128_032*_D_OffCrab1_E.root"; //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc"; // const char *onfile = "~magican/ct1test/wittek/mkn421_00-01"; const char *onfile = "/data/MAGIC/rootdata/2003_11_29/20031128_031*_D_Crab-Nebula_E.root"; const char *mcfile = "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root"; //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_4/gh/0/0/G_M0_00_0_550*.root"; //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_5/gh/0/0/G_M0_00_0_550*.root"; //----------------------------------------------- // path for input for Mars TString inPath = "/data/MAGIC/scratch/wittek/"; // path for output from Mars TString outPath = "/data/MAGIC/scratch/wittek/"; //----------------------------------------------- //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 = kTRUE; Bool_t GPad = kFALSE; // generate padding histograms? Bool_t WPad = kFALSE; // write out padding histograms ? Bool_t RPad = kFALSE; // read in padding histograms ? Bool_t Wout = kTRUE; // 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 CTrainRF = TRUE : create matrices of training events // and root files of training and test events // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas // - 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 // and root files of training and test 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 = kTRUE; Bool_t CMatrix = kTRUE; // 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 = kFALSE; Bool_t CCollArea= kFALSE; // calculate eff. collection areas Bool_t OEEst = kFALSE; // optimize energy estimator Bool_t WEX = kFALSE; // update root file ? Bool_t WRobert = kFALSE; // 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 MagicAnalysis : Start of Job A" << endl; gLog << "" << endl; gLog << "Macro MagicAnalysis : 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 MPadONOFF pad; pad.SetName("MPadONOFF"); pad.SetPadFlag(1); pad.SetDataType(typeInput); // generate the padding histograms if (GPad) { gLog << "=====================================================" << endl; gLog << "Start generating the padding histograms" << endl; //----------------------------------------- // ON events gLog << "-----------" << endl; gLog << "ON events :" << endl; gLog << "-----------" << endl; MTaskList tliston; MParList pliston; MReadMarsFile readON("Events", fileON); read.DisableAutoScheme(); MCT1ReadPreProc readON(fileON); //MFSelBasic selthetaon; //selthetaon.SetCuts(-100.0, 29.5, 35.5); //MContinue contthetaon(&selthetaon); MBlindPixelCalc blindon; blindon.SetUseBlindPixels(); MFSelBasic 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; MReadMarsFile readOFF("Events", fileOFF); read.DisableAutoScheme(); // MCT1ReadPreProc readOFF(fileOFF); MFSelBasic selthetaoff; selthetaoff.SetCuts(-100.0, 29.5, 35.5); MContinue contthetaoff(&selthetaoff); MBlindPixelCalc blindoff; blindoff.SetUseBlindPixels(); MFSelBasic 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("MGeomCamMagic", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenamein); read.DisableAutoScheme(); MGeomApply apply; MMcPedestalCopy pcopy; MMcPedestalNSBAdd pnsb; MPedestalWorkaround waround; // a way to find out whether one is dealing with MC : MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5); // MC f1.SetName("Select MC"); MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5); // data f2.SetName("Select Data"); MCerPhotCalc ncalc; ncalc.SetFilter(&f1); MCerPhotAnal2 nanal; nanal.SetFilter(&f2); //if (typeInput == "ON") //{ // MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", // "MCT1PointingCorrCalc"); //} MBlindPixelCalc blindbeforepad; //blindbeforepad.SetUseBlindPixels(); blindbeforepad.SetName("BlindBeforePadding"); MBlindPixelCalc blind; //blind.SetUseBlindPixels(); blind.SetName("BlindAfterPadding"); MFSelBasic 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"); // -------------------------------------------------- MFSelStandard 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("MMcRunHeader", "RunHeaders", kFALSE); //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(&f1); tliston.AddToList(&f2); tliston.AddToList(&apply); tliston.AddToList(&pcopy); //tliston.AddToList(&waround); tliston.AddToList(&pnsb); tliston.AddToList(&ncalc); tliston.AddToList(&nanal); 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 MagicAnalysis : 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 MagicAnalysis : Start of Job B_RF_UP" << endl; gLog << "" << endl; gLog << "Macro MagicAnalysis : 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"; TString extin = "1.root"; TString extout = "2.root"; //-------------------------------------------- // for the analysis using ON data only set typeMatrixHadrons = "ON" // ON and OFF data = "OFF" TString typeMatrixHadrons = "OFF"; gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; // 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 NameData = outPath; NameData += typeInput; TString inNameData(NameData); inNameData += extin; gLog << "inNameData = " << inNameData << endl; // name of output root file TString outNameData(NameData); outNameData += extout; gLog << "outNameData = " << outNameData << endl; //-------------------------------------------- // files to be read for generating // - the matrices of training events // - and the root files of training and test events // "hadrons" : TString filenameHad = outPath; filenameHad += typeMatrixHadrons; filenameHad += extin; Int_t howManyHadronsTrain = 12000; Int_t howManyHadronsTest = 12000; gLog << "filenameHad = " << filenameHad << ", howManyHadronsTrain = " << howManyHadronsTrain << ", howManyHadronsTest = " << howManyHadronsTest << endl; // "gammas" : TString filenameMC = outPath; filenameMC += "MC"; filenameMC += extin; Int_t howManyGammasTrain = 12000; Int_t howManyGammasTest = 12000; gLog << "filenameMC = " << filenameMC << ", howManyGammasTrain = " << howManyGammasTrain << ", howManyGammasTest = " << howManyGammasTest << endl; //-------------------------------------------- // files for the matrices of training events TString NameGammas = outPath; NameGammas += "RFmatrix_gammas_Train_"; NameGammas += "MC"; NameGammas += extin; TString NameHadrons = outPath; NameHadrons += "RFmatrix_hadrons_Train_"; NameHadrons += typeMatrixHadrons; NameHadrons += extin; //-------------------------------------------- // root files for the training events TString NameGammasTrain = outPath; NameGammasTrain += "RF_gammas_Train_"; NameGammasTrain += "MC"; TString inNameGammasTrain(NameGammasTrain); inNameGammasTrain += extin; TString outNameGammasTrain(NameGammasTrain); outNameGammasTrain += extout; TString NameHadronsTrain = outPath; NameHadronsTrain += "RF_hadrons_Train_"; NameHadronsTrain += typeMatrixHadrons; TString inNameHadronsTrain(NameHadronsTrain); inNameHadronsTrain += extin; TString outNameHadronsTrain(NameHadronsTrain); outNameHadronsTrain += extout; //-------------------------------------------- // root files for the test events TString NameGammasTest = outPath; NameGammasTest += "RF_gammas_Test_"; NameGammasTest += "MC"; TString inNameGammasTest(NameGammasTest); inNameGammasTest += extin; TString outNameGammasTest(NameGammasTest); outNameGammasTest += extout; TString NameHadronsTest = outPath; NameHadronsTest += "RF_hadrons_Test_"; NameHadronsTest += typeMatrixHadrons; TString inNameHadronsTest(NameHadronsTest); inNameHadronsTest += extin; TString outNameHadronsTest(NameHadronsTest); outNameHadronsTest += extout; //-------------------------------------------------------------------- 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 = " << NameGammas << endl; gLog << "" << endl; // read in the object with the name 'mtxName' from file 'NameGammas' // TFile fileg(NameGammas); 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 = " << NameHadrons << endl; gLog << "" << endl; // read in the object with the name 'mtxName' from file 'NameHadrons' // TFile fileh(NameHadrons); matrixh.Read(mtxName); matrixh.Print("SizeCols"); } //************************************************************************* // create matrices of training events // and root files of training and test events if (CTrainRF) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << " Create matrices of training events and root files of training and test events" << endl; gLog << " Gammas :" << endl; gLog << "---------" << endl; MParList plistg; MTaskList tlistg; MReadMarsFile readg("Events", filenameMC); readg.DisableAutoScheme(); 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(howManyGammasTrain+howManyGammasTest); selectorg.SetName("selectGammasTrainTest"); selectorg.SetInverted(); //selectorg.SetUseOrigDistribution(kTRUE); MContinue contg(&selectorg); contg.SetName("ContGammas"); Double_t probg = ( (Double_t) howManyGammasTrain ) / ( (Double_t)(howManyGammasTrain+howManyGammasTest) ); MFRandomSplit splitg(probg); MFillH fillmatg("MatrixGammas"); fillmatg.SetFilter(&splitg); fillmatg.SetName("fillGammas"); //----------------------- // for writing the root files of training and test events // for gammas MWriteRootFile writetraing(inNameGammasTrain, "RECREATE"); writetraing.SetName("WriteGammasTrain"); writetraing.SetFilter(&splitg); writetraing.AddContainer("MRawRunHeader", "RunHeaders"); writetraing.AddContainer("MTime", "Events"); writetraing.AddContainer("MMcEvt", "Events"); writetraing.AddContainer("ThetaOrig", "Events"); writetraing.AddContainer("MSrcPosCam", "Events"); writetraing.AddContainer("MSigmabar", "Events"); writetraing.AddContainer("MHillas", "Events"); writetraing.AddContainer("MHillasExt", "Events"); writetraing.AddContainer("MHillasSrc", "Events"); writetraing.AddContainer("MNewImagePar", "Events"); MContinue contgtrain(&splitg); contgtrain.SetName("ContGammaTrain"); MWriteRootFile writetestg(inNameGammasTest, "RECREATE"); writetestg.SetName("WriteGammasTest"); writetestg.AddContainer("MRawRunHeader", "RunHeaders"); writetestg.AddContainer("MTime", "Events"); writetestg.AddContainer("MMcEvt", "Events"); writetestg.AddContainer("ThetaOrig", "Events"); writetestg.AddContainer("MSrcPosCam", "Events"); writetestg.AddContainer("MSigmabar", "Events"); writetestg.AddContainer("MHillas", "Events"); writetestg.AddContainer("MHillasExt", "Events"); writetestg.AddContainer("MHillasSrc", "Events"); writetestg.AddContainer("MNewImagePar", "Events"); //----------------------- //***************************** fill gammas *** // entries in MParList plistg.AddToList(&tlistg); InitBinnings(&plistg); plistg.AddToList(&matrixg); //***************************** // entries in MTaskList tlistg.AddToList(&readg); tlistg.AddToList(&contg); tlistg.AddToList(&splitg); tlistg.AddToList(&fillmatg); tlistg.AddToList(&writetraing); tlistg.AddToList(&contgtrain); tlistg.AddToList(&writetestg); //***************************** 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); matrixg.Print("SizeCols"); Int_t generatedgTrain = matrixg.GetM().GetNrows(); if ( fabs(generatedgTrain-howManyGammasTrain) > 3.0*sqrt(howManyGammasTrain) ) { gLog << "ONOFFAnalysis.C : no.of generated gamma training events (" << generatedgTrain << ") is incompatible with the no.of requested events (" << howManyGammasTrain << ")" << endl; } Int_t generatedgTest = writetestg.GetNumExecutions(); if ( fabs(generatedgTest-howManyGammasTest) > 3.0*sqrt(howManyGammasTest) ) { gLog << "ONOFFAnalysis.C : no.of generated gamma test events (" << generatedgTest << ") is incompatible with the no.of requested events (" << howManyGammasTest << ")" << endl; } //***************************** fill hadrons *** gLog << "---------------------------------------------------------------" << endl; gLog << " Hadrons :" << endl; gLog << "----------" << endl; MParList plisth; MTaskList tlisth; MReadMarsFile readh("Events", filenameHad); readh.DisableAutoScheme(); 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); //use the original distribution from the gammas MH3 &href = *(selectorg.GetHistOrig()); MFEventSelector2 selectorh(href); selectorh.SetNumMax(howManyHadronsTrain+howManyHadronsTest); selectorh.SetName("selectHadronsTrainTest"); selectorh.SetInverted(); MContinue conth(&selectorh); conth.SetName("ContHadrons"); Double_t probh = ( (Double_t) howManyHadronsTrain ) / ( (Double_t)(howManyHadronsTrain+howManyHadronsTest) ); MFRandomSplit splith(probh); MFillH fillmath("MatrixHadrons"); fillmath.SetFilter(&splith); fillmath.SetName("fillHadrons"); //----------------------- // for writing the root files of training and test events // for hadrons MWriteRootFile writetrainh(inNameHadronsTrain, "RECREATE"); writetrainh.SetName("WriteHadronsTrain"); writetrainh.SetFilter(&splith); writetrainh.AddContainer("MRawRunHeader", "RunHeaders"); writetrainh.AddContainer("MTime", "Events"); writetrainh.AddContainer("MMcEvt", "Events"); writetrainh.AddContainer("ThetaOrig", "Events"); writetrainh.AddContainer("MSrcPosCam", "Events"); writetrainh.AddContainer("MSigmabar", "Events"); writetrainh.AddContainer("MHillas", "Events"); writetrainh.AddContainer("MHillasExt", "Events"); writetrainh.AddContainer("MHillasSrc", "Events"); writetrainh.AddContainer("MNewImagePar", "Events"); MContinue conthtrain(&splith); MWriteRootFile writetesth(inNameHadronsTest, "RECREATE"); writetesth.SetName("WriteHadronsTest"); writetesth.AddContainer("MRawRunHeader", "RunHeaders"); writetesth.AddContainer("MTime", "Events"); writetesth.AddContainer("MMcEvt", "Events"); writetesth.AddContainer("ThetaOrig", "Events"); writetesth.AddContainer("MSrcPosCam", "Events"); writetesth.AddContainer("MSigmabar", "Events"); writetesth.AddContainer("MHillas", "Events"); writetesth.AddContainer("MHillasExt", "Events"); writetesth.AddContainer("MHillasSrc", "Events"); writetesth.AddContainer("MNewImagePar", "Events"); //***************************** // entries in MParList plisth.AddToList(&tlisth); InitBinnings(&plisth); plisth.AddToList(&matrixh); //***************************** // entries in MTaskList tlisth.AddToList(&readh); tlisth.AddToList(&conth); tlisth.AddToList(&splith); tlisth.AddToList(&fillmath); tlisth.AddToList(&writetrainh); tlisth.AddToList(&conthtrain); tlisth.AddToList(&writetesth); //***************************** 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); matrixh.Print("SizeCols"); Int_t generatedhTrain = matrixh.GetM().GetNrows(); if ( fabs(generatedhTrain-howManyHadronsTrain) > 3.0*sqrt(howManyHadronsTrain) ) { gLog << "ONOFFAnalysis.C : no.of generated hadron training events (" << generatedhTrain << ") is incompatible with the no.of requested events (" << howManyHadronsTrain << ")" << endl; } Int_t generatedhTest = writetesth.GetNumExecutions(); if ( fabs(generatedhTest-howManyHadronsTest) > 3.0*sqrt(howManyHadronsTest) ) { gLog << "ONOFFAnalysis.C : no.of generated gamma test events (" << generatedhTest << ") is incompatible with the no.of requested events (" << howManyHadronsTest << ")" << endl; } //***************************************************** // 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(NameGammas, "RECREATE", ""); matrixg.Write(); gLog << "" << endl; gLog << "Macro MagicAnalysis : matrix of training events for gammas written onto file " << NameGammas << endl; //------------------------------------------- // "hadrons" gLog << "Hadrons :" << endl; matrixh.Print("SizeCols"); TFile writeh(NameHadrons, "RECREATE", ""); matrixh.Write(); gLog << "" << endl; gLog << "Macro MagicAnalysis : matrix of training events for hadrons written onto file " << NameHadrons << 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); gLog << "ONOFFAnalysis : RF trees were read in from file " << outRF << endl; // 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 MagicAnalysis : 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 root files with the RF hadronness // if (WRF) { //TString fileName(inNameHadronsTrain); //TString outName(outNameHadronsTrain); //TString fileName(inNameHadronsTest); //TString outName(outNameHadronsTest); //TString fileName(inNameGammasTrain); //TString outName(outNameGammasTrain); //TString fileName(inNameGammasTest); //TString outName(outNameGammasTest); TString fileName(inNameData); TString outName(outNameData); gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Update root file '" << fileName << "' with the RF hadronness; ==> " << outName << endl; MTaskList tliston; MParList pliston; // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", fileName); read.DisableAutoScheme(); //....................................................................... // calculate hadronnes for method of RANDOM FOREST MRanForestCalc rfcalc; rfcalc.SetHadronnessName(hadRFName); //....................................................................... //MWriteRootFile write(outName, "UPDATE"); MWriteRootFile write(outName, "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"); //----------------------------------------------------------------- MFSelFinal 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"); MFSelFinal 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("hadRF", "MHHadronness")->Print(); 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 '" << fileName << "' : " << 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 MagicAnalysis : 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 MagicAnalysis : Start of Job B_SC_UP" << endl; gLog << "" << endl; gLog << "Macro MagicAnalysis : 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 // MSupercuts 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 MFindSupercuts findsuper; findsuper.SetFilenameParam(parSCfile); findsuper.SetHadronnessName("HadSC"); //findsuper.SetUseOrigDistribution(kTRUE); //-------------------------- // create matrices and write them onto files if (CMatrix) { TString mname("costheta"); MBinning bin("Binning"+mname); bin.SetEdges(10, 0., 1.0); MH3 mh3("cos(MMcEvt.fTelescopeTheta)"); mh3.SetName(mname); MH::SetBinning(&mh3.GetHist(), &bin); //for (Int_t i=1; i<=mh3.GetNbins(); i++) // mh3.GetHist().SetBinContent(i, 1.0); if (filenameTrain == filenameTest) { if ( !findsuper.DefineTrainTestMatrix( filenameTrain, mh3, howManyTrain, howManyTest, fileMatrixTrain, fileMatrixTest) ) { *fLog << "MagicAnalysis.C : DefineTrainTestMatrix failed" << endl; return; } } else { if ( !findsuper.DefineTrainMatrix(filenameTrain, mh3, howManyTrain, fileMatrixTrain) ) { *fLog << "MagicAnalysis.C : DefineTrainMatrix failed" << endl; return; } if ( !findsuper.DefineTestMatrix( filenameTest, mh3, howManyTest, fileMatrixTest) ) { *fLog << "MagicAnalysis.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 MSupercuts constructor if (WOptimize) { gLog << "MagicAnalysis.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 << "MagicAnalysis.C : optimization of supercuts failed" << endl; return; } } //-------------------------------------- // test the supercuts on the test sample // if (RTest) { gLog << "MagicAnalysis.C : test the supercuts on the test matrix" << endl; Bool_t rt = findsuper.TestParams(); if (!rt) { gLog << "MagicAnalysis.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); MSupercuts scin; scin.Read("MSupercuts"); 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("MGeomCamMagic", "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"; MSupercutsCalc 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("MGeomCamMagic", "MGeomCam"); Float_t maxhadronness = 0.40; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; MFSelFinal 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"); MFSelFinal 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 MagicAnalysis : 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 MagicAnalysis : Start of Job C" << endl; gLog << "" << endl; gLog << "Macro MagicAnalysis : JobC = " << (JobC ? "kTRUE" : "kFALSE") << endl; TString ext("2.root"); TString extout("2.root"); TString typeHadrons("OFF"); TString typeGammas("MC"); //-------------------------------------------- // name of input data file TString NameData = outPath; NameData += typeHadrons; TString inNameData(NameData); inNameData += ext; gLog << "inNameData = " << inNameData << endl; // name of input MC file TString NameMC = outPath; NameMC += typeGammas; TString inNameMC(NameMC); inNameMC += ext; gLog << "inNameMC = " << inNameMC << endl; //-------------------------------------------- // root files for the training events TString NameGammasTrain = outPath; NameGammasTrain += "RF_gammas_Train_"; NameGammasTrain += typeGammas; TString outNameGammasTrain(NameGammasTrain); outNameGammasTrain += extout; TString NameHadronsTrain = outPath; NameHadronsTrain += "RF_hadrons_Train_"; NameHadronsTrain += typeHadrons; TString outNameHadronsTrain(NameHadronsTrain); outNameHadronsTrain += extout; //-------------------------------------------- // root files for the test events TString NameGammasTest = outPath; NameGammasTest += "RF_gammas_Test_"; NameGammasTest += typeGammas; TString outNameGammasTest(NameGammasTest); outNameGammasTest += extout; TString NameHadronsTest = outPath; NameHadronsTest += "RF_hadrons_Test_"; NameHadronsTest += typeHadrons; TString outNameHadronsTest(NameHadronsTest); outNameHadronsTest += extout; //-------------------------------------------------------------------- //TString filenameData(inNameData); //TString filenameMC(inNameMC); //TString filenameData(outNameHadronsTrain); //TString filenameMC(outNameGammasTrain); TString filenameData(outNameHadronsTest); TString filenameMC(outNameGammasTest); gLog << "filenameData = " << filenameData << endl; gLog << "filenameMC = " << filenameMC << endl; //----------------------------------------------------------------- MTaskList tliston; MParList pliston; // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "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; MFSelFinal 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"); MFSelFinal 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 MagicAnalysis : 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 MagicAnalysis : Start of Job D" << endl; gLog << "" << endl; gLog << "Macro MagicAnalysis : 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("MGeomCamMagic", "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("MGeomCamMagic", "MGeomCam"); MFSelFinal 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"); MFSelFinal 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 MagicAnalysis : 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 MagicAnalysis : Start of Job E_XX" << endl; gLog << "" << endl; gLog << "Macro MagicAnalysis : JobE_XX, CCollArea, OEEst, WEX = " << (JobE_XX ? "kTRUE" : "kFALSE") << ", " << (CCollArea?"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; //==================================================================== MHMcCollectionArea collarea; collarea.SetEaxis(MHMcCollectionArea::kLinear); MParList parlist; InitBinnings(&parlist); if (CCollArea) { gLog << "-----------------------------------------------" << endl; gLog << "Start calculation of effective collection areas" << endl; MTaskList tasklist; //--------------------------------------- // Setup the tasks to be executed // MReadMarsFile reader("Events", filenameArea); reader.DisableAutoScheme(); MFSelFinal cuthadrons; cuthadrons.SetHadronnessName(fhadronnessName); cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); MContinue conthadrons(&cuthadrons); MFillH filler("MHMcCollectionArea", "MMcEvt"); filler.SetName("CollectionArea"); //******************************** // entries in MParList parlist.AddToList(&tasklist); 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 // //MHMcCollectionArea *collarea = // (MHMcCollectionArea*)parlist.FindObject("MHMcCollectionArea"); collarea.CalcEfficiency(); collarea.DrawClone(); //--------------------------------------------- // Write histograms to a file // TFile f(collareaName, "RECREATE"); //collarea.GetHist()->Write(); //collarea.GetHAll()->Write(); //collarea.GetHSel()->Write(); collarea.Write(); f.Close(); gLog << "Collection area plots written onto file " << collareaName << endl; gLog << "Calculation of effective collection areas done" << endl; gLog << "-----------------------------------------------" << endl; //------------------------------------------------------------------ } if (!CCollArea) { gLog << "-----------------------------------------------" << endl; gLog << "Read in effective collection areas from file " << collareaName << endl; TFile collfile(collareaName); collfile.ls(); collarea.Read("MHMcCollectionArea"); collarea.DrawClone(); gLog << "Effective collection areas were read in from file " << collareaName << endl; gLog << "-----------------------------------------------" << endl; } // save binnings for call to MagicEEst 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; } //------------------------------------- TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; if (OEEst) { //=========================================================== // // Optimization of energy estimator // gLog << "Macro MagicAnalysis.C : calling MagicEEst" << endl; TString inpath(""); TString outpath(""); Int_t howMany = 2000; MagicEEst(inpath, filenameOpt, outpath, energyParName, fHilName, fHilNameSrc, fhadronnessName, howMany, maxhadronness, maxalpha, maxdist, binsE, binsTheta); gLog << "Macro MagicAnalysis.C : returning from MagicEEst" << endl; } if (WEX) { //----------------------------------------------------------- // // Read in parameters of energy estimator ("MMcEnergyEst") // and migration matrix ("MHMcEnergyMigration") // gLog << "================================================================" << endl; gLog << "Macro MagicAnalysis.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("MGeomCamMagic", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); //....................................................................... gLog << "MagicAnalysis.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("MGeomCamMagic", "MGeomCam"); MFSelFinal 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"); //--------------------------- MFSelFinal 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 Analysis : End of Job E_XX" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- }