
#include "CT1EEst.C"

void InitBinnings(MParList *plist)
    {
      cout << "InitBinnings" << endl;

      MBinning *binssize = new MBinning("BinningSize");
      binssize->SetEdgesLog(50, 10, 1.0e5);
      plist->AddToList(binssize);

      MBinning *binsdistc = new MBinning("BinningDist");
      binsdistc->SetEdges(50, 0, 1.4);
      plist->AddToList(binsdistc);

      MBinning *binswidth = new MBinning("BinningWidth");
      binswidth->SetEdges(50, 0, 1.0);
      plist->AddToList(binswidth);

      MBinning *binslength = new MBinning("BinningLength");
      binslength->SetEdges(50, 0, 1.0);
      plist->AddToList(binslength);

      MBinning *binsalpha = new MBinning("BinningAlpha");
      binsalpha->SetEdges(100, -100, 100);
      plist->AddToList(binsalpha);

      MBinning *binsasym = new MBinning("BinningAsym");
      binsasym->SetEdges(50, -1.5, 1.5);
      plist->AddToList(binsasym);

      MBinning *binsht = new MBinning("BinningHeadTail");
      binsht->SetEdges(50, -1.5, 1.5);
      plist->AddToList(binsht);

      MBinning *binsm3l = new MBinning("BinningM3Long");
      binsm3l->SetEdges(50, -1.5, 1.5);
      plist->AddToList(binsm3l);

      MBinning *binsm3t = new MBinning("BinningM3Trans");
      binsm3t->SetEdges(50, -1.5, 1.5);
      plist->AddToList(binsm3t);

   
      //.....
      MBinning *binsb = new MBinning("BinningSigmabar");
      binsb->SetEdges( 100,  0.0,  5.0);
      plist->AddToList(binsb);

      Int_t nbins = 8;
      TArrayD xbins(nbins+1);
      xbins[0] =  0.0;
      xbins[1] =  7.5;
      xbins[2] = 17.5; 
      xbins[3] = 23.5;
      xbins[4] = 29.5;
      xbins[5] = 35.5;
      xbins[6] = 42.0;
      xbins[7] = 55.0;
      xbins[8] = 68.0;

      MBinning *binth = new MBinning("BinningTheta");
      binth->SetEdges(xbins);
      plist->AddToList(binth);

      MBinning *binsdiff = new MBinning("BinningDiffsigma2");
      binsdiff->SetEdges(100, -5.0, 20.0);
      plist->AddToList(binsdiff);
    }

void DeleteBinnings(MParList *plist)
    {
      cout << "DeleteBinnings" << endl;

      TObject *bin;

      bin = plist->FindObject("BinningSize");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningDist");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningWidth");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningLength");
      if (!bin) delete bin;

      bin = plist->FindObject("BinningAlpha");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningAsym");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningHeadTail");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningM3Long");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningM3Trans");
      if (!bin) delete bin;

      //.....
      bin = plist->FindObject("BinningSigmabar");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningTheta");
      if (!bin) delete bin;
      bin = plist->FindObject("BinningDiffsigma2");
      if (!bin) delete bin;
    }

void CT1Analysis()
{
    //-----------------------------------------------
    const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; 

    //const char *onfile  = "~magican/ct1test/wittek/mkn421_on.preproc"; 
    const char *onfile  = "~magican/ct1test/wittek/mkn421_00-01"; 

    const char *mcfile  = "~magican/ct1test/wittek/mkn421_mc_pedrms_0.001.preproc";
    //-----------------------------------------------

    // path for input for Mars
    TString inPath = "~magican/ct1test/wittek/marsoutput/optionB/";

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



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

    // Job A_ON : read ON data
    //  - generate sigmabar vs. Theta plot; 
    //  - write root file for ON data (ON1.root);
    //  - generate hadron matrix for g/h separation

    Bool_t JobA_ON = kFALSE;  
    Bool_t WHistON = kFALSE;   // write out histogram sigmabar vs. Theta ?
    Bool_t WLookON = kFALSE;   // write out look-alike events for hadrons ?
    Bool_t WON1    = kFALSE;   // write out root file ON1.root ?


    // Job A_MC : read MC gamma data, 
    //  - read sigmabar vs. Theta plot from ON data  
    //  - do padding; 
    //  - write root file for MC gammas (MC1.root);
    //  - generate gamma matrix for g/h separation)

    Bool_t JobA_MC  = kFALSE;  
    Bool_t WLookMC  = kFALSE;  // write out look-alike events for gammas ?
    Bool_t WMC1     = kFALSE;  // write out root file MC1.root ?


    // Job B_NN_UP : read ON1.root (or MC1.root) file 
    //  - read hadron and gamma matrix
    //  - calculate hadroness for method of Nearest Neighbors
    //    and for the supercuts
    //  - update the input files with the hadronesses (ON2.root or MC2.root)

    Bool_t JobB_NN_UP  = kFALSE;  
    Bool_t WNN         = kFALSE;  // write root file ?


    // Job B_EST_UP : 
    //  - read MC2.root file 
    //  - select g/h separation method XX
    //  - optimize energy estimation for events passing the final cuts
    //  - write parameters of energy estimator onto file
    //  - update ON2.root and MC2.root files with estimated energy
    //    (ON_XX2.root and MC_XX2.root)

    Bool_t JobB_EST_UP  = kFALSE;  
    Bool_t WESTUP       = kFALSE;  // update root file ?


    // Job C_NN : read ON1.root and MC1.root files  
    //  - read look-alike events for hadrons and gammas 
    //  - calculate hadronness for method of NEAREST NEIGHBOR
    //  - produce Neyman-Pearson plot

    Bool_t JobC_NN  = kFALSE;  


    // Job C_RF : read ON1.root and MC1.root files  
    //  - read look-alike events for hadrons and gammas 
    //  - calculate hadronness for method of RANDOM FOREST
    //  - produce Neyman-Pearson plot

    Bool_t JobC_RF  = kFALSE;  


    // Job D_XX :  
    //  - select g/h separation method XX
    //  - read MC_XX2.root file 
    //  - calculate eff. collection area
    //  - read ON_XX2.root file 
    //  - apply final cuts
    //  - calculate flux
    //  - write root file for ON data after final cuts (ON3.root))
    Bool_t JobD_XX  = kTRUE;  
    Bool_t WXXD     = kFALSE;  // write out root file ON3.root ?


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

    
  //---------------------------------------------------------------------
  // Job A_ON
  //=========
    // read ON data file 

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

    //  - write a file of look-alike events (matrix_hadrons_ON.root)
    //    (to be used later together with the corresponding MC gamma file
    //     for the g/h separation in the analysis of the ON data)

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


 if (JobA_ON)
 {
    cout << "=====================================================" << endl;
    cout << "Macro CT1Analysis : Start of Job A_ON" << endl;
    cout << "" << endl;
    cout << "Macro CT1Analysis : JobA_ON, WHistON, WLookON, WON1 = " 
         << JobA_ON  << ",  " << WHistON << ",  " << WLookON << ",  " 
         << WON1 << endl;

    TString typeInput = "ON";
    cout << "typeInput = " << typeInput << endl;

    // name of input root file
    TString filenamein(onfile);

    // name of output root file
    TString outNameImage = outPath;
    outNameImage += typeInput;
    outNameImage += "1.root";


    //-----------------------------------------------------------
    MTaskList tliston;
    MParList pliston;
    MFilterList flist;

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


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

    MCT1ReadPreProc read(filenamein);

    MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", 
                                               "MCT1PointingCorrCalc");

    MBlindPixelCalc blind;

    MFCT1SelBasic selbasic;
    MContinue contbasic(&selbasic);

    MSigmabarCalc sigbarcalc;

    MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
    fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");

    MImgCleanStd    clean; 

    // calculate also extended image parameters
    TString fHilName = "MHillas";
    MHillasExt hext(fHilName);

    MHillasCalc    hcalc(fHilName);

    TString fHilSrcName = "MHillasSrc";
    MHillasSrc hilsrc(fHilSrcName);

    MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", "");
    csrc1.SetInput(fHilName);

    TString fnewparName = "MNewImagePar";
    MNewImagePar newpar(fnewparName);

    MNewImageParCalc newimage(sourceName, fnewparName);
    newimage.SetInput(fHilName);

    MFCT1SelStandard selstandard(fHilName, fHilSrcName);
    MContinue contstandard(&selstandard);

    MFillH hfill1("MHHillas",    fHilName);
    hfill1.SetTitle("Task to plot the source independent Hillas parameters");

    MFillH hfill2("MHStarMap",   fHilName);
    hfill2.SetTitle("Task to plot the star map (points along main axis of shower)");

    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);

    MFillH hfill3("MHHillasExt",   fHilSrcName);
    hfill3.SetTitle("Task to plot the extended Hillas parameters");

    MFillH hfill4("MHHillasSrc",   fHilSrcName);
    hfill4.SetTitle("Task to plot the source dependent Hillas parameters");

    MFillH hfill5("MHNewImagePar", fnewparName);
    hfill5.SetTitle("Task to plot the new image parameters");

    //+++++++++++++++++++++++++++++++++++++++++++++++++++
    // prepare writing of look-alike events (for hadrons)

    const char* mtxName = "MatrixHadrons";
    Int_t howMany = 500000;
    TString outName = outPath;
    outName += "matrix_hadrons_";
    outName += typeInput;
    outName += ".root";


    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Matrix for (hadrons)" << endl;
    cout << "input file                    = " << filenamein<< endl;
    cout << "matrix name                   = " << mtxName << endl;
    cout << "max. no.of look-alike events  = " << howMany << endl;
    cout << "name of output root file      = " << outName << endl;
    cout << "" << endl;


    // matrix limitation for look-alike events (approximate number)
    MFEventSelector selector("", "");
    selector.SetNumSelectEvts(howMany);

    //
    // --- setup the matrix and add it to the parlist
    //
    MHMatrix *pmatrix = new MHMatrix(mtxName);
    MHMatrix& matrix = *pmatrix;

    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
    matrix.AddColumn("MSigmabar.fSigmabar");
    matrix.AddColumn("log10(MHillas.fSize)");
    matrix.AddColumn("MHillasSrc.fDist");
    matrix.AddColumn("MHillas.fWidth");
    matrix.AddColumn("MHillas.fLength");
    matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
    matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)");
    matrix.AddColumn("MHillas.fConc");
    matrix.AddColumn("MNewImagePar.fLeakage1");

    matrix.Print("cols");

    MFillH fillmatg(mtxName);
    fillmatg.SetFilter(&flist);
    //+++++++++++++++++++++++++++++++++++++++++++++++++++

    if (WON1)
    {    
      MWriteRootFile write(outNameImage);

      write.AddContainer("MRawRunHeader", "RunHeaders");
      write.AddContainer("MTime",         "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("MSrcPosCam",    "Events");
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("MHillas",       "Events");
      write.AddContainer("MHillasSrc",    "Events");
      write.AddContainer("MNewImagePar",  "Events");
    }

    //*****************************
    // entries in MFilterList

    flist.AddToList(&selector);


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

    pliston.AddToList(&source);
    pliston.AddToList(&sigbarcalc);

    pliston.AddToList(&hext);
    pliston.AddToList(&hilsrc);
    pliston.AddToList(&newpar);
    pliston.AddToList(&hhilext);

    pliston.AddToList(pmatrix);


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

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

    tliston.AddToList(&hcalc);
    tliston.AddToList(&csrc1);
    tliston.AddToList(&newimage);

    tliston.AddToList(&flist);
    tliston.AddToList(&fillmatg);    

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

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

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

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

    //Int_t maxevents = 1000000000;
    Int_t maxevents = 20000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);
    DeleteBinnings(&pliston);

    //-------------------------------------------
    // Display the histograms
    //
    TObject *c;
    c = ((MHSigmaTheta*)pliston.FindObject("SigmaTheta", "MHSigmaTheta"))
                                                            ->DrawClone();

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


    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    //
    // ----------------------------------------------------------
    //  Definition of the reference sample of look-alike events
    //  (this is a subsample of the original sample)
    // ----------------------------------------------------------
    //
    cout << "" << endl;
    cout << "========================================================" << endl;
    // Select a maximum of nmaxevts events from the sample of look-alike 
    // events. They will form the reference sample.
    Int_t nmaxevts  = 2000;

    // target distribution for the variable in column refcolumn (start at 0);
    Int_t   refcolumn = 0;
    Int_t   nbins   =   5;
    Float_t frombin = 0.5;
    Float_t tobin   = 1.0;
    TH1F *thsh = new TH1F("thsh","target distribution", 
                           nbins, frombin, tobin);
    thsh->SetXTitle("cos( \\Theta)");
    thsh->SetYTitle("Counts");
    Float_t dbin = (tobin-frombin)/nbins;
    Float_t lbin = frombin +dbin*0.5;
    for (Int_t j=1; j<=nbins; j++) 
    {
      thsh->Fill(lbin,1.0);
      lbin += dbin;
    }

    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Macro CT1Analysis : define reference sample for hadrons" << endl; 
    cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
         << refcolumn << ",  " << nmaxevts << endl;    

    matrix.EnableGraphicalOutput();
    Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts);

    if ( !matrixok ) 
    {
      cout << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
      return;
    }
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    //-------------------------------------------
    // write out look-alike events (for hadrons)
    //
    if (WLookON)
    {
      TFile *writejens = new TFile(outName, "RECREATE", "");
      matrix.Write();

      cout << "" << endl;
      cout << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
           << outName << endl;

      writejens->Close();
      delete writejens;
    }


    //-------------------------------------------
    // Write histograms onto a file
  if (WHistON)
  {
      MHSigmaTheta *sigtheta = 
            (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
      if (!sigtheta)
	{
          cout << "Object 'SigmaTheta' not found" << endl;
          return;
	}
      TH2D *fHSigTh    = sigtheta->GetSigmaTheta();
      TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta();
      TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();

      TString outNameSigTh = outPath;
      outNameSigTh += "SigmaTheta_";
      outNameSigTh += typeInput;
      outNameSigTh += ".root";

      TFile outfile(outNameSigTh, "RECREATE");
      fHSigTh->Write();
      fHSigPixTh->Write();
      fHDifPixTh->Write();
      outfile.Close();
     
      cout << "" << endl;
      cout << "File " << outNameSigTh << " was written out" << endl;
  }


    cout << "Macro CT1Analysis : End of Job A_ON" << endl;
    cout << "===================================================" << endl;
 }


  //---------------------------------------------------------------------
   // Job A_MC
   //=========

    // read MC gamma data  

    //    - to pad them
    //      (using the 2D-histogram "sigmabar versus Theta" 
    //       (SigmaTheta_ON.root)  of the ON data)

    //    - to write a file of look-alike events (matrix_gammas_MC.root)
    //      (to be used later together with the corresponding hadron file
    //       for the g/h separation in the analysis of the ON data)

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


 if (JobA_MC)
 {
    cout << "=====================================================" << endl;
    cout << "Macro CT1Analysis : Start of Job A_MC" << endl;

    cout << "" << endl;
    cout << "Macro CT1Analysis : JobA_MC, WLookMC, WMC1 = " 
         << JobA_MC  << ",  " << WLookMC << ",  " << WMC1 << endl;

    TString typeInput = "MC";
    cout << "typeInput = " << typeInput << endl;


    // name of input root file
    TString filenamein(mcfile);

    // name of output root file
    TString outNameImage = outPath;
    outNameImage += typeInput;
    outNameImage += "1.root";


    // use for padding sigmabar vs. Theta from ON data
    TString typeHist = "ON";
    cout << "typeHist = " << typeHist << endl;

    //------------------------------------
    // Get the histograms "2D-ThetaSigmabar"
    // and                "3D-ThetaPixSigma"
    // and                "3D-ThetaPixDiff"


      TString outNameSigTh = outPath;
      outNameSigTh += "SigmaTheta_";
      outNameSigTh += typeHist;
      outNameSigTh += ".root";


      cout << "Reading in file " << outNameSigTh << endl;

      TFile *infile = new TFile(outNameSigTh);
      infile->ls();

      TH2D *fHSigmaTheta = 
      (TH2D*) gROOT.FindObject("2D-ThetaSigmabar");
      if (!fHSigmaTheta)
	{
          cout << "Object '2D-ThetaSigmabar' not found on root file" << endl;
          return;
	}
      cout << "Object '2D-ThetaSigmabar' was read in" << endl;

      TH3D *fHSigmaPixTheta = 
      (TH3D*) gROOT.FindObject("3D-ThetaPixSigma");
      if (!fHSigmaPixTheta)
	{
          cout << "Object '3D-ThetaPixSigma' not found on root file" << endl;
          return;
	}
      cout << "Object '3D-ThetaPixSigma' was read in" << endl;

      TH3D *fHDiffPixTheta = 
      (TH3D*) gROOT.FindObject("3D-ThetaPixDiff");
      if (!fHSigmaPixTheta)
	{
          cout << "Object '3D-ThetaPixDiff' not found on root file" << endl;
          return;
	}
      cout << "Object '3D-ThetaPixDiff' was read in" << endl;

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

    MTaskList tlist;
    MParList plist;
    MFilterList flist;

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

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

    MCT1ReadPreProc read(filenamein);

    MBlindPixelCalc blind;

    MFCT1SelBasic selbasic;
    MContinue contbasic(&selbasic);


    // There are 2 options for Thomas Schweizer's padding
    //     fPadFlag = 1   get Sigmabar from fHSigmaTheta
    //                    and Sigma    from fHDiffPixTheta
    //     fPadFlag = 2   get Sigma    from fHSigmaPixTheta
    
    MPadSchweizer padthomas("MPadSchweizer","Task for the padding (Schweizer)",
                             fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta);
    padthomas.SetPadFlag(1);

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

    MSigmabarCalc sigbarcalc;

    MFillH fillsigtheta ("MCSigmaTheta[MHSigmaTheta]", "MMcEvt");
    fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");

    MImgCleanStd    clean; 

    // calculate also extended image parameters
    TString fHilName = "MHillas";
    MHillasExt hext(fHilName);

    MHillasCalc    hcalc(fHilName);

    TString fHilSrcName = "MHillasSrc";
    MHillasSrc hilsrc(fHilSrcName);

    MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", "");
    csrc1.SetInput(fHilName);

    TString fnewparName = "MNewImagePar";
    MNewImagePar newpar(fnewparName);

    MNewImageParCalc newimage(sourceName, fnewparName);
    newimage.SetInput(fHilName);

    MFCT1SelStandard selstandard(fHilName, fHilSrcName);
    MContinue contstandard(&selstandard);

    MFillH hfill1("MHHillas",    fHilName);
    hfill1.SetTitle("Task to plot the source independent Hillas parameters");

    MFillH hfill2("MHStarMap",   fHilName);
    hfill2.SetTitle("Task to plot the star map (points along main axis of shower)");

    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);

    MFillH hfill3("MHHillasExt",   fHilSrcName);
    hfill3.SetTitle("Task to plot the extended Hillas parameters");

    MFillH hfill4("MHHillasSrc",   fHilSrcName);
    hfill4.SetTitle("Task to plot the source dependent Hillas parameters");

    MFillH hfill5("MHNewImagePar", fnewparName);
    hfill5.SetTitle("Task to plot the new image parameters");


    //+++++++++++++++++++++++++++++++++++++++++++++++++++
    // prepare writing of look-alike events (for gammas)

    const char* mtxName = "MatrixGammas";
    Int_t howMany = 500000;

    TString outName = outPath;
    outName += "matrix_gammas_";
    outName += typeInput;
    outName += ".root";


    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Matrix for (gammas)" << endl;
    cout << "input file                    = " << filenamein<< endl;
    cout << "matrix name                   = " << mtxName << endl;
    cout << "max. no.of look-alike events  = " << howMany << endl;
    cout << "name of output root file      = " << outName << endl;
    cout << "" << endl;


    // matrix limitation for look-alike events (approximate number)
    MFEventSelector selector("", "");
    selector.SetNumSelectEvts(howMany);

    //
    // --- setup the matrix and add it to the parlist
    //
    MHMatrix *pmatrix = new MHMatrix(mtxName);
    MHMatrix& matrix = *pmatrix;

    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
    matrix.AddColumn("MSigmabar.fSigmabar");
    matrix.AddColumn("log10(MHillas.fSize)");
    matrix.AddColumn("MHillasSrc.fDist");
    matrix.AddColumn("MHillas.fWidth");
    matrix.AddColumn("MHillas.fLength");
    matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
    matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)");
    matrix.AddColumn("MHillas.fConc");
    matrix.AddColumn("MNewImagePar.fLeakage1");

    matrix.Print("cols");

    MFillH fillmatg(mtxName);
    fillmatg.SetFilter(&flist);
    //+++++++++++++++++++++++++++++++++++++++++++++++++++


    if (WMC1)
    {    
      MWriteRootFile write(outNameImage);

      write.AddContainer("MRawRunHeader", "RunHeaders");
      write.AddContainer("MTime",         "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("MSrcPosCam",    "Events");
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("MHillas",       "Events");
      write.AddContainer("MHillasSrc",    "Events");
      write.AddContainer("MNewImagePar",  "Events");
    }

    //*****************************
    // entries in MFilterList

    flist.AddToList(&selector);


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

    plist.AddToList(&tlist);
    InitBinnings(&plist);

    plist.AddToList(&source);
    plist.AddToList(&sigbarcalc);

    plist.AddToList(&hext);
    plist.AddToList(&hilsrc);
    plist.AddToList(&newpar);
    plist.AddToList(&hhilext);

    plist.AddToList(pmatrix); 


    //*****************************
    // entries in MTaskList

    tlist.AddToList(&read);
    tlist.AddToList(&blind);
    tlist.AddToList(&padthomas);

    tlist.AddToList(&contbasic);
    tlist.AddToList(&sigbarcalc);
    tlist.AddToList(&fillsigtheta);
    tlist.AddToList(&clean);

    tlist.AddToList(&hcalc);
    tlist.AddToList(&csrc1);
    tlist.AddToList(&newimage);

    tlist.AddToList(&flist);
    tlist.AddToList(&fillmatg);

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

    tlist.AddToList(&contstandard);
    if (WMC1)
      tlist.AddToList(&write);

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


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

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

    tlist.PrintStatistics(0, kTRUE);
    DeleteBinnings(&plist);

    //-------------------------------------------
    // Display the histograms
    //
    TObject *c;
    c = ((MHSigmaTheta*)plist.FindObject("MCSigmaTheta", "MHSigmaTheta"))
                                                            ->DrawClone();
    c = plist.FindObject("MHHillas")->DrawClone();
    c = plist.FindObject("MHHillasExt")->DrawClone();
    c = plist.FindObject("MHHillasSrc")->DrawClone();
    c = plist.FindObject("MHNewImagePar")->DrawClone();
    c = plist.FindObject("MHStarMap")->DrawClone();


    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    //
    // ----------------------------------------------------------
    //  Definition of the reference sample of look-alike events
    //  (this is a subsample of the original sample)
    // ----------------------------------------------------------
    //
    cout << "" << endl;
    cout << "========================================================" << endl;
    // Select a maximum of nmaxevts events from the sample of look-alike 
    // events. They will form the reference sample.
    Int_t nmaxevts  = 2000;

    // target distribution for the variable in column refcolumn (start at 0);
    Int_t   refcolumn = 0;
    Int_t   nbins   =   5;
    Float_t frombin = 0.5;
    Float_t tobin   = 1.0;
    TH1F *thsh = new TH1F("thsh","target distribution", 
                           nbins, frombin, tobin);
    Float_t dbin = (tobin-frombin)/nbins;
    Float_t lbin = frombin +dbin*0.5;
    for (Int_t j=1; j<=nbins; j++) 
    {
      thsh->Fill(lbin,1.0);
      lbin += dbin;
    }

    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Macro CT1Analysis : define reference sample for gammas" << endl; 
    cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
         << refcolumn << ",  " << nmaxevts << endl;    

    matrix.EnableGraphicalOutput();
    Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts);
    if ( !matrixok ) 
    {
      cout << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
      return;
    }
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    //-------------------------------------------
    // write out look-alike events (for gammas)
    //
    if (WLookMC)
    {
      TFile *writejens = new TFile(outName, "RECREATE", "");
      matrix.Write();
      cout << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
           << outName << endl;

      writejens->Close();
      delete writejens;
    }

    cout << "Macro CT1Analysis : End of Job A_MC" 
         << endl;
    cout << "=========================================================" 
         << endl;
 }



  //---------------------------------------------------------------------
  // Job B_NN_UP
  //============

    //  - read in the matrices of look-alike events for gammas and hadrons

    // then read ON1.root (or MC1.root) file 
    //
    //  - calculate the hadroness for the method of nearest neighbors
    //
    //  - update input root file, including the hadroness


 if (JobB_NN_UP)
 {
    cout << "=====================================================" << endl;
    cout << "Macro CT1Analysis : Start of Job B_NN_UP" << endl;

    cout << "" << endl;
    cout << "Macro CT1Analysis : JobB_NN_UP, WNN = " 
         << JobB_NN_UP  << ",  " << WNN << endl;


    //TString typeInput = "ON";
    TString typeInput = "MC";
    cout << "typeInput = " << typeInput << endl;

    // name of input root file
    TString filenameData = outPath;
    filenameData += typeInput;
    filenameData += "1.root";
    cout << "filenameData = " << filenameData << endl; 

    // name of output root file
    TString outNameImage = outPath;
    outNameImage += typeInput;
    outNameImage += "2.root";

    //$$$$$$$$$$$$$$$$$$$$$$$$$$
    outNameImage = filenameData;
    //$$$$$$$$$$$$$$$$$$$$$$$$$$

    cout << "outNameImage = " << outNameImage << endl; 

    TString typeMatrixHadrons = "ON";
    cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;

   //*************************************************************************
   // read in matrices of look-alike events

    const char* mtxName = "MatrixGammas";

    TString outName = outPath;
    outName += "matrix_gammas_";
    outName += "MC";
    outName += ".root";


    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Get matrix for (gammas)" << endl;
    cout << "matrix name        = " << mtxName << endl;
    cout << "name of root file  = " << outName << endl;
    cout << "" << endl;


    // read in the object with the name 'mtxName' from file 'outName'
    //
    TFile *fileg = new TFile(outName); 

    MHMatrix matrixg;
    matrixg.Read(mtxName);
    matrixg.Print("cols");

    delete fileg;

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

    const char* mtxName = "MatrixHadrons";

    TString outName = outPath;
    outName += "matrix_hadrons_";
    outName += typeMatrixHadrons;
    outName += ".root";


    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << " Get matrix for (hadrons)" << endl;
    cout << "matrix name        = " << mtxName << endl;
    cout << "name of root file  = " << outName << endl;
    cout << "" << endl;


    // read in the object with the name mtxName
    //
    TFile *fileh = new TFile(outName); 

    MHMatrix matrixh;
    matrixh.Read(mtxName);
    matrixh.Print("cols");

    delete fileh;

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


    MTaskList tliston;
    MParList pliston;

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

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

    TString fHilName    = "MHillas"; 
    TString fHilSrcName = "MHillasSrc"; 
    TString fnewparName = "MNewImagePar"; 


    //.......................................................................
    // calculation of hadroness for method of Nearest Neighbors

    TString hadNNName = "HadNN";
    MMultiDimDistCalc nncalc;
    nncalc.SetUseNumRows(25);
    nncalc.SetUseKernelMethod(kFALSE);
    nncalc.SetHadronnessName(hadNNName);

    //.......................................................................
    // calculation of hadroness for the supercuts
    // (=0.25 if fullfilled, =0.75 otherwise)

    TString hadSCName = "HadSC";
    MCT1SupercutsCalc sccalc(fHilName, fHilSrcName);
    sccalc.SetHadronnessName(hadSCName);

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

    if (WNN)
    {    
      //MWriteRootFile write(outNameImage);
      MWriteRootFile write(outNameImage, "UPDATE");

      //write.AddContainer("MRawRunHeader", "RunHeaders");
      //write.AddContainer("MTime",         "Events");
      //write.AddContainer("MMcEvt",        "Events");
      //write.AddContainer("MSrcPosCam",    "Events");
      //write.AddContainer("MSigmabar",     "Events");
      //write.AddContainer("MHillas",       "Events");
      //write.AddContainer("MHillasSrc",    "Events");
      //write.AddContainer("MNewImagePar",  "Events");

      write.AddContainer(hadNNName,       "Events");
      write.AddContainer(hadSCName,       "Events");
    }

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


    Float_t hadcut   =  0.40;
    Float_t alphacut = 100.0;
    MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
    selfinalgh.SetCuts(hadcut, alphacut);
    selfinalgh.SetHadronnessName(hadNNName);
    MContinue contfinalgh(&selfinalgh);

    MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);

    Float_t hadcut   = 0.40;
    Float_t alphacut = 20.0;
    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetCuts(hadcut, alphacut);
    selfinal.SetHadronnessName(hadNNName);
    MContinue contfinal(&selfinal);


    MFillH hfill1("MHHillas",    fHilName);
    MFillH hfill2("MHStarMap",   fHilName);

    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);

    MFillH hfill3("MHHillasExt",   fHilSrcName);
    MFillH hfill4("MHHillasSrc",   fHilSrcName);
    MFillH hfill5("MHNewImagePar", fnewparName);


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

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

    pliston.AddToList(&matrixg);
    pliston.AddToList(&matrixh);

    pliston.AddToList(&hhilext); 

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

    tliston.AddToList(&nncalc);
    tliston.AddToList(&sccalc);
    tliston.AddToList(&fillhadnn);
    tliston.AddToList(&fillhadsc);

    if (WNN)
      tliston.AddToList(&write);

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

    tliston.AddToList(&contfinal);

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

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

    if ( !evtloop.Eventloop() )
      //if ( !evtloop.Eventloop(1000) )
        return;

    tliston.PrintStatistics(0, kTRUE);
    DeleteBinnings(&pliston);

    //-------------------------------------------
    // Display the histograms
    //
    TObject *c;

    c = pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
    c = pliston.FindObject("hadSC", "MHHadronness")->DrawClone();

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


    cout << "Macro CT1Analysis : End of Job B_NN_UP" << endl;
    cout << "=======================================================" << endl;
 }
  //---------------------------------------------------------------------



  //---------------------------------------------------------------------
  // Job B_EST_UP
  //============

    //  - read MC2.root file 
    //  - select g/h separation method XX
    //  - optimize energy estimator for events passing final cuts
    //  - write parameters of energy estimator onto file "energyest_XX.root"
    //
    //  - read ON2.root and MC2.root files 
    //  - update input root file with the estimated energy
    //    (ON_XX2.root, MC_XX2.root)


 if (JobB_EST_UP)
 {
    cout << "=====================================================" << endl;
    cout << "Macro CT1Analysis : Start of Job B_EST_UP" << endl;

    cout << "" << endl;
    cout << "Macro CT1Analysis : JobB_EST_UP, WESTUP = " 
         << JobB_EST_UP  << ",  " << WESTUP << endl;


    TString typeON = "ON";
    TString typeMC = "MC";
    TString ext    = "3.root";

    //------------------------------
    // name of MC file to be used for optimizing the energy estimator
    TString filenameOpt(outPath);
    filenameOpt += typeMC;
    filenameOpt += ext; 
    cout << "filenameOpt = " << filenameOpt << endl;

    //------------------------------
    // selection of g/h separation method
    // and definition of final selections

    TString XX("NN");
    //TString XX("SC");
    //TString XX("RF");
    TString fhadronnessName("Had");
    fhadronnessName += XX;
    cout << "fhadronnessName = " << fhadronnessName << endl;

    // maximum values of the hadronness and of |alpha|
    Float_t maxhadronness   = 0.40;
    Float_t maxalpha        = 20.0;
    cout << "Maximum values of hadronness and |alpha| = "
         << maxhadronness << ",  " << maxalpha << endl;

    // name of file containing the parameters of the energy estimator
    TString energyParName(outPath);
    energyParName += "energyest_";
    energyParName += XX;
    energyParName += ".root";
    cout << "energyParName = " << energyParName << endl;


    //------------------------------
    // name of ON file to be updated
    TString filenameON(outPath);
    filenameON += typeON;
    filenameON += ext;
    cout << "filenameON = " << filenameON << endl;

    // name of MC file to be updated
    TString filenameMC(outPath);
    filenameMC += typeMC;
    filenameMC += ext;
    cout << "filenameMC = " << filenameMC << endl;

    //------------------------------
    // name of updated ON file 
    TString filenameONup(outPath);
    filenameONup += typeON;
    filenameONup += "_";
    filenameONup += XX;
    filenameONup += ext;
    cout << "filenameONup = " << filenameONup << endl;

    // name of updated MC file 
    TString filenameMCup(outPath);
    filenameMCup += typeMC;
    filenameMCup += "_";
    filenameMCup += XX;
    filenameMCup += ext;
    cout << "filenameMCup = " << filenameMCup << endl;

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

    TString fHilName    = "MHillas"; 
    TString fHilSrcName = "MHillasSrc"; 
    TString fnewparName = "MNewImagePar"; 

    //===========================================================
    //
    // Optimization of energy estimator
    //

    TString inpath("");
    TString outpath("");
    Int_t howMany = 2000;
    CT1EEst(inpath,   filenameOpt,   outpath, energyParName, 
            fHilName, fHilSrcName,   fhadronnessName,
            howMany,  maxhadronness, maxalpha);

    //-----------------------------------------------------------
    //
    // Read in parameters of energy estimator
    //

    TFile enparam(energyParName);
 
    TArrayD parA(5);
    TArrayD parB(7);
    for (Int_t i=0; i<parA.GetSize(); i++)
      parA[i] = EnergyParams(i);
    for (Int_t i=0; i<parB.GetSize(); i++)
      parB[i] = EnergyParams( i+parA.GetSize() );


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

    //---------------------------------------------------
    // Update ON data
    //
    MTaskList tliston;
    MParList pliston;

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

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

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

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

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

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

      MWriteRootFile write(filenameONup);

      write.AddContainer("MRawRunHeader", "RunHeaders");
      write.AddContainer("MTime",         "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("MSrcPosCam",    "Events");
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("MHillas",       "Events");
      write.AddContainer("MHillasSrc",    "Events");
      write.AddContainer("MNewImagePar",  "Events");

      write.AddContainer("HadNN",         "Events");
      write.AddContainer("HadSC",         "Events");

      write.AddContainer("MEnergyEst",    "Events");

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

    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetCuts(maxhadronness, maxalpha);
    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);
    //evtloop.Write();

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

    tliston.PrintStatistics(0, kTRUE);
    DeleteBinnings(&pliston);

    //---------------------------------------------------
    //---------------------------------------------------
    // Update MC data
    //
    MTaskList tliston;
    MParList pliston;

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

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

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

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

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

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

      MWriteRootFile write(filenameMCup);

      write.AddContainer("MRawRunHeader", "RunHeaders");
      write.AddContainer("MTime",         "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("MSrcPosCam",    "Events");
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("MHillas",       "Events");
      write.AddContainer("MHillasSrc",    "Events");
      write.AddContainer("MNewImagePar",  "Events");

      write.AddContainer("HadNN",         "Events");
      write.AddContainer("HadSC",         "Events");

      write.AddContainer("MEnergyEst",    "Events");

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

    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetCuts(maxhadronness, maxalpha);
    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);
    //evtloop.Write();

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

    tliston.PrintStatistics(0, kTRUE);
    DeleteBinnings(&pliston);


    //==========   end update   ============================================
   }
    
    enparam.Close();

    cout << "Macro CT1Analysis : End of Job B_EST_UP" << endl;
    cout << "=======================================================" << endl;
 }
  //---------------------------------------------------------------------


  //---------------------------------------------------------------------
  // Job C_NN
  //=========

    // read ON1 and MC1 data files  
    //
    //  - calculate hadronness for method of Nearest Neighbors
    //  - produce Neyman-Pearson plot
 
 if (JobC_NN)
 {
    cout << "=====================================================" << endl;
    cout << "Macro CT1Analysis : Start of Job C_NN" << endl;

    cout << "" << endl;
    cout << "Macro CT1Analysis : JobC_NN = " << JobC_NN  << endl;


    TString typeMC = "MC";
    cout << "typeMC = " << typeMC << endl;

    TString typeData = "ON";
    cout << "typeData = " << typeData << endl;

    TString typeMatrixHadrons = "ON";
    cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;



   //*************************************************************************
   // read in matrices of look-alike events

    const char* mtxName = "MatrixGammas";

    TString outName = outPath;
    outName += "matrix_gammas_";
    outName += "MC";
    outName += ".root";

    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Get matrix for (gammas)" << endl;
    cout << "matrix name        = " << mtxName << endl;
    cout << "name of root file  = " << outName << endl;
    cout << "" << endl;


    // read in the object with the name 'mtxName' from file 'outName'
    //
    TFile *fileg = new TFile(outName); 

    MHMatrix matrixg;
    matrixg.Read(mtxName);
    matrixg.Print("cols");

    delete fileg;


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

    const char* mtxName = "MatrixHadrons";

    TString outName = outPath;
    outName += "matrix_hadrons_";
    outName += typeMatrixHadrons;
    outName += ".root";

    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << " Get matrix for (hadrons)" << endl;
    cout << "matrix name        = " << mtxName << endl;
    cout << "name of root file  = " << outName << endl;
    cout << "" << endl;


    // read in the object with the name mtxName
    //
    TFile *fileh = new TFile(outName); 

    MHMatrix matrixh;
    matrixh.Read(mtxName);
    matrixh.Print("cols");

    delete fileh;

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


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

    MTaskList tliston;
    MParList pliston;

    //-------------------------------------------
    // create the tasks which should be executed 
    //
    TString filenameData = outPath;
    filenameData += typeData;
    filenameData += "1.root";

    TString filenameMC = outPath;
    filenameMC += typeMC;
    filenameMC += "1.root";
   
    cout << "filenameData = " << filenameData << endl;
    cout << "filenameMC   = " << filenameMC   << endl;

    MReadMarsFile read("Events", filenameMC);
    read.AddFile(filenameData);
    read.DisableAutoScheme();

    MFParticleId fgamma ("MMcEvt", '=', kGAMMA);
    MFParticleId fhadron("MMcEvt", '!', kGAMMA);

    //.......................................................................
    // calculate hadronnes for method of Nearest Neighbors

    TString hadName = "HadNN";
    MMultiDimDistCalc nncalc;
    nncalc.SetUseNumRows(25);
    nncalc.SetUseKernelMethod(kFALSE);
    nncalc.SetHadronnessName(hadName);

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

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

    TString fHilName    = "MHillas"; 
    TString fHilSrcName = "MHillasSrc"; 
    TString fnewparName = "MNewImagePar"; 

    Float_t hadcut   =  0.36;
    Float_t alphacut = 100.0;
    MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
    selfinalgh.SetCuts(hadcut, alphacut);
    selfinalgh.SetHadronnessName(hadName);
    MContinue contfinalgh(&selfinalgh);

    MFillH fillhadnn("MHHadronness", hadName);

    Float_t hadcut   = 0.36;
    Float_t alphacut = 20.0;
    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetCuts(hadcut, alphacut);
    selfinal.SetHadronnessName(hadName);
    MContinue contfinal(&selfinal);

    MFillH hfill1("MHHillas",    fHilName);
    MFillH hfill2("MHStarMap",   fHilName);

    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);

    MFillH hfill3("MHHillasExt",   fHilSrcName);
    MFillH hfill4("MHHillasSrc",   fHilSrcName);
    MFillH hfill5("MHNewImagePar", fnewparName);


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

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

    pliston.AddToList(&matrixg);
    pliston.AddToList(&matrixh);

    pliston.AddToList(&hhilext); 


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

    tliston.AddToList(&fgamma);
    tliston.AddToList(&fhadron);
    tliston.AddToList(&nncalc);
    tliston.AddToList(&fillhadnn);
   
    tliston.AddToList(&contfinalgh);
    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    tliston.AddToList(&contfinal);

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

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

    Int_t maxevents = 1000000000;
    //Int_t maxevents = 35000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);
    DeleteBinnings(&pliston);

    //-------------------------------------------
    // Display the histograms
    //
    TObject *c;
    c = pliston.FindObject("MHHadronness");
    if (c)
    {
      c->DrawClone();
      c->Print();
    }
    c = pliston.FindObject("MHHillas")->DrawClone();
    c = pliston.FindObject("MHHillasExt")->DrawClone();
    c = pliston.FindObject("MHHillasSrc")->DrawClone();
    c = pliston.FindObject("MHNewImagePar")->DrawClone();
    c = pliston.FindObject("MHStarMap")->DrawClone();



    cout << "Macro CT1Analysis : End of Job C_NN" << endl;
    cout << "===================================================" << endl;
 }



  //---------------------------------------------------------------------
  // Job C_RF
  //=========

    // read ON1 and MC1 data files  
    //
    //  - calculate hadronness for method of Random Forest
    //  - produce Neyman-Pearson plot
 
 if (JobC_RF)
 {
    cout << "=====================================================" << endl;
    cout << "Macro CT1Analysis : Start of Job C_RF" << endl;

    cout << "" << endl;
    cout << "Macro CT1Analysis : JobC_RF = " << JobC_RF  << endl;


    TString typeMC = "MC";
    cout << "typeMC = " << typeMC << endl;

    TString typeData = "ON";
    cout << "typeData = " << typeData << endl;

    TString typeMatrixHadrons = "ON";
    cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;

   //*************************************************************************
   // read in matrices of look-alike events

    const char* mtxName = "MatrixGammas";

    TString outName = outPath;
    outName += "matrix_gammas_";
    outName += "MC";
    outName += ".root";

    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Get matrix for (gammas)" << endl;
    cout << "matrix name        = " << mtxName << endl;
    cout << "name of root file  = " << outName << endl;
    cout << "" << endl;


    // read in the object with the name 'mtxName' from file 'outName'
    //
    TFile *fileg = new TFile(outName); 

    MHMatrix matrixg;
    matrixg.Read(mtxName);
    matrixg.Print("cols");

    delete fileg;


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

    const char* mtxName = "MatrixHadrons";

    TString outName = outPath;
    outName += "matrix_hadrons_";
    outName += typeMatrixHadrons;
    outName += ".root";

    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << " Get matrix for (hadrons)" << endl;
    cout << "matrix name        = " << mtxName << endl;
    cout << "name of root file  = " << outName << endl;
    cout << "" << endl;


    // read in the object with the name mtxName
    //
    TFile *fileh = new TFile(outName); 

    MHMatrix matrixh;
    matrixh.Read(mtxName);
    matrixh.Print("cols");

    delete fileh;

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

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

    cout << "Macro CT1Analysis : start growing trees" << endl;

    MTaskList tlist2;
    MParList plist2;
    plist2.AddToList(&tlist2);

    plist2.AddToList(&matrixg);
    plist2.AddToList(&matrixh);

    MRanForestGrow rfgrow2;
    rfgrow2.SetNumTrees(100);
    rfgrow2.SetNumTry(4);
    rfgrow2.SetNdSize(10);

    TString outRF = outPath;
    outRF += "RF.root";
    MWriteRootFile rfwrite2(outRF);
    rfwrite2.AddContainer("MRanTree", "TREE");

    // list of tasks for the loop over the trees
    
    tlist2.AddToList(&rfgrow2);
    tlist2.AddToList(&rfwrite2);

    //-------------------
    // Execute tree loop
    //
    MEvtLoop treeloop;
    treeloop.SetParList(&plist2);

    if ( !treeloop.Eventloop() )
        return;

    tlist2.PrintStatistics(0, kTRUE);


    // get adresses of objects which are used in the next eventloop
    MRanForest *fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
    if (!fRanForest)
    {
        *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
        return kFALSE;
    }

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

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

    MTaskList tliston;
    MParList pliston;

    //-------------------------------------------
    // create the tasks which should be executed 
    //
    TString filenameData = outPath;
    filenameData += typeData;
    filenameData += "1.root";

    TString filenameMC = outPath;
    filenameMC += typeMC;
    filenameMC += "1.root";
   
    cout << "filenameData = " << filenameData << endl;
    cout << "filenameMC   = " << filenameMC   << endl;

    MReadMarsFile read("Events", filenameMC);
    read.AddFile(filenameData);
    read.DisableAutoScheme();

    MFParticleId fgamma ("MMcEvt", '=', kGAMMA);
    MFParticleId fhadron("MMcEvt", '!', kGAMMA);

    //.......................................................................
    // calculate hadronnes for method of Random Forest

    TString hadName = "HadRF";
    MRanForestCalc rfcalc;
    rfcalc.SetHadronnessName(hadName);

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

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

    TString fHilName    = "MHillas"; 
    TString fHilSrcName = "MHillasSrc"; 
    TString fnewparName = "MNewImagePar"; 

    Float_t hadcut   =  0.36;
    Float_t alphacut = 100.0;
    MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
    selfinalgh.SetCuts(hadcut, alphacut);
    selfinalgh.SetHadronnessName(hadName);
    MContinue contfinalgh(&selfinalgh);

    MFillH fillhadrf("MHHadronness", hadName);
    MFillH fillhrf("MHRanForest");

    Float_t hadcut   = 0.36;
    Float_t alphacut = 20.0;
    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetCuts(hadcut, alphacut);
    selfinal.SetHadronnessName(hadName);
    MContinue contfinal(&selfinal);

    MFillH hfill1("MHHillas",    fHilName);
    MFillH hfill2("MHStarMap",   fHilName);

    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);

    MFillH hfill3("MHHillasExt",   fHilSrcName);
    MFillH hfill4("MHHillasSrc",   fHilSrcName);
    MFillH hfill5("MHNewImagePar", fnewparName);


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

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

    pliston.AddToList(fRanForest);
    pliston.AddToList(fRanTree);

    pliston.AddToList(&hhilext); 


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

    tliston.AddToList(&fgamma);
    tliston.AddToList(&fhadron);
    tliston.AddToList(&rfcalc);
    tliston.AddToList(&fillhadrf);
    tliston.AddToList(&fillhrf);
   
    tliston.AddToList(&contfinalgh);
    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    tliston.AddToList(&contfinal);

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

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

    Int_t maxevents = 1000000000;
    //Int_t maxevents = 35000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);
    DeleteBinnings(&pliston);

    //-------------------------------------------
    // Display the histograms
    //
    TObject *c;
    c = pliston.FindObject("MHRanForest")->DrawClone();    
    c = pliston.FindObject("MHHadronness");
    if (c)
    {
      c->DrawClone();
      c->Print();
    }

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



    cout << "Macro CT1Analysis : End of Job C_RF" << endl;
    cout << "===================================================" << endl;
 }


  //---------------------------------------------------------------------
  // Job D_XX
  //=========

    //  - select g/h separation method XX
    //  - read MC_XX2.root file 
    //  - calculate eff. collection area
    //  - read ON_XX2.root file 
    //  - apply final cuts
    //  - calculate flux
    //  - write root file for ON data after final cuts (ON_XX3.root))


 if (JobD_XX)
 {
    cout << "=====================================================" << endl;
    cout << "Macro CT1Analysis : Start of Job D_XX" << endl;

    cout << "" << endl;
    cout << "Macro CT1Analysis : JobD_XX, WXXD = " 
         << JobD_XX  << ",  " << WXXD << endl;

    // type of data to be analysed
    TString typeData = "ON";
    //TString typeData = "OFF";
    //TString typeData = "MC";
    cout << "typeData = " << typeData << endl;

    TString typeMC   = "MC";
    TString ext      = "3.root";
    TString extout   = "4.root";

    //------------------------------
    // selection of g/h separation method
    // and definition of final selections

    TString XX("NN");
    //TString XX("SC");
    //TString XX("RF");
    TString fhadronnessName("Had");
    fhadronnessName += XX;
    cout << "fhadronnessName = " << fhadronnessName << endl;

    // maximum values of the hadronness and of |alpha|
    Float_t maxhadronness   = 0.40;
    Float_t maxalpha        = 20.0;
    cout << "Maximum values of hadronness and |alpha| = "
         << maxhadronness << ",  " << maxalpha << endl;


    //------------------------------
    // name of MC file to be used for calculating the eff. collection areas
    TString filenameArea(outPath);
    filenameArea += typeMC;
    filenameArea += "_";
    filenameArea += XX;
    filenameArea += ext; 
    cout << "filenameArea = " << filenameArea << endl;

    //------------------------------
    // name of file containing the eff. collection areas
    TString collareaName(outPath);
    collareaName += "area_";
    collareaName += XX;
    collareaName += ".root";
    cout << "collareaName = " << collareaName << endl;

    //------------------------------
    // name of data file to be analysed
    TString filenameData(outPath);
    filenameData += typeData;
    filenameData += "_";
    filenameData += XX;
    filenameData += ext;
    cout << "filenameData = " << filenameData << endl;

    //------------------------------
    // name of output data file (after the final cuts)
    TString filenameDataout(outPath);
    filenameDataout += typeData;
    filenameDataout += "_";
    filenameDataout += XX;
    filenameDataout += extout;
    cout << "filenameDataout = " << filenameDataout << endl;


    //====================================================================
    cout << "-----------------------------------------------" << endl;
    cout << "Start calculation of effective collection areas" << endl;
    MParList  parlist;
    MTaskList tasklist;

    //---------------------------------------
    // Setup the tasks to be executed
    //
    MReadMarsFile reader("Events", filenameArea);

    MFCT1SelFinal cuthadrons;
    cuthadrons.SetHadronnessName(fhadronnessName);
    cuthadrons.SetCuts(maxhadronness, maxalpha);

    MContinue conthadrons(&cuthadrons);

    MHMcCT1CollectionArea* collarea = 
              new MHMcCT1CollectionArea("MHMcCT1CollectionArea","",30,2.,5.);
    MMcCT1CollectionAreaCalc effi;

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

    parlist.AddToList(&tasklist);
    parlist.AddToList(collarea);

    //********************************
    // entries in MTaskList

    tasklist.AddToList(&reader);   
    tasklist.AddToList(&conthadrons);
    tasklist.AddToList(&effi);

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

    //-----------------------------------------
    // Execute event loop
    //
    MEvtLoop magic;
    magic.SetParList(&parlist);

    MProgressBar bar;
    magic.SetProgressBar(&bar);
    if (!magic.Eventloop())
        return;

    tasklist.PrintStatistics(0, kTRUE);

    // Display the histograms
    //
    parlist.FindObject("MHMcCT1CollectionArea")->DrawClone("lego");

    //---------------------------------------------
    // Write histograms to a file 
    //

    TFile f(collareaName, "RECREATE");
    collarea->GetHist()->Write();
    collarea->GetHAll()->Write();
    collarea->GetHSel()->Write();
    f.Close();


    cout << "Calculation of effective collection areas done" << endl;
    cout << "-----------------------------------------------" << endl;    
    //------------------------------------------------------------------


    //*************************************************************************
    //
    // Analyse the data
    //

    MTaskList tliston;
    MParList pliston;

    TString hadName("Had");
    hadName += XX;

    TString fHilName    = "MHillas"; 
    TString fHilSrcName = "MHillasSrc"; 
    TString fnewparName = "MNewImagePar"; 

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

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

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

    if (WXXD)
    {    
      MWriteRootFile write(filenameDataout);

      write.AddContainer("MRawRunHeader", "RunHeaders");
      write.AddContainer("MTime",         "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("MSrcPosCam",    "Events");
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("MHillas",       "Events");
      write.AddContainer("MHillasSrc",    "Events");
      write.AddContainer("MNewImagePar",  "Events");

      write.AddContainer("HadNN",         "Events");
      write.AddContainer("HadSC",         "Events");

      write.AddContainer("MEnergyEst",    "Events");
    }

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

    MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
    selfinalgh.SetCuts(maxhadronness, 100.0);
    selfinalgh.SetHadronnessName(fhadronnessName);
    MContinue contfinalgh(&selfinalgh);

    MFillH fillhadnn("MHadNN[MHHadronness]", "HadNN");
    MFillH fillhadsc("MHadSC[MHHadronness]", "HadSC");
    //MFillH fillhadrf("MHadRF[MHHadronness]", "HadRF");

    MFillH hfill1("MHHillas",    fHilName);
    MFillH hfill2("MHStarMap",   fHilName);

    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);

    MFillH hfill3("MHHillasExt",   fHilSrcName);
    MFillH hfill4("MHHillasSrc",   fHilSrcName);
    MFillH hfill5("MHNewImagePar", fnewparName);

    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetCuts(maxhadronness, maxalpha);
    selfinal.SetHadronnessName(fhadronnessName);
    MContinue contfinal(&selfinal);


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

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

    pliston.AddToList(&hhilext); 

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

    tliston.AddToList(&contfinalgh);

    if (WXXD)
      tliston.AddToList(&write);

    tliston.AddToList(&contfinal);
    tliston.AddToList(&fillhadnn);
    tliston.AddToList(&fillhadsc);
    //tliston.AddToList(&fillhadrf);
    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);



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

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

    Int_t maxevents = 1000000000;
    //Int_t maxevents = 10000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);
    DeleteBinnings(&pliston);

    //-------------------------------------------
    // Display the histograms
    //
    TObject *c;

    c = pliston.FindObject("MHadNN", "MHHadronness");
    if (c)
    {
      c->DrawClone();
      c->Print();
    }
    c = pliston.FindObject("MHadRF", "MHHadronness");
    if (c)
    {
      c->DrawClone();
      c->Print();
    }
    c = pliston.FindObject("MHadSC", "MHHadronness");
    if (c)
    {
      c->DrawClone();
      c->Print();
    }

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


    cout << "Macro CT1Analysis : End of Job D_XX" << endl;
    cout << "=======================================================" << endl;
 }
  //---------------------------------------------------------------------

}







