#include "CT1EEst.C" void InitBinnings(MParList *plist) { cout << "InitBinnings" << endl; 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 *binsht = new MBinning("BinningHeadTail"); binsht->SetEdges(50, -1.5, 1.5); plist->AddToList(binsht); 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); Int_t nbins = 8; TArrayD xbins(nbins+1); xbins[0] = 0.0; xbins[1] = 7.5; xbins[2] = 17.5; xbins[3] = 23.5; xbins[4] = 29.5; xbins[5] = 35.5; xbins[6] = 42.0; xbins[7] = 55.0; xbins[8] = 68.0; MBinning *binth = new MBinning("BinningTheta"); binth->SetEdges(xbins); plist->AddToList(binth); MBinning *binsdiff = new MBinning("BinningDiffsigma2"); binsdiff->SetEdges(100, -5.0, 20.0); plist->AddToList(binsdiff); } void DeleteBinnings(MParList *plist) { cout << "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("BinningHeadTail"); 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("BinningDiffsigma2"); if (!bin) delete bin; } void CT1Analysis() { //----------------------------------------------- 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/"; //************************************************************************ // Job A_ON : read ON data // - generate sigmabar vs. Theta plot; // - write root file for ON data (ON1.root); // - generate hadron matrix for g/h separation Bool_t JobA_ON = kFALSE; Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ? Bool_t WLookON = kFALSE; // write out look-alike events for hadrons ? 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); // - generate gamma matrix for g/h separation) Bool_t JobA_MC = kFALSE; Bool_t WLookMC = kFALSE; // write out look-alike events for gammas ? Bool_t WMC1 = kFALSE; // write out root file MC1.root ? // Job B_NN_UP : read ON1.root (or MC1.root) file // - read hadron and gamma matrix // - calculate hadroness for method of Nearest Neighbors // and for the supercuts // - update the input files with the hadronesses (ON2.root or MC2.root) Bool_t JobB_NN_UP = kFALSE; Bool_t WNN = kFALSE; // write root file ? // Job B_EST_UP : // - read MC2.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 ON2.root and MC2.root files with estimated energy // (ON_XX2.root and MC_XX2.root) Bool_t JobB_EST_UP = kFALSE; Bool_t WESTUP = kFALSE; // update root file ? // Job C_NN : read ON1.root and MC1.root files // - read look-alike events for hadrons and gammas // - calculate hadronness for method of NEAREST NEIGHBOR // - produce Neyman-Pearson plot Bool_t JobC_NN = kFALSE; // Job C_RF : read ON1.root and MC1.root files // - read look-alike events for hadrons and gammas // - calculate hadronness for method of RANDOM FOREST // - produce Neyman-Pearson plot Bool_t JobC_RF = kFALSE; // Job D_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 JobD_XX = kTRUE; Bool_t WXXD = 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 look-alike events (matrix_hadrons_ON.root) // (to be used later together with the corresponding MC gamma file // for the g/h separation in the analysis of the ON 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) { cout << "=====================================================" << endl; cout << "Macro CT1Analysis : Start of Job A_ON" << endl; cout << "" << endl; cout << "Macro CT1Analysis : JobA_ON, WHistON, WLookON, WON1 = " << JobA_ON << ", " << WHistON << ", " << WLookON << ", " << WON1 << endl; TString typeInput = "ON"; cout << "typeInput = " << typeInput << endl; // name of input root file TString filenamein(onfile); // name of output root file TString outNameImage = outPath; outNameImage += typeInput; outNameImage += "1.root"; //----------------------------------------------------------- MTaskList tliston; MParList pliston; MFilterList flist; char *sourceName = "MSrcPosCam"; MSrcPosCam source(sourceName); //------------------------------------------- // create the tasks which should be executed // MCT1ReadPreProc read(filenamein); MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", "MCT1PointingCorrCalc"); MBlindPixelCalc blind; MFCT1SelBasic selbasic; MContinue contbasic(&selbasic); MSigmabarCalc sigbarcalc; MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); MImgCleanStd clean; // calculate also extended image parameters TString fHilName = "MHillas"; MHillasExt hext(fHilName); MHillasCalc hcalc(fHilName); TString fHilSrcName = "MHillasSrc"; MHillasSrc hilsrc(fHilSrcName); MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", ""); csrc1.SetInput(fHilName); TString fnewparName = "MNewImagePar"; MNewImagePar newpar(fnewparName); MNewImageParCalc newimage(sourceName, fnewparName); newimage.SetInput(fHilName); MFCT1SelStandard selstandard(fHilName, fHilSrcName); MContinue contstandard(&selstandard); MFillH hfill1("MHHillas", fHilName); hfill1.SetTitle("Task to plot the source independent Hillas parameters"); MFillH hfill2("MHStarMap", fHilName); hfill2.SetTitle("Task to plot the star map (points along main axis of shower)"); MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); MFillH hfill3("MHHillasExt", fHilSrcName); hfill3.SetTitle("Task to plot the extended Hillas parameters"); MFillH hfill4("MHHillasSrc", fHilSrcName); hfill4.SetTitle("Task to plot the source dependent Hillas parameters"); MFillH hfill5("MHNewImagePar", fnewparName); hfill5.SetTitle("Task to plot the new image parameters"); //+++++++++++++++++++++++++++++++++++++++++++++++++++ // prepare writing of look-alike events (for hadrons) const char* mtxName = "MatrixHadrons"; Int_t howMany = 500000; TString outName = outPath; outName += "matrix_hadrons_"; outName += typeInput; outName += ".root"; cout << "" << endl; cout << "========================================================" << endl; cout << "Matrix for (hadrons)" << endl; cout << "input file = " << filenamein<< endl; cout << "matrix name = " << mtxName << endl; cout << "max. no.of look-alike events = " << howMany << endl; cout << "name of output root file = " << outName << endl; cout << "" << endl; // matrix limitation for look-alike events (approximate number) MFEventSelector selector("", ""); selector.SetNumSelectEvts(howMany); // // --- setup the matrix and add it to the parlist // MHMatrix *pmatrix = new MHMatrix(mtxName); MHMatrix& matrix = *pmatrix; matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); matrix.AddColumn("MSigmabar.fSigmabar"); matrix.AddColumn("log10(MHillas.fSize)"); matrix.AddColumn("MHillasSrc.fDist"); matrix.AddColumn("MHillas.fWidth"); matrix.AddColumn("MHillas.fLength"); matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)"); matrix.AddColumn("MHillas.fConc"); matrix.AddColumn("MNewImagePar.fLeakage1"); matrix.Print("cols"); MFillH fillmatg(mtxName); fillmatg.SetFilter(&flist); //+++++++++++++++++++++++++++++++++++++++++++++++++++ if (WON1) { MWriteRootFile write(outNameImage); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); write.AddContainer("MMcEvt", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); } //***************************** // entries in MFilterList flist.AddToList(&selector); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&source); pliston.AddToList(&sigbarcalc); pliston.AddToList(&hext); pliston.AddToList(&hilsrc); pliston.AddToList(&newpar); pliston.AddToList(&hhilext); pliston.AddToList(pmatrix); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&pointcorr); tliston.AddToList(&blind); tliston.AddToList(&contbasic); tliston.AddToList(&sigbarcalc); tliston.AddToList(&fillsigtheta); tliston.AddToList(&clean); tliston.AddToList(&hcalc); tliston.AddToList(&csrc1); tliston.AddToList(&newimage); tliston.AddToList(&flist); tliston.AddToList(&fillmatg); 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.SetProgressBar(&bar); //evtloop.Write(); //Int_t maxevents = 1000000000; Int_t maxevents = 20000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); DeleteBinnings(&pliston); //------------------------------------------- // Display the histograms // TObject *c; c = ((MHSigmaTheta*)pliston.FindObject("SigmaTheta", "MHSigmaTheta")) ->DrawClone(); c = pliston.FindObject("MHHillas")->DrawClone(); c = pliston.FindObject("MHHillasExt")->DrawClone(); c = pliston.FindObject("MHHillasSrc")->DrawClone(); c = pliston.FindObject("MHNewImagePar")->DrawClone(); c = pliston.FindObject("MHStarMap")->DrawClone(); //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ // // ---------------------------------------------------------- // Definition of the reference sample of look-alike events // (this is a subsample of the original sample) // ---------------------------------------------------------- // cout << "" << endl; cout << "========================================================" << 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->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; } cout << "" << endl; cout << "========================================================" << endl; cout << "Macro CT1Analysis : define reference sample for hadrons" << endl; cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " << refcolumn << ", " << nmaxevts << endl; matrix.EnableGraphicalOutput(); Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts); if ( !matrixok ) { cout << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; return; } //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ //------------------------------------------- // write out look-alike events (for hadrons) // if (WLookON) { TFile *writejens = new TFile(outName, "RECREATE", ""); matrix.Write(); cout << "" << endl; cout << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " << outName << endl; writejens->Close(); delete writejens; } //------------------------------------------- // Write histograms onto a file if (WHistON) { MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); if (!sigtheta) { cout << "Object 'SigmaTheta' not found" << endl; return; } TH2D *fHSigTh = sigtheta->GetSigmaTheta(); TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta(); TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); TString outNameSigTh = outPath; outNameSigTh += "SigmaTheta_"; outNameSigTh += typeInput; outNameSigTh += ".root"; TFile outfile(outNameSigTh, "RECREATE"); fHSigTh->Write(); fHSigPixTh->Write(); fHDifPixTh->Write(); outfile.Close(); cout << "" << endl; cout << "File " << outNameSigTh << " was written out" << endl; } cout << "Macro CT1Analysis : End of Job A_ON" << endl; cout << "===================================================" << 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 look-alike events (matrix_gammas_MC.root) // (to be used later together with the corresponding hadron file // for the g/h separation in the analysis 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) { cout << "=====================================================" << endl; cout << "Macro CT1Analysis : Start of Job A_MC" << endl; cout << "" << endl; cout << "Macro CT1Analysis : JobA_MC, WLookMC, WMC1 = " << JobA_MC << ", " << WLookMC << ", " << WMC1 << endl; TString typeInput = "MC"; cout << "typeInput = " << typeInput << endl; // name of input root file TString filenamein(mcfile); // name of output root file TString outNameImage = outPath; outNameImage += typeInput; outNameImage += "1.root"; // use for padding sigmabar vs. Theta from ON data TString typeHist = "ON"; cout << "typeHist = " << typeHist << endl; //------------------------------------ // Get the histograms "2D-ThetaSigmabar" // and "3D-ThetaPixSigma" // and "3D-ThetaPixDiff" TString outNameSigTh = outPath; outNameSigTh += "SigmaTheta_"; outNameSigTh += typeHist; outNameSigTh += ".root"; cout << "Reading in file " << outNameSigTh << endl; TFile *infile = new TFile(outNameSigTh); infile->ls(); TH2D *fHSigmaTheta = (TH2D*) gROOT.FindObject("2D-ThetaSigmabar"); if (!fHSigmaTheta) { cout << "Object '2D-ThetaSigmabar' not found on root file" << endl; return; } cout << "Object '2D-ThetaSigmabar' was read in" << endl; TH3D *fHSigmaPixTheta = (TH3D*) gROOT.FindObject("3D-ThetaPixSigma"); if (!fHSigmaPixTheta) { cout << "Object '3D-ThetaPixSigma' not found on root file" << endl; return; } cout << "Object '3D-ThetaPixSigma' was read in" << endl; TH3D *fHDiffPixTheta = (TH3D*) gROOT.FindObject("3D-ThetaPixDiff"); if (!fHSigmaPixTheta) { cout << "Object '3D-ThetaPixDiff' not found on root file" << endl; return; } cout << "Object '3D-ThetaPixDiff' was read in" << endl; //------------------------------------ MTaskList tlist; MParList plist; MFilterList flist; char *sourceName = "MSrcPosCam"; MSrcPosCam source(sourceName); //------------------------------------------- // create the tasks which should be executed // MCT1ReadPreProc read(filenamein); MBlindPixelCalc blind; MFCT1SelBasic selbasic; MContinue contbasic(&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 MPadSchweizer padthomas("MPadSchweizer","Task for the padding (Schweizer)", fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta); padthomas.SetPadFlag(1); //........................................... MSigmabarCalc sigbarcalc; MFillH fillsigtheta ("MCSigmaTheta[MHSigmaTheta]", "MMcEvt"); fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); MImgCleanStd clean; // calculate also extended image parameters TString fHilName = "MHillas"; MHillasExt hext(fHilName); MHillasCalc hcalc(fHilName); TString fHilSrcName = "MHillasSrc"; MHillasSrc hilsrc(fHilSrcName); MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", ""); csrc1.SetInput(fHilName); TString fnewparName = "MNewImagePar"; MNewImagePar newpar(fnewparName); MNewImageParCalc newimage(sourceName, fnewparName); newimage.SetInput(fHilName); MFCT1SelStandard selstandard(fHilName, fHilSrcName); MContinue contstandard(&selstandard); MFillH hfill1("MHHillas", fHilName); hfill1.SetTitle("Task to plot the source independent Hillas parameters"); MFillH hfill2("MHStarMap", fHilName); hfill2.SetTitle("Task to plot the star map (points along main axis of shower)"); MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); MFillH hfill3("MHHillasExt", fHilSrcName); hfill3.SetTitle("Task to plot the extended Hillas parameters"); MFillH hfill4("MHHillasSrc", fHilSrcName); hfill4.SetTitle("Task to plot the source dependent Hillas parameters"); MFillH hfill5("MHNewImagePar", fnewparName); hfill5.SetTitle("Task to plot the new image parameters"); //+++++++++++++++++++++++++++++++++++++++++++++++++++ // prepare writing of look-alike events (for gammas) const char* mtxName = "MatrixGammas"; Int_t howMany = 500000; TString outName = outPath; outName += "matrix_gammas_"; outName += typeInput; outName += ".root"; cout << "" << endl; cout << "========================================================" << endl; cout << "Matrix for (gammas)" << endl; cout << "input file = " << filenamein<< endl; cout << "matrix name = " << mtxName << endl; cout << "max. no.of look-alike events = " << howMany << endl; cout << "name of output root file = " << outName << endl; cout << "" << endl; // matrix limitation for look-alike events (approximate number) MFEventSelector selector("", ""); selector.SetNumSelectEvts(howMany); // // --- setup the matrix and add it to the parlist // MHMatrix *pmatrix = new MHMatrix(mtxName); MHMatrix& matrix = *pmatrix; matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); matrix.AddColumn("MSigmabar.fSigmabar"); matrix.AddColumn("log10(MHillas.fSize)"); matrix.AddColumn("MHillasSrc.fDist"); matrix.AddColumn("MHillas.fWidth"); matrix.AddColumn("MHillas.fLength"); matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)"); matrix.AddColumn("MHillas.fConc"); matrix.AddColumn("MNewImagePar.fLeakage1"); matrix.Print("cols"); MFillH fillmatg(mtxName); fillmatg.SetFilter(&flist); //+++++++++++++++++++++++++++++++++++++++++++++++++++ if (WMC1) { MWriteRootFile write(outNameImage); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); write.AddContainer("MMcEvt", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); } //***************************** // entries in MFilterList flist.AddToList(&selector); //***************************** // entries in MParList plist.AddToList(&tlist); InitBinnings(&plist); plist.AddToList(&source); plist.AddToList(&sigbarcalc); plist.AddToList(&hext); plist.AddToList(&hilsrc); plist.AddToList(&newpar); plist.AddToList(&hhilext); plist.AddToList(pmatrix); //***************************** // entries in MTaskList tlist.AddToList(&read); tlist.AddToList(&blind); tlist.AddToList(&padthomas); tlist.AddToList(&contbasic); tlist.AddToList(&sigbarcalc); tlist.AddToList(&fillsigtheta); tlist.AddToList(&clean); tlist.AddToList(&hcalc); tlist.AddToList(&csrc1); tlist.AddToList(&newimage); tlist.AddToList(&flist); tlist.AddToList(&fillmatg); 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.SetProgressBar(&bar); //evtloop.Write(); //Int_t maxevents = 1000000000; Int_t maxevents = 1000; if ( !evtloop.Eventloop(maxevents) ) return; tlist.PrintStatistics(0, kTRUE); DeleteBinnings(&plist); //------------------------------------------- // Display the histograms // TObject *c; c = ((MHSigmaTheta*)plist.FindObject("MCSigmaTheta", "MHSigmaTheta")) ->DrawClone(); c = plist.FindObject("MHHillas")->DrawClone(); c = plist.FindObject("MHHillasExt")->DrawClone(); c = plist.FindObject("MHHillasSrc")->DrawClone(); c = plist.FindObject("MHNewImagePar")->DrawClone(); c = plist.FindObject("MHStarMap")->DrawClone(); //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ // // ---------------------------------------------------------- // Definition of the reference sample of look-alike events // (this is a subsample of the original sample) // ---------------------------------------------------------- // cout << "" << endl; cout << "========================================================" << 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); 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; } cout << "" << endl; cout << "========================================================" << endl; cout << "Macro CT1Analysis : define reference sample for gammas" << endl; cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " << refcolumn << ", " << nmaxevts << endl; matrix.EnableGraphicalOutput(); Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts); if ( !matrixok ) { cout << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; return; } //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ //------------------------------------------- // write out look-alike events (for gammas) // if (WLookMC) { TFile *writejens = new TFile(outName, "RECREATE", ""); matrix.Write(); cout << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " << outName << endl; writejens->Close(); delete writejens; } cout << "Macro CT1Analysis : End of Job A_MC" << endl; cout << "=========================================================" << 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) { cout << "=====================================================" << endl; cout << "Macro CT1Analysis : Start of Job B_NN_UP" << endl; cout << "" << endl; cout << "Macro CT1Analysis : JobB_NN_UP, WNN = " << JobB_NN_UP << ", " << WNN << endl; //TString typeInput = "ON"; TString typeInput = "MC"; cout << "typeInput = " << typeInput << endl; // name of input root file TString filenameData = outPath; filenameData += typeInput; filenameData += "1.root"; cout << "filenameData = " << filenameData << endl; // name of output root file TString outNameImage = outPath; outNameImage += typeInput; outNameImage += "2.root"; //$$$$$$$$$$$$$$$$$$$$$$$$$$ outNameImage = filenameData; //$$$$$$$$$$$$$$$$$$$$$$$$$$ cout << "outNameImage = " << outNameImage << endl; TString typeMatrixHadrons = "ON"; cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl; //************************************************************************* // read in matrices of look-alike events const char* mtxName = "MatrixGammas"; TString outName = outPath; outName += "matrix_gammas_"; outName += "MC"; outName += ".root"; cout << "" << endl; cout << "========================================================" << endl; cout << "Get matrix for (gammas)" << endl; cout << "matrix name = " << mtxName << endl; cout << "name of root file = " << outName << endl; cout << "" << endl; // read in the object with the name 'mtxName' from file 'outName' // TFile *fileg = new TFile(outName); MHMatrix matrixg; matrixg.Read(mtxName); matrixg.Print("cols"); delete fileg; //***************************************************************** const char* mtxName = "MatrixHadrons"; TString outName = outPath; outName += "matrix_hadrons_"; outName += typeMatrixHadrons; outName += ".root"; cout << "" << endl; cout << "========================================================" << endl; cout << " Get matrix for (hadrons)" << endl; cout << "matrix name = " << mtxName << endl; cout << "name of root file = " << outName << endl; cout << "" << endl; // read in the object with the name mtxName // TFile *fileh = new TFile(outName); MHMatrix matrixh; matrixh.Read(mtxName); matrixh.Print("cols"); delete fileh; //************************************************************************* MTaskList tliston; MParList pliston; //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); TString fHilName = "MHillas"; TString fHilSrcName = "MHillasSrc"; TString fnewparName = "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, fHilSrcName); sccalc.SetHadronnessName(hadSCName); //....................................................................... if (WNN) { //MWriteRootFile write(outNameImage); MWriteRootFile write(outNameImage, "UPDATE"); //write.AddContainer("MRawRunHeader", "RunHeaders"); //write.AddContainer("MTime", "Events"); //write.AddContainer("MMcEvt", "Events"); //write.AddContainer("MSrcPosCam", "Events"); //write.AddContainer("MSigmabar", "Events"); //write.AddContainer("MHillas", "Events"); //write.AddContainer("MHillasSrc", "Events"); //write.AddContainer("MNewImagePar", "Events"); write.AddContainer(hadNNName, "Events"); write.AddContainer(hadSCName, "Events"); } //----------------------------------------------------------------- // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); Float_t hadcut = 0.40; Float_t alphacut = 100.0; MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); selfinalgh.SetCuts(hadcut, alphacut); selfinalgh.SetHadronnessName(hadNNName); MContinue contfinalgh(&selfinalgh); MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); Float_t hadcut = 0.40; Float_t alphacut = 20.0; MFCT1SelFinal selfinal(fHilName, fHilSrcName); selfinal.SetCuts(hadcut, alphacut); selfinal.SetHadronnessName(hadNNName); MContinue contfinal(&selfinal); MFillH hfill1("MHHillas", fHilName); MFillH hfill2("MHStarMap", fHilName); MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); MFillH hfill3("MHHillasExt", fHilSrcName); MFillH hfill4("MHHillasSrc", fHilSrcName); MFillH hfill5("MHNewImagePar", fnewparName); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&matrixg); pliston.AddToList(&matrixh); pliston.AddToList(&hhilext); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&nncalc); tliston.AddToList(&sccalc); tliston.AddToList(&fillhadnn); tliston.AddToList(&fillhadsc); if (WNN) tliston.AddToList(&write); 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); //evtloop.Write(); if ( !evtloop.Eventloop() ) //if ( !evtloop.Eventloop(1000) ) return; tliston.PrintStatistics(0, kTRUE); DeleteBinnings(&pliston); //------------------------------------------- // Display the histograms // TObject *c; c = pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); c = pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); c = pliston.FindObject("MHHillas")->DrawClone(); c = pliston.FindObject("MHHillasExt")->DrawClone(); c = pliston.FindObject("MHHillasSrc")->DrawClone(); c = pliston.FindObject("MHNewImagePar")->DrawClone(); c = pliston.FindObject("MHStarMap")->DrawClone(); cout << "Macro CT1Analysis : End of Job B_NN_UP" << endl; cout << "=======================================================" << endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- // Job B_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 (JobB_EST_UP) { cout << "=====================================================" << endl; cout << "Macro CT1Analysis : Start of Job B_EST_UP" << endl; cout << "" << endl; cout << "Macro CT1Analysis : JobB_EST_UP, WESTUP = " << JobB_EST_UP << ", " << WESTUP << endl; TString typeON = "ON"; TString typeMC = "MC"; TString ext = "3.root"; //------------------------------ // name of MC file to be used for optimizing the energy estimator TString filenameOpt(outPath); filenameOpt += typeMC; filenameOpt += ext; cout << "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; cout << "fhadronnessName = " << fhadronnessName << endl; // maximum values of the hadronness and of |alpha| Float_t maxhadronness = 0.40; Float_t maxalpha = 20.0; cout << "Maximum values of hadronness and |alpha| = " << maxhadronness << ", " << maxalpha << endl; // name of file containing the parameters of the energy estimator TString energyParName(outPath); energyParName += "energyest_"; energyParName += XX; energyParName += ".root"; cout << "energyParName = " << energyParName << endl; //------------------------------ // name of ON file to be updated TString filenameON(outPath); filenameON += typeON; filenameON += ext; cout << "filenameON = " << filenameON << endl; // name of MC file to be updated TString filenameMC(outPath); filenameMC += typeMC; filenameMC += ext; cout << "filenameMC = " << filenameMC << endl; //------------------------------ // name of updated ON file TString filenameONup(outPath); filenameONup += typeON; filenameONup += "_"; filenameONup += XX; filenameONup += ext; cout << "filenameONup = " << filenameONup << endl; // name of updated MC file TString filenameMCup(outPath); filenameMCup += typeMC; filenameMCup += "_"; filenameMCup += XX; filenameMCup += ext; cout << "filenameMCup = " << filenameMCup << endl; //----------------------------------------------------------- TString fHilName = "MHillas"; TString fHilSrcName = "MHillasSrc"; TString fnewparName = "MNewImagePar"; //=========================================================== // // Optimization of energy estimator // TString inpath(""); TString outpath(""); Int_t howMany = 2000; CT1EEst(inpath, filenameOpt, outpath, energyParName, fHilName, fHilSrcName, fhadronnessName, howMany, maxhadronness, maxalpha); //----------------------------------------------------------- // // Read in parameters of energy estimator // TFile enparam(energyParName); TArrayD parA(5); TArrayD parB(7); for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); TString fHilName = "MHillas"; TString fHilSrcName = "MHillasSrc"; TString fnewparName = "MNewImagePar"; Float_t hadcut = 0.36; Float_t alphacut = 100.0; MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); selfinalgh.SetCuts(hadcut, alphacut); selfinalgh.SetHadronnessName(hadName); MContinue contfinalgh(&selfinalgh); MFillH fillhadnn("MHHadronness", hadName); Float_t hadcut = 0.36; Float_t alphacut = 20.0; MFCT1SelFinal selfinal(fHilName, fHilSrcName); selfinal.SetCuts(hadcut, alphacut); selfinal.SetHadronnessName(hadName); MContinue contfinal(&selfinal); MFillH hfill1("MHHillas", fHilName); MFillH hfill2("MHStarMap", fHilName); MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); MFillH hfill3("MHHillasExt", fHilSrcName); MFillH hfill4("MHHillasSrc", fHilSrcName); MFillH hfill5("MHNewImagePar", fnewparName); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&matrixg); pliston.AddToList(&matrixh); pliston.AddToList(&hhilext); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&fgamma); tliston.AddToList(&fhadron); tliston.AddToList(&nncalc); tliston.AddToList(&fillhadnn); 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); //evtloop.Write(); Int_t maxevents = 1000000000; //Int_t maxevents = 35000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); DeleteBinnings(&pliston); //------------------------------------------- // Display the histograms // TObject *c; c = pliston.FindObject("MHHadronness"); if (c) { c->DrawClone(); c->Print(); } c = pliston.FindObject("MHHillas")->DrawClone(); c = pliston.FindObject("MHHillasExt")->DrawClone(); c = pliston.FindObject("MHHillasSrc")->DrawClone(); c = pliston.FindObject("MHNewImagePar")->DrawClone(); c = pliston.FindObject("MHStarMap")->DrawClone(); cout << "Macro CT1Analysis : End of Job C_NN" << endl; cout << "===================================================" << endl; } //--------------------------------------------------------------------- // Job C_RF //========= // read ON1 and MC1 data files // // - calculate hadronness for method of Random Forest // - produce Neyman-Pearson plot if (JobC_RF) { cout << "=====================================================" << endl; cout << "Macro CT1Analysis : Start of Job C_RF" << endl; cout << "" << endl; cout << "Macro CT1Analysis : JobC_RF = " << JobC_RF << endl; TString typeMC = "MC"; cout << "typeMC = " << typeMC << endl; TString typeData = "ON"; cout << "typeData = " << typeData << endl; TString typeMatrixHadrons = "ON"; cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl; //************************************************************************* // read in matrices of look-alike events const char* mtxName = "MatrixGammas"; TString outName = outPath; outName += "matrix_gammas_"; outName += "MC"; outName += ".root"; cout << "" << endl; cout << "========================================================" << endl; cout << "Get matrix for (gammas)" << endl; cout << "matrix name = " << mtxName << endl; cout << "name of root file = " << outName << endl; cout << "" << endl; // read in the object with the name 'mtxName' from file 'outName' // TFile *fileg = new TFile(outName); MHMatrix matrixg; matrixg.Read(mtxName); matrixg.Print("cols"); delete fileg; //***************************************************************** const char* mtxName = "MatrixHadrons"; TString outName = outPath; outName += "matrix_hadrons_"; outName += typeMatrixHadrons; outName += ".root"; cout << "" << endl; cout << "========================================================" << endl; cout << " Get matrix for (hadrons)" << endl; cout << "matrix name = " << mtxName << endl; cout << "name of root file = " << outName << endl; cout << "" << endl; // read in the object with the name mtxName // TFile *fileh = new TFile(outName); MHMatrix matrixh; matrixh.Read(mtxName); matrixh.Print("cols"); delete fileh; //***************************************************************** //----------------------------------------------------------------- // grow the trees of the random forest (event loop = tree loop) cout << "Macro CT1Analysis : start growing trees" << endl; MTaskList tlist2; MParList plist2; plist2.AddToList(&tlist2); plist2.AddToList(&matrixg); plist2.AddToList(&matrixh); MRanForestGrow rfgrow2; rfgrow2.SetNumTrees(100); rfgrow2.SetNumTry(4); rfgrow2.SetNdSize(10); TString outRF = outPath; outRF += "RF.root"; 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 MRanForest *fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); if (!fRanForest) { *fLog << err << dbginf << "MRanForest not found... aborting." << endl; return kFALSE; } MRanTree *fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); if (!fRanTree) { *fLog << err << dbginf << "MRanTree not found... aborting." << endl; return kFALSE; } //----------------------------------------------------------------- MTaskList tliston; MParList pliston; //------------------------------------------- // create the tasks which should be executed // TString filenameData = outPath; filenameData += typeData; filenameData += "1.root"; TString filenameMC = outPath; filenameMC += typeMC; filenameMC += "1.root"; cout << "filenameData = " << filenameData << endl; cout << "filenameMC = " << filenameMC << endl; MReadMarsFile read("Events", filenameMC); read.AddFile(filenameData); read.DisableAutoScheme(); MFParticleId fgamma ("MMcEvt", '=', kGAMMA); MFParticleId fhadron("MMcEvt", '!', kGAMMA); //....................................................................... // calculate hadronnes for method of Random Forest TString hadName = "HadRF"; MRanForestCalc rfcalc; rfcalc.SetHadronnessName(hadName); //....................................................................... // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); TString fHilName = "MHillas"; TString fHilSrcName = "MHillasSrc"; TString fnewparName = "MNewImagePar"; Float_t hadcut = 0.36; Float_t alphacut = 100.0; MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); selfinalgh.SetCuts(hadcut, alphacut); selfinalgh.SetHadronnessName(hadName); MContinue contfinalgh(&selfinalgh); MFillH fillhadrf("MHHadronness", hadName); MFillH fillhrf("MHRanForest"); Float_t hadcut = 0.36; Float_t alphacut = 20.0; MFCT1SelFinal selfinal(fHilName, fHilSrcName); selfinal.SetCuts(hadcut, alphacut); selfinal.SetHadronnessName(hadName); MContinue contfinal(&selfinal); MFillH hfill1("MHHillas", fHilName); MFillH hfill2("MHStarMap", fHilName); MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); MFillH hfill3("MHHillasExt", fHilSrcName); MFillH hfill4("MHHillasSrc", fHilSrcName); MFillH hfill5("MHNewImagePar", fnewparName); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(fRanForest); pliston.AddToList(fRanTree); pliston.AddToList(&hhilext); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&fgamma); tliston.AddToList(&fhadron); tliston.AddToList(&rfcalc); tliston.AddToList(&fillhadrf); tliston.AddToList(&fillhrf); 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); //evtloop.Write(); Int_t maxevents = 1000000000; //Int_t maxevents = 35000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); DeleteBinnings(&pliston); //------------------------------------------- // Display the histograms // TObject *c; c = pliston.FindObject("MHRanForest")->DrawClone(); c = pliston.FindObject("MHHadronness"); if (c) { c->DrawClone(); c->Print(); } c = pliston.FindObject("MHHillas")->DrawClone(); c = pliston.FindObject("MHHillasExt")->DrawClone(); c = pliston.FindObject("MHHillasSrc")->DrawClone(); c = pliston.FindObject("MHNewImagePar")->DrawClone(); c = pliston.FindObject("MHStarMap")->DrawClone(); cout << "Macro CT1Analysis : End of Job C_RF" << endl; cout << "===================================================" << endl; } //--------------------------------------------------------------------- // Job D_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 (JobD_XX) { cout << "=====================================================" << endl; cout << "Macro CT1Analysis : Start of Job D_XX" << endl; cout << "" << endl; cout << "Macro CT1Analysis : JobD_XX, WXXD = " << JobD_XX << ", " << WXXD << endl; // type of data to be analysed TString typeData = "ON"; //TString typeData = "OFF"; //TString typeData = "MC"; cout << "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("NN"); //TString XX("SC"); //TString XX("RF"); TString fhadronnessName("Had"); fhadronnessName += XX; cout << "fhadronnessName = " << fhadronnessName << endl; // maximum values of the hadronness and of |alpha| Float_t maxhadronness = 0.40; Float_t maxalpha = 20.0; cout << "Maximum values of hadronness and |alpha| = " << maxhadronness << ", " << maxalpha << endl; //------------------------------ // name of MC file to be used for calculating the eff. collection areas TString filenameArea(outPath); filenameArea += typeMC; filenameArea += "_"; filenameArea += XX; filenameArea += ext; cout << "filenameArea = " << filenameArea << endl; //------------------------------ // name of file containing the eff. collection areas TString collareaName(outPath); collareaName += "area_"; collareaName += XX; collareaName += ".root"; cout << "collareaName = " << collareaName << endl; //------------------------------ // name of data file to be analysed TString filenameData(outPath); filenameData += typeData; filenameData += "_"; filenameData += XX; filenameData += ext; cout << "filenameData = " << filenameData << endl; //------------------------------ // name of output data file (after the final cuts) TString filenameDataout(outPath); filenameDataout += typeData; filenameDataout += "_"; filenameDataout += XX; filenameDataout += extout; cout << "filenameDataout = " << filenameDataout << endl; //==================================================================== cout << "-----------------------------------------------" << endl; cout << "Start calculation of effective collection areas" << endl; MParList parlist; MTaskList tasklist; //--------------------------------------- // Setup the tasks to be executed // MReadMarsFile reader("Events", filenameArea); MFCT1SelFinal cuthadrons; cuthadrons.SetHadronnessName(fhadronnessName); cuthadrons.SetCuts(maxhadronness, maxalpha); MContinue conthadrons(&cuthadrons); MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea("MHMcCT1CollectionArea","",30,2.,5.); MMcCT1CollectionAreaCalc effi; //******************************** // entries in MParList parlist.AddToList(&tasklist); parlist.AddToList(collarea); //******************************** // entries in MTaskList tasklist.AddToList(&reader); tasklist.AddToList(&conthadrons); tasklist.AddToList(&effi); //******************************** //----------------------------------------- // Execute event loop // MEvtLoop magic; magic.SetParList(&parlist); MProgressBar bar; magic.SetProgressBar(&bar); if (!magic.Eventloop()) return; tasklist.PrintStatistics(0, kTRUE); // Display the histograms // parlist.FindObject("MHMcCT1CollectionArea")->DrawClone("lego"); //--------------------------------------------- // Write histograms to a file // TFile f(collareaName, "RECREATE"); collarea->GetHist()->Write(); collarea->GetHAll()->Write(); collarea->GetHSel()->Write(); f.Close(); cout << "Calculation of effective collection areas done" << endl; cout << "-----------------------------------------------" << endl; //------------------------------------------------------------------ //************************************************************************* // // Analyse the data // MTaskList tliston; MParList pliston; TString hadName("Had"); hadName += XX; TString fHilName = "MHillas"; TString fHilSrcName = "MHillasSrc"; TString fnewparName = "MNewImagePar"; //------------------------------------------- // create the tasks which should be executed // MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); //....................................................................... if (WXXD) { MWriteRootFile write(filenameDataout); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); write.AddContainer("MMcEvt", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); write.AddContainer("HadNN", "Events"); write.AddContainer("HadSC", "Events"); write.AddContainer("MEnergyEst", "Events"); } //----------------------------------------------------------------- // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); selfinalgh.SetCuts(maxhadronness, 100.0); selfinalgh.SetHadronnessName(fhadronnessName); MContinue contfinalgh(&selfinalgh); MFillH fillhadnn("MHadNN[MHHadronness]", "HadNN"); MFillH fillhadsc("MHadSC[MHHadronness]", "HadSC"); //MFillH fillhadrf("MHadRF[MHHadronness]", "HadRF"); MFillH hfill1("MHHillas", fHilName); MFillH hfill2("MHStarMap", fHilName); MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); MFillH hfill3("MHHillasExt", fHilSrcName); MFillH hfill4("MHHillasSrc", fHilSrcName); MFillH hfill5("MHNewImagePar", fnewparName); MFCT1SelFinal selfinal(fHilName, fHilSrcName); selfinal.SetCuts(maxhadronness, maxalpha); selfinal.SetHadronnessName(fhadronnessName); MContinue contfinal(&selfinal); //***************************** // entries in MParList pliston.AddToList(&tliston); InitBinnings(&pliston); pliston.AddToList(&hhilext); //***************************** // entries in MTaskList tliston.AddToList(&read); tliston.AddToList(&contfinalgh); if (WXXD) tliston.AddToList(&write); 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); //evtloop.Write(); Int_t maxevents = 1000000000; //Int_t maxevents = 10000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); DeleteBinnings(&pliston); //------------------------------------------- // Display the histograms // TObject *c; c = pliston.FindObject("MHadNN", "MHHadronness"); if (c) { c->DrawClone(); c->Print(); } c = pliston.FindObject("MHadRF", "MHHadronness"); if (c) { c->DrawClone(); c->Print(); } c = pliston.FindObject("MHadSC", "MHHadronness"); if (c) { c->DrawClone(); c->Print(); } c = pliston.FindObject("MHHillas")->DrawClone(); c = pliston.FindObject("MHHillasExt")->DrawClone(); c = pliston.FindObject("MHHillasSrc")->DrawClone(); c = pliston.FindObject("MHNewImagePar")->DrawClone(); c = pliston.FindObject("MHStarMap")->DrawClone(); cout << "Macro CT1Analysis : End of Job D_XX" << endl; cout << "=======================================================" << endl; } //--------------------------------------------------------------------- }