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