Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2074) +++ trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2074) @@ -0,0 +1,2388 @@ + +#include "CT1EEst.C" + +void InitBinnings(MParList *plist) +{ + cout << "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 *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); + + 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 *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() +{ + 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/optionC/"; + + // path for output from Mars + TString outPath = "~magican/ct1test/wittek/marsoutput/optionC/"; + + + + //************************************************************************ + + // 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 = kTRUE; + Bool_t WHistON = kTRUE; // write out histogram sigmabar vs. Theta ? + Bool_t WLookON = kTRUE; // write out look-alike events for hadrons ? + Bool_t WON1 = kTRUE; // 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 = kFALSE; + 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; + blind.SetUseBlindPixels(); + + MFCT1SelBasic selbasic; + MContinue contbasic(&selbasic); + + MFillH fillblind("MHBlindPixels", "MBlindPixels"); + + MSigmabarCalc sigbarcalc; + + MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); + fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); + + 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); + + + MHHillas hh; + hh.SetMmScale(kFALSE); + 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 hext; + hext.SetMmScale(kFALSE); + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetTitle("Task to plot the extended Hillas parameters"); + + MHHillasSrc hsrc; + hsrc.SetMmScale(kFALSE); + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetTitle("Task to plot the source dependent Hillas parameters"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetTitle("Task to plot the new image parameters"); + // -------------------------------------------------- + + 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); + + + //+++++++++++++++++++++++++++++++++++++++++++++++++++ + // 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)*(MHillasExt.fM3Long)"); + matrix.AddColumn("MNewImagePar.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("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + } + + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + //MF daniel( "(MRawRunHeader.fRunNumber<13167)||(MRawRunHeader.fRunNumber>13167)" ); + //MContinue contdaniel(&daniel); + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + + //***************************** + // entries in MFilterList + + flist.AddToList(&selector); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&source); + pliston.AddToList(pmatrix); + + + //***************************** + // 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(&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 = 1000; + 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("MHBlindPixels")->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->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; + } + + 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); + delete thsh; + 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; + + //Int_t null = 0; + //if (null == 0) goto skip0; + + //------------------------------------ + // 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; + + //------------------------------------ + skip0: + + MTaskList tlist; + MParList plist; + MFilterList flist; + + char *sourceName = "MSrcPosCam"; + MSrcPosCam source(sourceName); + + //------------------------------------------- + // create the tasks which should be executed + // + + MCT1ReadPreProc read(filenamein); + + MBlindPixelCalc blind; + blind.SetUseBlindPixels(); + + MFCT1SelBasic selbasic; + MContinue contbasic(&selbasic); + + MFillH fillblind("MHBlindPixels", "MBlindPixels"); + + // 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)"); + padthomas.SetHistograms(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; + + // 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); + + + MHHillas hh; + hh.SetMmScale(kFALSE); + 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 hext; + hext.SetMmScale(kFALSE); + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetTitle("Task to plot the extended Hillas parameters"); + + MHHillasSrc hsrc; + hsrc.SetMmScale(kFALSE); + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetTitle("Task to plot the source dependent Hillas parameters"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetTitle("Task to plot the new image parameters"); + // -------------------------------------------------- + + 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); + + + //+++++++++++++++++++++++++++++++++++++++++++++++++++ + // 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)*(MHillasExt.fM3Long)"); + matrix.AddColumn("MNewImagePar.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("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 MFilterList + + flist.AddToList(&selector); + + skip: + + //***************************** + // entries in MParList + + plist.AddToList(&tlist); + InitBinnings(&plist); + + plist.AddToList(&source); + plist.AddToList(pmatrix); + + + //***************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&blind); + tlist.AddToList(&padthomas); + + tlist.AddToList(&contbasic); + tlist.AddToList(&fillblind); + tlist.AddToList(&sigbarcalc); + tlist.AddToList(&fillsigtheta); + tlist.AddToList(&clean); + + tlist.AddToList(&hcalc); + tlist.AddToList(&hsrccalc); + + 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 + // + + MHSigmaTheta* cc = (MHSigmaTheta*)plist.FindObject("MCSigmaTheta", "MHSigmaTheta"); + if (!cc) + { + cout << "MCSigmaTheta[MHSigmaTheta] object not found" << endl; + return; + } + cc ->DrawClone(); + + TObject *c; + c = plist.FindObject("MHBlindPixels")->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); + 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; + } + + 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); + delete thsh; + 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 maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadNNName); + MContinue contfinalgh(&selfinalgh); + + MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + + MFCT1SelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadNNName); + MContinue contfinal(&selfinal); + + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt"); + hhilext.SetHillasName(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; + Float_t maxdist = 10.0; + cout << "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"; + 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 maxhadronness = 0.36; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadName); + MContinue contfinalgh(&selfinalgh); + + MFillH fillhadnn("MHHadronness", hadName); + + MFCT1SelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadName); + MContinue contfinal(&selfinal); + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt"); + hhilext.SetHillasName(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 maxhadronness = 0.36; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadName); + MContinue contfinalgh(&selfinalgh); + + MFillH fillhadrf("MHHadronness", hadName); + MFillH fillhrf("MHRanForest"); + + MFCT1SelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadName); + MContinue contfinal(&selfinal); + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt"); + hhilext.SetHillasName(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, |ALPHA| and DIST + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + cout << "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; + 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); + reader.DisableAutoScheme(); + + MFCT1SelFinal cuthadrons; + cuthadrons.SetHadronnessName(fhadronnessName); + cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); + + MContinue conthadrons(&cuthadrons); + + MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea(); + + MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); + + //******************************** + // 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(); + + + 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, maxdist); + 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"); + hhilext.SetHillasName(fHilName); + + MFillH hfill3("MHHillasExt", fHilSrcName); + MFillH hfill4("MHHillasSrc", fHilSrcName); + MFillH hfill5("MHNewImagePar", fnewparName); + + MFCT1SelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + 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 + // + + MHHadronness *c; + c = (MHHadronness*)pliston.FindObject("MHadNN", "MHHadronness"); + if (c) + { + c->DrawClone(); + c->Print(); + } + c = (MHHadronness*)pliston.FindObject("MHadRF", "MHHadronness"); + if (c) + { + c->DrawClone(); + c->Print(); + } + c = (MHHadronness*)pliston.FindObject("MHadSC", "MHHadronness"); + if (c) + { + c->DrawClone(); + c->Print(); + } + + TObject *cc; + cc = pliston.FindObject("MHHillas")->DrawClone(); + cc = pliston.FindObject("MHHillasExt")->DrawClone(); + cc = pliston.FindObject("MHHillasSrc")->DrawClone(); + cc = pliston.FindObject("MHNewImagePar")->DrawClone(); + cc = pliston.FindObject("MHStarMap")->DrawClone(); + + + cout << "Macro CT1Analysis : End of Job D_XX" << endl; + cout << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + +} +