
#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] = 
                       {0.0, 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;

        if (!plist)
	{
	  gLog << "Deletebinnins : MParlist no longer existing" << endl;
          return;
	}

        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;

        gLog << "exit DeleteBinnings" << endl;
}

//************************************************************************
void ONOFFCT1Analysis()
{
      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/optionC/";

      // path for output from Mars
      TString outPath = "~magican/ct1test/wittek/marsoutput/optionC/";

      //-----------------------------------------------

      TEnv env("macros/CT1env.rc");
      Bool_t printEnv = kFALSE;

    //************************************************************************

    // Job A : 
    //  - produce MHSigmaTheta plots for ON and OFF data
    //  - write out (or read in) these MHSigmaTheta plots
    //  - read ON (or OFF or MC) data
    //  - pad the events; 
    //  - write root file for ON (or OFF or MC) data (ON1.root, ...);

    Bool_t JobA    = kTRUE;  
    Bool_t WPad    = kTRUE;   // write out padding histograms ?
    Bool_t RPad    = kFALSE;   // read in padding histograms ?
    Bool_t Wout    = kFALSE;   // write out root file ON1.root 
                               // (or OFF1.root or MC1.root)?

    // Job B_NN_UP : read ON1.root (OFF1.root or MC1.root) file 
    //  - depending on RlookNN : create (or read in) hadron and gamma matrix
    //  - calculate hadroness for method of NEAREST NEIGHBORS
    //  - update the input files with the hadroness (ON1.root, OFF1.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 (OFF1.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
    //    and for the SUPERCUTS
    //  - update the input files with the hadronesses (ON1.root, OFF1.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, OFF1.root and MC1.root files with estimated energy
    //    (ON_XX1.root, OFF_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
  //=========
    // read ON data file 

    //  - produce the 2D-histogram "sigmabar versus Theta" 
    //    (SigmaTheta_ON.root) for ON data
    //    (to be used for the padding of the MC gamma data)

    //  - write a file of ON events (ON1.root) 
    //    (after the standard cuts, before the g/h separation)
    //    (to be used together with the corresponding MC gamma file (MC1.root)
    //     for the optimization of the g/h separation)


 if (JobA)
 {
    gLog << "=====================================================" << endl;
    gLog << "Macro CT1Analysis : Start of Job A" << endl;
    gLog << "" << endl;
    gLog << "Macro CT1Analysis : JobA, WPad, RPad, Wout = " 
         << JobA  << ",  " << WPad << ",  " << RPad << ",  " << Wout << endl;


    //--------------------------------------------------
    // names of ON and OFF files to be read
    // for generating the histograms to be used in the padding 
    TString fileON  = onfile;
    TString fileOFF = offfile;
    gLog << "fileON, fileOFF = " << fileON << ",  " << fileOFF << endl;

    // name of file to conatin the histograms for the padding
    TString outNameSigTh = outPath;
    outNameSigTh += "SigmaTheta";
    outNameSigTh += ".root";

    //--------------------------------------------------
    // type of data to be padded 
    TString typeInput = "ON";
    //TString typeInput = "OFF";
    //TString typeInput = "MC";
    gLog << "typeInput = " << typeInput << endl;


    // name of input root file
    if (typeInput == "ON")
      TString filenamein(onfile);
    else if (typeInput == "OFF")
      TString filenamein(offfile);
    else if (typeInput == "MC")
      TString filenamein(mcfile);
    gLog << "data to be padded : " << filenamein << endl;

    // name of output root file
    TString outNameImage = outPath;
    outNameImage += typeInput;
    outNameImage += "1.root";
    gLog << "padded data to be written onto : " << outNameImage << endl;

    //--------------------------------------------------

    //************************************************************
    // generate histograms to be used in the padding
    // 
    // read ON and OFF data files
    // generate (or read in) the padding histograms for ON and OFF data
    //                       and merge these histograms

    MCT1PadONOFF pad;
    pad.SetName("MCT1PadONOFF");
    pad.SetPadFlag(1);
    pad.SetDataType(typeInput);

    // generate the padding histograms
    if (!RPad)
    {
      gLog << "=====================================================" << endl;
      gLog << "Start generating the padding histograms" << endl;
      //-----------------------------------------
      // ON events

      gLog << "-----------" << endl;
      gLog << "ON events :" << endl;
      gLog << "-----------" << endl;

      MTaskList tliston;
      MParList pliston;

      MCT1ReadPreProc readON(fileON);

      MFCT1SelBasic selthetaon;
      selthetaon.SetCuts(-100.0, 29.5, 35.5);
      MContinue contthetaon(&selthetaon);

      MBlindPixelCalc blindon;
      blindon.SetUseBlindPixels();

      MFCT1SelBasic selbasicon;
      MContinue contbasicon(&selbasicon);

      MHBlindPixels blindON("BlindPixelsON");
      MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels");

      MSigmabarCalc sigbarcalcon;

      MHSigmaTheta sigthON("SigmaThetaON");
      MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt");
     
      //*****************************
      // entries in MParList
    
      pliston.AddToList(&tliston);
      InitBinnings(&pliston);
      pliston.AddToList(&blindON);
      pliston.AddToList(&sigthON);


      //*****************************
      // entries in MTaskList
    
      tliston.AddToList(&readON);
      //tliston.AddToList(&contthetaon);

      tliston.AddToList(&blindon);

      tliston.AddToList(&contbasicon);
      tliston.AddToList(&fillblindON);
      tliston.AddToList(&sigbarcalcon);
      tliston.AddToList(&fillsigthetaON);

      MProgressBar baron;
      MEvtLoop evtloopon;
      evtloopon.SetParList(&pliston);
      evtloopon.SetProgressBar(&baron);

      Int_t maxeventson = -1;
      //Int_t maxeventson = 20000;
      if ( !evtloopon.Eventloop(maxeventson) )
          return;

      tliston.PrintStatistics(0, kTRUE);

      blindON.DrawClone();
      sigthON.DrawClone();

      // save the histograms for the padding
      TH2D *sigthon     = sigthON.GetSigmaTheta();
      TH3D *sigpixthon  = sigthON.GetSigmaPixTheta();
      TH3D *diffpixthon = sigthON.GetDiffPixTheta();

      TH2D *blindidthon = blindON.GetBlindId();
      TH2D *blindnthon  = blindON.GetBlindN();

      //-----------------------------------------
      // OFF events

      gLog << "------------" << endl;
      gLog << "OFF events :" << endl;
      gLog << "------------" << endl;

      MTaskList tlistoff;
      MParList plistoff;

      MCT1ReadPreProc readOFF(fileOFF);

      MFCT1SelBasic selthetaoff;
      selthetaoff.SetCuts(-100.0, 29.5, 35.5);
      MContinue contthetaoff(&selthetaoff);

      MBlindPixelCalc blindoff;
      blindoff.SetUseBlindPixels();

      MFCT1SelBasic selbasicoff;
      MContinue contbasicoff(&selbasicoff);

      MHBlindPixels blindOFF("BlindPixelsOFF");
      MHBlindPixels *pblindOFF = &blindOFF;
      gLog << "pblindOFF = " << pblindOFF << endl;

      MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels");

      MSigmabarCalc sigbarcalcoff;

      MHSigmaTheta sigthOFF("SigmaThetaOFF");
      MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt");
     
      //*****************************
      // entries in MParList
    
      plistoff.AddToList(&tlistoff);
      InitBinnings(&plistoff);
      plistoff.AddToList(&blindOFF);
      plistoff.AddToList(&sigthOFF);


      //*****************************
      // entries in MTaskList
    
      tlistoff.AddToList(&readOFF);
      //tlistoff.AddToList(&contthetaoff);

      tlistoff.AddToList(&blindoff);

      tlistoff.AddToList(&contbasicoff);
      tlistoff.AddToList(&fillblindOFF);
      tlistoff.AddToList(&sigbarcalcoff);
      tlistoff.AddToList(&fillsigthetaOFF);

      MProgressBar baroff;
      MEvtLoop evtloopoff;
      evtloopoff.SetParList(&plistoff);
      evtloopoff.SetProgressBar(&baroff);

      Int_t maxeventsoff = -1;
      //Int_t maxeventsoff = 20000;
      if ( !evtloopoff.Eventloop(maxeventsoff) )
          return;

      tlistoff.PrintStatistics(0, kTRUE);

      blindOFF.DrawClone();
      sigthOFF.DrawClone();

      // save the histograms for the padding
      TH2D *sigthoff     = sigthOFF.GetSigmaTheta();
      TH3D *sigpixthoff  = sigthOFF.GetSigmaPixTheta();
      TH3D *diffpixthoff = sigthOFF.GetDiffPixTheta();

      TH2D *blindidthoff = blindOFF.GetBlindId();
      TH2D *blindnthoff  = blindOFF.GetBlindN();


      //-----------------------------------------

      gLog << "End of generating the padding histograms" << endl;
      gLog << "=====================================================" << endl;

      pad.MergeHistograms(sigthon,     sigthoff,
                          sigpixthon,  sigpixthoff,
                          diffpixthon, diffpixthoff,
                          blindidthon, blindidthoff,
                          blindnthon,  blindnthoff);

      if (WPad)
      {
        // write the padding histograms onto a file  ---------
        pad.WriteTargetDist(outNameSigTh);     
      }
    }

    // read the padding histograms ---------------------------
    if (RPad)
    {
      pad.ReadTargetDist(outNameSigTh);
    }


    //************************************************************

  if (Wout)
  {
    gLog << "=====================================================" << endl;
    gLog << "Start the padding" << endl;

    //-----------------------------------------------------------
    MTaskList tliston;
    MParList pliston;

    char *sourceName = "MSrcPosCam";
    MSrcPosCam source(sourceName);

    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");

    //-------------------------------------------
    // create the tasks which should be executed 
    //

    MCT1ReadPreProc read(filenamein);

    MFCT1SelBasic seltheta;
    seltheta.SetCuts(-100.0, 29.5, 35.5);
    MContinue conttheta(&seltheta);

    if (typeInput ==  "ON")
    {
      MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", 
                                                 "MCT1PointingCorrCalc");
    }

    MBlindPixelCalc blindbeforepad;
    blindbeforepad.SetUseBlindPixels();
    blindbeforepad.SetName("BlindBeforePadding");

    MBlindPixelCalc blind;
    blind.SetUseBlindPixels();
    blind.SetName("BlindAfterPadding");

    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");


    //*****************************
    // entries in MParList
    
    pliston.AddToList(&tliston);
    InitBinnings(&pliston);

    pliston.AddToList(&source);


    //*****************************
    // entries in MTaskList
    
    tliston.AddToList(&read);
    //tliston.AddToList(&conttheta);

    tliston.AddToList(&blindbeforepad);
    tliston.AddToList(&pad);
    if (typeInput ==  "ON")
      tliston.AddToList(&pointcorr);
    tliston.AddToList(&blind);

    tliston.AddToList(&contbasic);
    tliston.AddToList(&fillblind);
    tliston.AddToList(&sigbarcalc);
    tliston.AddToList(&fillsigtheta);
    tliston.AddToList(&clean);

    tliston.AddToList(&hcalc);
    tliston.AddToList(&hsrccalc);

    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    tliston.AddToList(&contstandard);
    tliston.AddToList(&write);

    //*****************************

    //-------------------------------------------
    // Execute event loop
    //
    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&pliston);
    evtloop.ReadEnv(env, "", printEnv);
    evtloop.SetProgressBar(&bar);
    //  evtloop.Write();

    Int_t maxevents = -1;
    //Int_t maxevents = 1000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // Display the histograms

    pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone();
    pliston.FindObject("BlindPixels", "MHBlindPixels")->DrawClone();

    pliston.FindObject("MHHillas")->DrawClone();
    pliston.FindObject("MHHillasExt")->DrawClone();
    pliston.FindObject("MHHillasSrc")->DrawClone();
    pliston.FindObject("MHNewImagePar")->DrawClone();
    pliston.FindObject("MHStarMap")->DrawClone();

    DeleteBinnings(&pliston);

      gLog << "End of padding" << endl;
      gLog << "=====================================================" << endl;
  }  


    gLog << "Macro CT1Analysis : End of Job A" << 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 (ON, OFF or MC)

    //TString typeInput = "ON";
    //TString typeInput = "OFF";
    TString typeInput = "MC";
    gLog << "typeInput = " << typeInput << endl;

    // name of input root file
    TString filenameData = outPath;
    filenameData += typeInput;
    filenameData += "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 filenameHad = outPath;
    filenameHad += "OFF";
    filenameHad += "1.root";
    Int_t howManyHadrons = 500000;
    gLog << "filenameHad = " << filenameHad << ",   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 = "OFF";
    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", filenameHad);
    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);

    //.......................................................................

      //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("HadRF",         "Events");
      write.AddContainer("HadSC",         "Events");
      write.AddContainer("HadNN",         "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");

    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(&fillhadnn);

    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("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 (ON, OFF or MC)

    TString typeInput = "ON";
    //TString typeInput = "OFF";
    //TString typeInput = "MC";
    gLog << "typeInput = " << typeInput << endl;

    // name of input root file
    TString filenameData = outPath;
    filenameData += typeInput;
    filenameData += "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 filenameHad = outPath;
    filenameHad += "OFF";
    filenameHad += "1.root";
    Int_t howManyHadrons = 1000000;
    gLog << "filenameHad = " << filenameHad << ",   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 = "OFF";
    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", filenameHad);
    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)
    {
        gLog << err << dbginf << "MRanForest not found... aborting." << endl;
        return kFALSE;
    }

    fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
    if (!fRanTree)                                  
    {                                                                          
        gLog << err << dbginf << "MRanTree not found... aborting." << endl;    
        return kFALSE;
    }

    }

    //-----------------------------------------------------------------
    // grow the trees of the random forest (event loop = tree loop)

    if (!RTree)
    {

    gLog << "" << endl;
    gLog << "========================================================" << endl;
    gLog << "Macro 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)
    {
        gLog << err << dbginf << "MRanForest not found... aborting." << endl;
        return kFALSE;
    }

    fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
    if (!fRanTree)                                  
    {                                                                          
        gLog << err << dbginf << "MRanTree not found... aborting." << endl;    
        return kFALSE;
    }

    }
    // end of growing the trees of the random forest
    //-----------------------------------------------------------------




    //-----------------------------------------------------------------
    // Update the 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);

    //.......................................................................
    // 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("HadRF",         "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(hadRFName);
    selfinalgh.SetName("SelFinalgh");
    MContinue contfinalgh(&selfinalgh);
    contfinalgh.SetName("ContSelFinalgh");

    MFillH fillranfor("MHRanForest");
    fillranfor.SetName("HRanForest");

    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
    fillhadrf.SetName("HhadRF");
    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
    fillhadsc.SetName("HhadSC");


    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(&sccalc);
    tliston.AddToList(&fillranfor);
    tliston.AddToList(&fillhadrf);
    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("MHRanForest")->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 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 += "OFF";
    filenameData += "2.root";
    gLog << "filenameData = " << filenameData << endl;

    // name of input MC file
    TString filenameMC = outPath;
    filenameMC += "MC";
    filenameMC += "2.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(hadRFName);
    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(hadRFName);
    selfinal.SetName("SelFinal");
    MContinue contfinal(&selfinal);
    contfinal.SetName("ContSelFinal");


    MFillH hfill1("MHHillas",    fHilName);
    hfill1.SetName("HHillas");

    MFillH hfill2("MHStarMap",   fHilName);
    hfill2.SetName("HStarMap");

    MFillH hfill3("MHHillasExt",    fHilNameSrc);
    hfill3.SetName("HHillasExt");
    
    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    hfill4.SetName("HHillasSrc");    

    MFillH hfill5("MHNewImagePar", fImgParName);
    hfill5.SetName("HNewImagePar");


    //*****************************
    // entries in MParList

    pliston.AddToList(&tliston);
    InitBinnings(&pliston);


    //*****************************
    // entries in MTaskList
    
    tliston.AddToList(&read);

    //tliston.AddToList(&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      = "2.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.30;
    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(&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();

    //-------------------------------------------
    // fit alpha distribution to get the number of excess events
    //

    MHHillasSrc* hillasSrc = 
      (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
    TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
  
    MHOnSubtraction onsub;
    onsub.Calc(alphaHist, &pliston, kTRUE, 13.1);
    //-------------------------------------------

    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 typeOFF = "OFF";
    TString typeMC  = "MC";
    TString ext    = "2.root";
    TString extout = "3.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 OFF file to be updated
    TString filenameOFF(outPath);
    filenameOFF += typeOFF;
    filenameOFF += ext;
    gLog << "filenameOFF = " << filenameOFF << 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 OFF file 
    TString filenameOFFup(outPath);
    filenameOFFup += typeOFF;
    filenameOFFup += "_";
    filenameOFFup += XX;
    filenameOFFup += extout;
    gLog << "filenameOFFup = " << filenameOFFup << 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; i<parA.GetSize(); i++)
      parA[i] = mcest.GetCoeff(i);
    for (Int_t i=0; i<parB.GetSize(); i++)
      parB[i] = mcest.GetCoeff( i+parA.GetSize() );


   if (WESTUP)
   {
    //==========   start update   ============================================
    //
    // Update ON, OFF and MC root files with the estimated energy

    //---------------------------------------------------
    // Update OFF data 
    //
    gLog << "============================================================"
         << endl;
    gLog << "Macro CT1Analysis.C : update file '" << filenameOFF
         << "'" << endl;

    MTaskList tlistoff;
    MParList plistoff;


    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)plistoff->FindCreateObj("MGeomCamCT1", "MGeomCam");

    //-------------------------------------------
    // create the tasks which should be executed 
    //

    MReadMarsFile read("Events", filenameOFF);
    read.DisableAutoScheme();

    //---------------------------
    // calculate estimated energy

    MEnergyEstParam eest2(fHilName);
    eest2.Add(fHilNameSrc);

    eest2.SetCoeffA(parA);
    eest2.SetCoeffB(parB);

    //.......................................................................

      MWriteRootFile write(filenameOFFup);

      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

    plistoff.AddToList(&tlistoff);
    InitBinnings(&plistoff);


    //*****************************
    // entries in MTaskList
    
    tlistoff.AddToList(&read);
    tlistoff.AddToList(&eest2);
    tlistoff.AddToList(&write);
    tlistoff.AddToList(&contfinal);

    //*****************************

    //-------------------------------------------
    // Execute event loop
    //
    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plistoff);
    evtloop.SetProgressBar(&bar);

    Int_t maxevents = -1;
    //Int_t maxevents = 1000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tlistoff.PrintStatistics(0, kTRUE);
    DeleteBinnings(&plistoff);

    //---------------------------------------------------

    //---------------------------------------------------
    // Update ON data
    //
    gLog << "============================================================"
         << endl;
    gLog << "Macro CT1Analysis.C : update file '" << filenameON
         << "'" << 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", 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      = "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;
    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();
    //MHMcCT1CollectionArea* collarea;

    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 = 
        (MHMcCT1CollectionArea*)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();

    //delete collarea;

    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;
 }
  //---------------------------------------------------------------------

}



