#include "CT1EgyEst.C" void InitBinnings(MParList *plist) { gLog << "InitBinnings" << endl; MBinning *binse = new MBinning("BinningE"); //binse->SetEdgesLog(30, 1.0e2, 1.0e5); binse->SetEdges(30, 2, 5); 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"); Double_t yedge[9] = {7.5, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; TArrayD yed(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(9,zedge); bincosth->SetEdges(zed); plist->AddToList(bincosth); MBinning *binsdiff = new MBinning("BinningDiffsigma2"); binsdiff->SetEdges(100, -5.0, 20.0); plist->AddToList(binsdiff); } void DeleteBinnings(MParList *plist) { gLog << "DeleteBinnings" << endl; TObject *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; } //************************************************************************ void CT1Analysis() { 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/optionB/"; // path for output from Mars TString outPath = "~magican/ct1test/wittek/marsoutput/optionB/"; //----------------------------------------------- TEnv env("macros/CT1env.rc"); Bool_t printEnv = kFALSE; //************************************************************************ // Job A_ON : read ON data // - generate sigmabar vs. Theta plot; // - write root file for ON data (ON1.root); Bool_t JobA_ON = kFALSE; Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ? Bool_t WON1 = kFALSE; // write out root file ON1.root ? // Job A_MC : read MC gamma data, // - read sigmabar vs. Theta plot from ON data // - do padding; // - write root file for MC gammas (MC1.root); Bool_t JobA_MC = kFALSE; Bool_t WMC1 = kFALSE; // write out root file MC1.root ? // Job B_NN_UP : read ON1.root (or MC1.root) file // - depending on RlookNN : create (or read in) hadron and gamma matrix // - calculate hadroness for method of NEAREST NEIGHBORS // and for the SUPERCUTS // - update the input files with the hadronesses (ON1.root or MC1.root) Bool_t JobB_NN_UP = kFALSE; Bool_t RLookNN = kFALSE; // read in look-alike events Bool_t WNN = kFALSE; // update input root file ? // Job B_RF_UP : read ON1.root (or MC1.root) file // - depending on RLook : create (or read in) hadron and gamma matrix // - depending on RTree : create (or read in) trees // - calculate hadroness for method of RANDOM FOREST // - update the input files with the hadroness (ON1.root or MC1.root) Bool_t JobB_RF_UP = kFALSE; Bool_t RLookRF = kFALSE; // read in look-alike events Bool_t RTree = kFALSE; // read in trees Bool_t WRF = kFALSE; // update input root file ? // Job C: // - read ON1.root and MC1.root files // which should have been updated to contain the hadronnesses // for the method of // NEAREST NEIGHBORS // SUPERCUTS and // RF // - produce Neyman-Pearson plots Bool_t JobC = kFALSE; // Job D : // - select g/h separation method XX // - read ON2 (or MC2) root file // - apply cuts in hadronness // - make plots Bool_t JobD = kFALSE; // Job E_EST_UP : // - read MC1.root file // - select g/h separation method XX // - optimize energy estimation for events passing the final cuts // - write parameters of energy estimator onto file // - update ON1.root and MC1.root files with estimated energy // (ON_XX1.root and MC_XX1.root) Bool_t JobE_EST_UP = kFALSE; Bool_t WESTUP = kFALSE; // update root files ? // Job F_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 (ON3.root)) Bool_t JobF_XX = kFALSE; Bool_t WXX = kFALSE; // write out root file ON3.root ? //************************************************************************ //--------------------------------------------------------------------- // Job A_ON //========= // 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_ON) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job A_ON" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobA_ON, WHistON, WON1 = " << JobA_ON << ", " << WHistON << ", " << WON1 << endl; // name of input root file TString filenamein(onfile); // name of output root file TString outNameImage = outPath; outNameImage += "ON"; outNameImage += "1.root"; //-------------------------------------------------- // use for padding sigmabar vs. Theta from ON data TString typeHist = "ON"; gLog << "typeHist = " << typeHist << endl; // name of file to conatin the histograms for the padding TString outNameSigTh = outPath; outNameSigTh += "SigmaTheta_"; outNameSigTh += typeHist; outNameSigTh += ".root"; //----------------------------------------------------------- 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); MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", "MCT1PointingCorrCalc"); MBlindPixelCalc blind; blind.SetUseBlindPixels(); 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"); //$$$$$$$$$$$$$$$$$$$$$$$$$$$$ //MF daniel( "(MRawRunHeader.fRunNumber<13167)||(MRawRunHeader.fRunNumber>13167)" ); //MContinue contdaniel(&daniel); //$$$$$$$$$$$$$$$$$$$$$$$$$$$$ //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&source); //***************************** // entries in MTaskList tliston.AddToList(&read); //$$$$$$$$$$$$$$$$ //tliston.AddToList(&contdaniel); //$$$$$$$$$$$$$$$$ 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); if (WON1) tliston.AddToList(&write); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&pliston); evtloop.ReadEnv(env, "", printEnv); evtloop.SetProgressBar(&bar); //if (WON1) // 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(); //------------------------------------------- // Write histograms onto a file if (WHistON) { MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); MHBlindPixels *blindpixels = (MHBlindPixels*)pliston.FindObject("BlindPixels"); if (!sigtheta || !blindpixels) { gLog << "Object 'SigmaTheta' or 'BlindPixels' not found" << endl; return; } TH2D *fHSigTh = sigtheta->GetSigmaTheta(); TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta(); TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); TH2D *fHBlindId = blindpixels->GetBlindId(); TH2D *fHBlindN = blindpixels->GetBlindN(); TFile outfile(outNameSigTh, "RECREATE"); fHSigTh->Write(); fHSigPixTh->Write(); fHDifPixTh->Write(); fHBlindId->Write(); fHBlindN->Write(); gLog << "" << endl; gLog << "File " << outNameSigTh << " was written out" << endl; } DeleteBinnings(&pliston); gLog << "Macro CT1Analysis : End of Job A_ON" << endl; gLog << "===================================================" << endl; } //--------------------------------------------------------------------- // Job A_MC //========= // read MC gamma data // - to pad them // (using the 2D-histogram "sigmabar versus Theta" // (SigmaTheta_ON.root) of the ON data) // - to write a file of padded MC gamma events (MC1.root) // (after the standard cuts, before the g/h separation) // (to be used together with the corresponding hadron file // for the optimization of the g/h separation) if (JobA_MC) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job A_MC" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobA_MC, WMC1 = " << JobA_MC << ", " << WMC1 << endl; // name of input root file TString filenamein(mcfile); // name of output root file TString outNameImage = outPath; outNameImage += "MC"; outNameImage += "1.root"; //------------------------------------------------ // use for padding sigmabar vs. Theta from ON data TString typeHist = "ON"; gLog << "typeHist = " << typeHist << endl; // name of file containing the histograms for the padding TString outNameSigTh = outPath; outNameSigTh += "SigmaTheta_"; outNameSigTh += typeHist; outNameSigTh += ".root"; //------------------------------------ // Get the histograms "2D-ThetaSigmabar" // and "3D-ThetaPixSigma" // and "3D-ThetaPixDiff" // and "2D-IdBlindPixels" // and "2D-NBlindPixels" gLog << "Reading in file " << outNameSigTh << endl; TFile *infile = new TFile(outNameSigTh); infile->ls(); TH2D *fHSigmaTheta = (TH2D*) gROOT->FindObject("2D-ThetaSigmabar"); if (!fHSigmaTheta) { gLog << "Object '2D-ThetaSigmabar' not found on root file" << endl; return; } gLog << "Object '2D-ThetaSigmabar' was read in" << endl; TH3D *fHSigmaPixTheta = (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); if (!fHSigmaPixTheta) { gLog << "Object '3D-ThetaPixSigma' not found on root file" << endl; return; } gLog << "Object '3D-ThetaPixSigma' was read in" << endl; TH3D *fHDiffPixTheta = (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); if (!fHDiffPixTheta) { gLog << "Object '3D-ThetaPixDiff' not found on root file" << endl; return; } gLog << "Object '3D-ThetaPixDiff' was read in" << endl; TH2D *fHIdBlindPixels = (TH2D*) gROOT->FindObject("2D-IdBlindPixels"); if (!fHIdBlindPixels) { gLog << "Object '2D-IdBlindPixels' not found on root file" << endl; return; } gLog << "Object '2D-IdBlindPixels' was read in" << endl; TH2D *fHNBlindPixels = (TH2D*) gROOT->FindObject("2D-NBlindPixels"); if (!fHNBlindPixels) { gLog << "Object '2D-NBlindPixels' not found on root file" << endl; return; } gLog << "Object '2D-NBlindPixels' was read in" << endl; //------------------------------------ MTaskList tlist; MParList plist; char *sourceName = "MSrcPosCam"; MSrcPosCam source(sourceName); // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)plist->FindCreateObj("MGeomCamCT1", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MCT1ReadPreProc read(filenamein); MBlindPixelCalc blind; blind.SetUseBlindPixels(); MFCT1SelBasic selbasic; MContinue contbasic(&selbasic); contbasic.SetName("SelBasic"); // There are 2 options for Thomas Schweizer's padding // fPadFlag = 1 get Sigmabar from fHSigmaTheta // and Sigma from fHDiffPixTheta // fPadFlag = 2 get Sigma from fHSigmaPixTheta MCT1PadSchweizer padthomas("MCT1PadSchweizer","Task for the padding (Schweizer)"); padthomas.SetHistograms(fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta, fHIdBlindPixels, fHNBlindPixels); padthomas.SetPadFlag(1); MFillH fillblind("MCBlindPixels[MHBlindPixels]", "MBlindPixels"); fillblind.SetName("HBlind"); //........................................... MSigmabarCalc sigbarcalc; MFillH fillsigtheta ("MCSigmaTheta[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 plist.AddToList(&tlist); InitBinnings(&plist); plist.AddToList(&source); //***************************** // entries in MTaskList tlist.AddToList(&read); tlist.AddToList(&padthomas); tlist.AddToList(&blind); tlist.AddToList(&contbasic); tlist.AddToList(&fillblind); tlist.AddToList(&sigbarcalc); tlist.AddToList(&fillsigtheta); tlist.AddToList(&clean); tlist.AddToList(&hcalc); tlist.AddToList(&hsrccalc); tlist.AddToList(&hfill1); tlist.AddToList(&hfill2); tlist.AddToList(&hfill3); tlist.AddToList(&hfill4); tlist.AddToList(&hfill5); tlist.AddToList(&contstandard); if (WMC1) tlist.AddToList(&write); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.ReadEnv(env, "", printEnv); evtloop.SetProgressBar(&bar); //if (WMC1) // evtloop.Write(); Int_t maxevents = -1; //Int_t maxevents = 1000; if ( !evtloop.Eventloop(maxevents) ) return; tlist.PrintStatistics(0, kTRUE); //------------------------------------------- // Display the histograms // plist.FindObject("MCSigmaTheta", "MHSigmaTheta")->DrawClone(); plist.FindObject("MCBlindPixels", "MHBlindPixels")->DrawClone(); plist.FindObject("MHHillas")->DrawClone(); plist.FindObject("MHHillasExt")->DrawClone(); plist.FindObject("MHHillasSrc")->DrawClone(); plist.FindObject("MHNewImagePar")->DrawClone(); plist.FindObject("MHStarMap")->DrawClone(); DeleteBinnings(&plist); gLog << "Macro CT1Analysis : End of Job A_MC" << endl; gLog << "=========================================================" << endl; } //--------------------------------------------------------------------- // Job B_NN_UP //============ // - read in the matrices of look-alike events for gammas and hadrons // then read ON1.root (or MC1.root) file // // - calculate the hadroness for the method of nearest neighbors // // - update input root file, including the hadroness if (JobB_NN_UP) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job B_NN_UP" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobB_NN_UP, RLookNN, WNN = " << JobB_NN_UP << ", " << RLookNN << ", " << WNN << endl; //-------------------------------------------- // file to be updated (either ON or MC) TString typeInput = "ON"; //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 look-alike events // "hadrons" : TString filenameON = outPath; filenameON += "ON"; filenameON += "1.root"; Int_t howManyHadrons = 500000; gLog << "filenameON = " << filenameON << ", howManyHadrons = " << howManyHadrons << endl; // "gammas" : TString filenameMC = outPath; filenameMC += "MC"; filenameMC += "1.root"; Int_t howManyGammas = 50000; gLog << "filenameMC = " << filenameMC << ", howManyGammas = " << howManyGammas << endl; //-------------------------------------------- // files of look-alike events TString outNameGammas = outPath; outNameGammas += "matrix_gammas_"; outNameGammas += "MC"; outNameGammas += ".root"; TString typeMatrixHadrons = "ON"; gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; TString outNameHadrons = outPath; outNameHadrons += "matrix_hadrons_"; outNameHadrons += typeMatrixHadrons; outNameHadrons += ".root"; MHMatrix matrixg("MatrixGammas"); MHMatrix matrixh("MatrixHadrons"); //************************************************************************* // read in matrices of look-alike events if (RLookNN) { 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("cols"); //***************************************************************** 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("cols"); } //************************************************************************* // create matrices of look-alike events if (!RLookNN) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << " Create matrices of look-alike 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", filenameON); readh.DisableAutoScheme(); MFParticleId fgamma("MMcEvt", '=', kGAMMA); fgamma.SetName("gammaID"); MFParticleId fhadrons("MMcEvt", '!', kGAMMA); fhadrons.SetName("hadronID)"); MFEventSelector selectorg; selectorg.SetNumSelectEvts(howManyGammas); selectorg.SetName("selectGammas"); MFEventSelector selectorh; selectorh.SetNumSelectEvts(howManyHadrons); selectorh.SetName("selectHadrons"); MFillH fillmatg("MatrixGammas"); fillmatg.SetFilter(&flistg); fillmatg.SetName("fillGammas"); MFillH fillmath("MatrixHadrons"); fillmath.SetFilter(&flisth); fillmath.SetName("fillHadrons"); //***************************** 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.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; // 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.SetParList(&plisth); evtlooph.ReadEnv(env, "", printEnv); evtlooph.SetProgressBar(&matrixbar); Int_t maxevents = -1; if (!evtlooph.Eventloop(maxevents)) return; tlisth.PrintStatistics(0, kTRUE); //***************************************************** //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ // // ---------------------------------------------------------- // Definition of the reference sample of look-alike events // (this is a subsample of the original sample) // ---------------------------------------------------------- // gLog << "" << endl; gLog << "========================================================" << endl; // Select a maximum of nmaxevts events from the sample of look-alike // events. They will form the reference sample. Int_t nmaxevts = 2000; // target distribution for the variable in column refcolumn (start at 0); //Int_t refcolumn = 0; //Int_t nbins = 5; //Float_t frombin = 0.5; //Float_t tobin = 1.0; //TH1F *thsh = new TH1F("thsh","target distribution", // nbins, frombin, tobin); //thsh->SetDirectory(NULL); //thsh->SetXTitle("cos( \\Theta)"); //thsh->SetYTitle("Counts"); //Float_t dbin = (tobin-frombin)/nbins; //Float_t lbin = frombin +dbin*0.5; //for (Int_t j=1; j<=nbins; j++) //{ // thsh->Fill(lbin,1.0); // lbin += dbin; //} Int_t refcolumn = 0; MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); TH1F *thsh = new TH1F(); thsh->SetNameTitle("thsh","target distribution"); MH::SetBinning(thsh, binscostheta); Int_t nbins = thsh->GetNbinsX(); Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; for (Int_t j=1; j<=nbins; j++) { thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); } gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " << refcolumn << ", " << nmaxevts << endl; matrixg.EnableGraphicalOutput(); Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); if ( !matrixok ) { gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; return; } gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " << refcolumn << ", " << nmaxevts << endl; matrixh.EnableGraphicalOutput(); matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); delete thsh; if ( !matrixok ) { gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; return; } //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ // write out look-alike events gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Write out look=alike events" << endl; //------------------------------------------- // "gammas" gLog << "Gammas :" << endl; matrixg.Print("cols"); TFile writeg(outNameGammas, "RECREATE", ""); matrixg.Write(); gLog << "" << endl; gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " << outNameGammas << endl; //------------------------------------------- // "hadrons" gLog << "Hadrons :" << endl; matrixh.Print("cols"); TFile writeh(outNameHadrons, "RECREATE", ""); matrixh.Write(); gLog << "" << endl; gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " << outNameHadrons << endl; } //********** end of creating matrices of look-alike events *********** //----------------------------------------------------------------- // Update the input files with the NN and SC hadronness // if (WNN) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Update input file '" << filenameData << "' with the NN and SC 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(); TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; //....................................................................... // calculation of hadroness for method of Nearest Neighbors TString hadNNName = "HadNN"; MMultiDimDistCalc nncalc; nncalc.SetUseNumRows(25); nncalc.SetUseKernelMethod(kFALSE); nncalc.SetHadronnessName(hadNNName); //....................................................................... // 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(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("HadNN", "Events"); write.AddContainer("HadSC", "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(hadNNName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); fillhadnn.SetName("HhadNN"); MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); fillhadsc.SetName("HhadSC"); MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(hadNNName); 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); pliston.AddToList(&matrixg); pliston.AddToList(&matrixh); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&nncalc); tliston.AddToList(&sccalc); tliston.AddToList(&fillhadnn); tliston.AddToList(&fillhadsc); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); tliston.AddToList(&write); tliston.AddToList(&contfinalgh); 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(); 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(); DeleteBinnings(&pliston); } gLog << "Macro CT1Analysis : End of Job B_NN_UP" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- // Job B_RF_UP //============ // - create (or read in) the matrices of look-alike 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, RLookRF, RTree, WRF = " << JobB_RF_UP << ", " << RLookRF << ", " << RTree << ", " << WRF << endl; //-------------------------------------------- // parameters for the random forest Int_t NumTrees = 100; Int_t NumTry = 3; Int_t NdSize = 3; //-------------------------------------------- // file to be updated (either ON or MC) TString typeInput = "ON"; //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 generating the look-alike events // "hadrons" : TString filenameON = outPath; filenameON += "ON"; filenameON += "1.root"; Int_t howManyHadrons = 1000000; gLog << "filenameON = " << filenameON << ", howManyHadrons = " << howManyHadrons << endl; // "gammas" : TString filenameMC = outPath; filenameMC += "MC"; filenameMC += "1.root"; Int_t howManyGammas = 50000; gLog << "filenameMC = " << filenameMC << ", howManyGammas = " << howManyGammas << endl; //-------------------------------------------- // files of look-alike events TString outNameGammas = outPath; outNameGammas += "RFmatrix_gammas_"; outNameGammas += "MC"; outNameGammas += ".root"; TString typeMatrixHadrons = "ON"; gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; TString outNameHadrons = outPath; outNameHadrons += "RFmatrix_hadrons_"; outNameHadrons += typeMatrixHadrons; outNameHadrons += ".root"; MHMatrix matrixg("MatrixGammas"); MHMatrix matrixh("MatrixHadrons"); //-------------------------------------------- // file of trees of the random forest TString outRF = outPath; outRF += "RF.root"; //************************************************************************* // read in matrices of look-alike events if (RLookRF) { 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("cols"); //***************************************************************** 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("cols"); } //************************************************************************* // create matrices of look-alike events if (!RLookRF) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << " Create matrices of look-alike 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", filenameON); readh.DisableAutoScheme(); MFParticleId fgamma("MMcEvt", '=', kGAMMA); fgamma.SetName("gammaID"); MFParticleId fhadrons("MMcEvt", '!', kGAMMA); fhadrons.SetName("hadronID)"); MFEventSelector selectorg; selectorg.SetNumSelectEvts(howManyGammas); selectorg.SetName("selectGammas"); MFEventSelector selectorh; selectorh.SetNumSelectEvts(howManyHadrons); selectorh.SetName("selectHadrons"); MFillH fillmatg("MatrixGammas"); fillmatg.SetFilter(&flistg); fillmatg.SetName("fillGammas"); MFillH fillmath("MatrixHadrons"); fillmath.SetFilter(&flisth); fillmath.SetName("fillHadrons"); //***************************** 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.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; // 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.SetParList(&plisth); evtlooph.ReadEnv(env, "", printEnv); evtlooph.SetProgressBar(&matrixbar); Int_t maxevents = -1; if (!evtlooph.Eventloop(maxevents)) return; tlisth.PrintStatistics(0, kTRUE); //***************************************************** //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ // // ---------------------------------------------------------- // Definition of the reference sample of look-alike events // (this is a subsample of the original sample) // ---------------------------------------------------------- // gLog << "" << endl; gLog << "========================================================" << endl; // Select a maximum of nmaxevts events from the sample of look-alike // events. They will form the reference sample. Int_t nmaxevts = 10000; // target distribution for the variable in column refcolumn (start at 0); //Int_t refcolumn = 0; //Int_t nbins = 5; //Float_t frombin = 0.5; //Float_t tobin = 1.0; //TH1F *thsh = new TH1F("thsh","target distribution", // nbins, frombin, tobin); //thsh->SetDirectory(NULL); //thsh->SetXTitle("cos( \\Theta)"); //thsh->SetYTitle("Counts"); //Float_t dbin = (tobin-frombin)/nbins; //Float_t lbin = frombin +dbin*0.5; //for (Int_t j=1; j<=nbins; j++) //{ // thsh->Fill(lbin,1.0); // lbin += dbin; //} Int_t refcolumn = 0; MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); TH1F *thsh = new TH1F(); thsh->SetNameTitle("thsh","target distribution"); MH::SetBinning(thsh, binscostheta); Int_t nbins = thsh->GetNbinsX(); Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; for (Int_t j=1; j<=nbins; j++) { thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); } gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " << refcolumn << ", " << nmaxevts << endl; matrixg.EnableGraphicalOutput(); Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); if ( !matrixok ) { gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; return; } gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " << refcolumn << ", " << nmaxevts << endl; matrixh.EnableGraphicalOutput(); matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); delete thsh; if ( !matrixok ) { gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; return; } //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ // write out look-alike events gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Write out look=alike events" << endl; //------------------------------------------- // "gammas" gLog << "Gammas :" << endl; matrixg.Print("cols"); TFile writeg(outNameGammas, "RECREATE", ""); matrixg.Write(); gLog << "" << endl; gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " << outNameGammas << endl; //------------------------------------------- // "hadrons" gLog << "Hadrons :" << endl; matrixh.Print("cols"); TFile writeh(outNameHadrons, "RECREATE", ""); matrixh.Write(); gLog << "" << endl; gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " << outNameHadrons << endl; } //********** end of creating matrices of look-alike 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.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) { *fLog << err << dbginf << "MRanForest not found... aborting." << endl; return kFALSE; } fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); if (!fRanTree) { *fLog << 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"); // list of tasks for the loop over the trees tlist2.AddToList(&rfgrow2); tlist2.AddToList(&rfwrite2); //------------------- // Execute tree loop // MEvtLoop treeloop; treeloop.SetParList(&plist2); if ( !treeloop.Eventloop() ) return; tlist2.PrintStatistics(0, kTRUE); // get adresses of objects which are used in the next eventloop fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); if (!fRanForest) { *fLog << err << dbginf << "MRanForest not found... aborting." << endl; return kFALSE; } fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); if (!fRanTree) { *fLog << 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(); TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; //....................................................................... // calculate hadronnes for method of RANDOM FOREST TString hadRFName = "HadRF"; 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("HadNN", "Events"); write.AddContainer("HadSC", "Events"); write.AddContainer(hadRFName, "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(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"); 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); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&rfcalc); tliston.AddToList(&fillranfor); tliston.AddToList(&fillhadrf); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); tliston.AddToList(&write); tliston.AddToList(&contfinalgh); 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("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(); DeleteBinnings(&pliston); } gLog << "Macro CT1Analysis : End of Job B_RF_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 NEAREST NEIGHBORS and for the SUOERCUTS // - produce Neyman-Pearson plots if (JobC) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job C" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobC = " << JobC << endl; // name of input data file TString filenameData = outPath; filenameData += "ON"; 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 hadNNName = "HadNN"; 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(hadNNName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); fillhadnn.SetName("HhadNN"); 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(hadNNName); 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(&fillhadnn); 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("hadNN", "MHHadronness")->DrawClone(); 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 << 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("NN"); //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.20; 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 fillhadnn("hadNN[MHHadronness]", "HadNN"); fillhadnn.SetName("HhadNN"); MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); fillhadsc.SetName("HhadSC"); MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); fillhadrf.SetName("HhadRF"); 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); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&contfinalgh); tliston.AddToList(&contfinal); tliston.AddToList(&fillhadnn); tliston.AddToList(&fillhadsc); tliston.AddToList(&fillhadrf); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); //***************************** //------------------------------------------- // 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("hadNN", "MHHadronness")->DrawClone(); 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(); DeleteBinnings(&pliston); gLog << "Macro CT1Analysis : End of Job D" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- // Job E_EST_UP //============ // - read MC2.root file // - select g/h separation method XX // - optimize energy estimator for events passing final cuts // - write parameters of energy estimator onto file "energyest_XX.root" // // - read ON2.root and MC2.root files // - update input root file with the estimated energy // (ON_XX2.root, MC_XX2.root) if (JobE_EST_UP) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = " << JobE_EST_UP << ", " << WESTUP << endl; TString typeON = "ON"; TString typeMC = "MC"; TString ext = "3.root"; TString extout = "4.root"; //------------------------------ // name of MC file to be used for optimizing the energy estimator TString filenameOpt(outPath); filenameOpt += typeMC; filenameOpt += ext; gLog << "filenameOpt = " << filenameOpt << endl; //------------------------------ // selection of g/h separation method // and definition of final selections TString XX("NN"); //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.40; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; gLog << "Maximum values of hadronness, |ALPHA| and DIST = " << maxhadronness << ", " << maxalpha << ", " << maxdist << 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 ON file to be updated TString filenameON(outPath); filenameON += typeON; filenameON += ext; gLog << "filenameON = " << filenameON << endl; // name of MC file to be updated TString filenameMC(outPath); filenameMC += typeMC; filenameMC += ext; gLog << "filenameMC = " << filenameMC << endl; //------------------------------ // name of updated ON file TString filenameONup(outPath); filenameONup += typeON; filenameONup += "_"; filenameONup += XX; filenameONup += extout; gLog << "filenameONup = " << filenameONup << endl; // name of updated MC file TString filenameMCup(outPath); filenameMCup += typeMC; filenameMCup += "_"; filenameMCup += XX; filenameMCup += extout; gLog << "filenameMCup = " << filenameMCup << endl; //----------------------------------------------------------- TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; TString fHilNameSrc = "MHillasSrc"; TString fImgParName = "MNewImagePar"; //=========================================================== // // Optimization of energy estimator // TString inpath(""); TString outpath(""); Int_t howMany = 2000; CT1EEst(inpath, filenameOpt, outpath, energyParName, fHilName, fHilNameSrc, fhadronnessName, howMany, maxhadronness, maxalpha, maxdist); //----------------------------------------------------------- // // Read in parameters of energy estimator // gLog << "================================================================" << endl; gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '" << energyParName << "'" << endl; TFile enparam(energyParName); MMcEnergyEst mcest("MMcEnergyEst"); mcest.Read("MMcEnergyEst"); enparam.Close(); TArrayD parA(5); TArrayD parB(7); for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameON); read.DisableAutoScheme(); //--------------------------- // calculate estimated energy MEnergyEstParam eest2(fHilName); eest2.Add(fHilNameSrc); eest2.SetCoeffA(parA); eest2.SetCoeffB(parB); //....................................................................... MWriteRootFile write(filenameONup); 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"); //----------------------------------------------------------------- MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(fhadronnessName); MContinue contfinal(&selfinal); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&eest2); tliston.AddToList(&write); tliston.AddToList(&contfinal); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); Int_t maxevents = -1; //Int_t maxevents = 1000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); DeleteBinnings(&pliston); //--------------------------------------------------- //--------------------------------------------------- // Update MC data // gLog << "============================================================" << endl; gLog << "Macro CT1Analysis.C : update file '" << filenameMC << "'" << endl; MTaskList tlistmc; MParList plistmc; //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameMC); read.DisableAutoScheme(); //--------------------------- // calculate estimated energy MEnergyEstParam eest2(fHilName); eest2.Add(fHilNameSrc); eest2.SetCoeffA(parA); eest2.SetCoeffB(parB); //....................................................................... MWriteRootFile write(filenameMCup); 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"); //----------------------------------------------------------------- MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(fhadronnessName); MContinue contfinal(&selfinal); //***************************** // entries in MParList plistmc.AddToList(&tlistmc); InitBinnings(&plistmc); //***************************** // entries in MTaskList tlistmc.AddToList(&read); tlistmc.AddToList(&eest2); tlistmc.AddToList(&write); tlistmc.AddToList(&contfinal); //***************************** //------------------------------------------- // Execute event loop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plistmc); evtloop.SetProgressBar(&bar); Int_t maxevents = -1; //Int_t maxevents = 1000; if ( !evtloop.Eventloop(maxevents) ) return; tlistmc.PrintStatistics(0, kTRUE); DeleteBinnings(&plistmc); //========== end update ============================================ } enparam.Close(); gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- // Job F_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 (JobF_XX) { gLog << "=====================================================" << endl; gLog << "Macro CT1Analysis : Start of Job F_XX" << endl; gLog << "" << endl; gLog << "Macro CT1Analysis : JobF_XX, WXX = " << JobF_XX << ", " << WXX << 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 = "2.root"; TString extout = "3.root"; //------------------------------ // selection of g/h separation method // and definition of final selections //TString XX("NN"); 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.40; 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 calculating the eff. collection areas TString filenameArea(outPath); filenameArea += typeMC; filenameArea += "_"; filenameArea += XX; 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 += "_"; filenameData += XX; 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; //==================================================================== 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 = new MHMcCT1CollectionArea(); 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 = parlist.FindObject("MHMcCT1CollectionArea"); collarea->CalcEfficiency(); collarea->DrawClone("lego"); //--------------------------------------------- // Write histograms to a file // TFile f(collareaName, "RECREATE"); collarea->GetHist()->Write(); collarea->GetHAll()->Write(); collarea->GetHSel()->Write(); f.Close(); gLog << "Calculation of effective collection areas done" << endl; gLog << "-----------------------------------------------" << 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(); //....................................................................... MWriteRootFile write(filenameDataout); 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"); 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); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&contfinalgh); if (WXX) 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); 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("hadNN", "MHHadronness")->DrawClone(); 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(); DeleteBinnings(&pliston); gLog << "Macro CT1Analysis : End of Job F_XX" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- }