void AnalyseCT1()
{
  //-----------------------------------------------
    //const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; 
    //const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; 
    //const char *mcfile = "/hd43/rwagner/CT1/PreprocData/mc_Gammas.preproc.0.6";
    //const char *mcfile = "/hd43/rwagner/mkn421_mc.preproc";
    //const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc";

  //-----------------------------------------------
  const char *offfile = "~/Mars200203/CT1data/offdata.preproc"; 

  const char *onfile  = "~/Mars200203/CT1data/mkn421_on.preproc"; 
  //const char *onfile  = "~/Mars200203/CT1data/Crab_preproc_daniel_0.6"; 

  //const char *mcfile  = "~/Mars200203/CT1data/mc_gammas.preproc"; //Version 0.5
  //const char *mcfile  = "~/Mars200203/CT1data/mc_Gammas.preproc.0.6";
  const char *mcfile = "~/Mars200203/CT1data/mkn421_mc_pedrms_0.001.preproc";

  //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6";
  //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6_040303";

  //-----------------------------------------------
  //const char *offfile = "~magican/home/ct1test/PreprocData/offdata.preproc"; 
  //const char *onfile  = "~magican/home/ct1test/PreprocData/mkn421_on.preproc"; 
  //const char *mcfile  = "~magican/home/ct1test/PreprocData/mc_gammas.preproc";


    const char *filename;

    //************************************************************************
    // Job 1 : read OFF data  
    // (generate sigmabar vs. Theta plot; write root file for OFF data (OFF1);
    //  generate hadron matrix for g/h separation)

    Bool_t Job1     = kFALSE;   
    Bool_t WHistOFF = kTRUE;   // write out histogram sigmabar vs. Theta ?
    Bool_t WLookOFF = kTRUE;   // write out look-alike events for hadrons ?
    Bool_t WOFF1    = kTRUE;   // write out root file OFF1 ?


    // Job 2 : read ON data
    // (generate sigmabar vs. Theta plot; write root file for ON data (ON1))

    Bool_t Job2    = 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 ?

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

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



    // Job 4 : read OFF1 and MC1 files, read hadron and gamma matrix
    // (produce the Neyman-Pearson plot)
    Bool_t Job4  = kTRUE;  


    // Job 5 : read MC1 file, read hadron and gamma matrix 
    // (do g/h separation; write root file for MC gammas after final cuts (MC2))
    Bool_t Job5   = kFALSE;  
    Bool_t WMC2   = kFALSE;  // write out root file MC2 ?


    // Job 6 : read ON data, read hadron and gamma matrix
    // (do g/h separation; write root file for ON data after final cuts (ON2))
    Bool_t Job6  = kFALSE;  
    Bool_t WON2  = kTRUE;  // write out root file ON2 ?
    //************************************************************************




  //---------------------------------------------------------------------
    // Job 1
    //======

    // read OFF data file 

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

    //  - write a file of look-alike events (hadron matrix)
    //    (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 OFF events (OFF1) 
    //    (after the standard cuts, before the g/h separation)
    //    (to be used together with the corresponding MC gamma file
    //     for the optimization of the g/h separation)

 if (Job1)
 {
    cout << "=====================================================" << endl;
    cout << "Macro AnalyseCT1 : Start of Job 1" << endl;

    cout << "" << endl;
    cout << "Macro AnalyseCT1 : Job1, WHistOFF, WLookOFF, WOFF1 = " 
         << Job1  << ",  " << WHistOFF << ",  " << WLookOFF << ",  " 
         << WOFF1 << endl;


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

    MTaskList tliston;
    MParList pliston;
    pliston.AddToList(&tliston);


    //::::::::::::::::::::::::::::::::::::::::::

    MBinning binssize("BinningSize");
    binssize.SetEdgesLog(50, 10, 1.0e5);
    pliston.AddToList(&binssize);

    MBinning binsdistc("BinningDist");
    binsdistc.SetEdges(50, 0, 1.4);
    pliston.AddToList(&binsdistc);

    MBinning binswidth("BinningWidth");
    binswidth.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binswidth);

    MBinning binslength("BinningLength");
    binslength.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binslength);

    MBinning binsalpha("BinningAlpha");
    binsalpha.SetEdges(100, -100, 100);
    pliston.AddToList(&binsalpha);

    MBinning binsasym("BinningAsym");
    binsasym.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsasym);

    MBinning binsht("BinningHeadTail");
    binsht.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsht);

    MBinning binsm3l("BinningM3Long");
    binsm3l.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3l);

    MBinning binsm3t("BinningM3Trans");
    binsm3t.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3t);

   
    //.....
    MBinning binsb("BinningSigmabar");
    binsb.SetEdges( 100,  0.0,  5.0);
    pliston.AddToList(&binsb);

    MBinning binth("BinningTheta");
    binth.SetEdges( 70, -0.5, 69.5);    
    pliston.AddToList(&binth);

    MBinning binsdiff("BinningDiffsigma2");
    binsdiff.SetEdges(100, -5.0, 20.0);
    pliston.AddToList(&binsdiff);
    //::::::::::::::::::::::::::::::::::::::::::


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

    filename = offfile;
    RName = "ReadCT1data";
    MCT1ReadPreProc read(filename, RName);

    MSelBasic selbasic;

    MBlindPixelCalc blind;
    blind.SetUseInterpolation();
    blind.SetUseBlindPixels();
    // blind.SetUseCetralPixel();

    // create an object of class MSigmabar, 
    // because class MHSigmaTheta will look for it
    MSigmabar sigbar;
    pliston.AddToList(&sigbar);
    MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
    fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");

    MImgCleanStd    clean; 

    // calculate also extended image parameters
    TString fHilName = "Hillas";
    MHillasExt hext(fHilName);
    pliston.AddToList(&hext);
    MHillasExt *fHillas = &hext;

    // name of output container for MHillasCalc
    //      = name of MHillas object
    MHillasCalc    hcalc(fHilName);

    // name of output container for MHillasSrcCalc
    //      = name of MHillasSrc object
    TString fHilSrcName = "HillasSrc";
    MHillasSrc hilsrc(fHilSrcName);
    MHillasSrc *fHillasSrc = &hilsrc;
    pliston.AddToList(&hilsrc);
    MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);

    MSelStandard selstand(fHilName, fHilSrcName);

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

    TString nxt = "HillasExt";
    MHHillasExt hhilext(nxt, "", fHilName);
    pliston.AddToList(&hhilext); 
    TString namext = nxt;
    namext += "[MHHillasExt]";
    MFillH hfill3(namext, fHilSrcName);

    MFillH hfill4("MHHillasSrc", fHilSrcName);

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

    const char* mtxName = "MatrixHadrons";
    Int_t howMany = 50000;
    TString outName = "matrix_hadrons_";
    outName += typeData;
    outName += ".root";

    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Matrix for (hadrons)" << endl;
    cout << "input file                    = " << filename<< 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("", "", RName);
    selector.SetNumSelectEvts(howMany);
    MFilterList flist;
    flist.AddToList(&selector);

    //
    // --- 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(Hillas.fSize)");
    matrix.AddColumn("HillasSrc.fDist");
    matrix.AddColumn("Hillas.fWidth");
    matrix.AddColumn("Hillas.fLength");
    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
    matrix.AddColumn("abs(Hillas.fM3Long)");
    matrix.AddColumn("Hillas.fConc");
    matrix.AddColumn("Hillas.fConc1");

    pliston.AddToList(pmatrix);

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


    if (WOFF1)
    {    
      TString outNameImage = typeData;
      outNameImage += "1.root";
      MWriteRootFile write(outNameImage);
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("Hillas",        "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("HillasSrc",     "Events");
      write.AddContainer("MRawRunHeader", "RunHeaders");
      //write.AddContainer("MMcRunHeader","RunHeaders");
      write.AddContainer("MSrcPosCam",    "RunHeaders");
    }

    //+++++++++++++++++++++++++++++++++++++++++++++++++++


    //*****************************
    // tasks to be executed
    
    tliston.AddToList(&read);

    tliston.AddToList(&selbasic);
    tliston.AddToList(&blind);
    tliston.AddToList(&fillsigtheta);
    tliston.AddToList(&clean);

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

    tliston.AddToList(&selstand);
    if (WOFF1)
      tliston.AddToList(&write);

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

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

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

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

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

    tliston.PrintStatistics(0, kTRUE);

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

    c = pliston.FindObject("MHHillasSrc")->DrawClone();

    //pliston.FindObject("MHHillasExt")->DrawClone();
    c = pliston.FindObject(nxt)->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 1);
    // set refcolumn negative if distribution is not to be changed
    Int_t   refcolumn = 1;
    Int_t   nbins   =  10;
    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 AnalyseCT1 : define reference sample for hadrons" << endl; 
    cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
         << refcolumn << ",  " << nmaxevts << endl;    

    if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) 
    {
      cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl;
      return;
    }
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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

      writejens->Close();
      delete writejens;
    }

    //-------------------------------------------
    // Write histograms onto a file
  if (WHistOFF)
  {
      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 = "SigmaTheta_";
      outNameSigTh += typeData;
      outNameSigTh += ".root";

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


    cout << "Macro AnalyseCT1 : End of Job 1" << endl;
    cout << "===================================================" << endl;
 }
  //---------------------------------------------------------------------

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

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

    //  - write a file of look-alike events (hadron matrix)
    //    (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) 
    //    (after the standard cuts, before the g/h separation)
    //    (to be used together with the corresponding MC gamma file
    //     for the optimization of the g/h separation)

 if (Job2)
 {
    cout << "=====================================================" << endl;
    cout << "Macro AnalyseCT1 : Start of Job 2" << endl;

    cout << "" << endl;
    cout << "Macro AnalyseCT1 : Job2, WHistON, WLookON, WON1 = " 
         << Job2  << ",  " << WHistON << ",  " << WLookON << ",  " 
         << WON1 << endl;

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

    MTaskList tliston;
    MParList pliston;
    pliston.AddToList(&tliston);


    //::::::::::::::::::::::::::::::::::::::::::

    MBinning binssize("BinningSize");
    binssize.SetEdgesLog(50, 10, 1.0e5);
    pliston.AddToList(&binssize);

    MBinning binsdistc("BinningDist");
    binsdistc.SetEdges(50, 0, 1.4);
    pliston.AddToList(&binsdistc);

    MBinning binswidth("BinningWidth");
    binswidth.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binswidth);

    MBinning binslength("BinningLength");
    binslength.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binslength);

    MBinning binsalpha("BinningAlpha");
    binsalpha.SetEdges(100, -100, 100);
    pliston.AddToList(&binsalpha);

    MBinning binsasym("BinningAsym");
    binsasym.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsasym);

    MBinning binsht("BinningHeadTail");
    binsht.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsht);

    MBinning binsm3l("BinningM3Long");
    binsm3l.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3l);

    MBinning binsm3t("BinningM3Trans");
    binsm3t.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3t);

   
    //.....
    MBinning binsb("BinningSigmabar");
    binsb.SetEdges( 100,  0.0,  5.0);
    pliston.AddToList(&binsb);

    MBinning binth("BinningTheta");
    binth.SetEdges( 70, -0.5, 69.5);    
    pliston.AddToList(&binth);

    MBinning binsdiff("BinningDiffsigma2");
    binsdiff.SetEdges(100, -5.0, 20.0);
    pliston.AddToList(&binsdiff);
    //::::::::::::::::::::::::::::::::::::::::::


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

    filename = onfile;
    RName = "ReadCT1data";
    MCT1ReadPreProc read(filename, RName);

    MSelBasic selbasic;

    MBlindPixelCalc blind;
    blind.SetUseInterpolation();
    blind.SetUseBlindPixels();
    // blind.SetUseCetralPixel();

    // create an object of class MSigmabar, 
    // because class MHSigmaTheta will look for it
    MSigmabar sigbar;
    pliston.AddToList(&sigbar);
    MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
    fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");

    MImgCleanStd    clean; 

    // calculate also extended image parameters
    TString fHilName = "Hillas";
    MHillasExt hext(fHilName);
    pliston.AddToList(&hext);
    MHillasExt *fHillas = &hext;

    // name of output container for MHillasCalc
    //      = name of MHillas object
    MHillasCalc    hcalc(fHilName);

    // name of output container for MHillasSrcCalc
    //      = name of MHillasSrc object
    TString fHilSrcName = "HillasSrc";
    MHillasSrc hilsrc(fHilSrcName);
    MHillasSrc *fHillasSrc = &hilsrc;
    pliston.AddToList(&hilsrc);
    MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);

    MSelStandard selstand(fHilName, fHilSrcName);

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

    TString nxt = "HillasExt";
    MHHillasExt hhilext(nxt, "", fHilName);
    pliston.AddToList(&hhilext); 
    TString namext = nxt;
    namext += "[MHHillasExt]";
    MFillH hfill3(namext, fHilSrcName);

    MFillH hfill4("MHHillasSrc", fHilSrcName);

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

    const char* mtxName = "MatrixHadrons";
    Int_t howMany = 50000;
    TString outName = "matrix_hadrons_";
    outName += typeData;
    outName += ".root";


    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Matrix for (hadrons)" << endl;
    cout << "input file                    = " << filename<< 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("", "", RName);
    selector.SetNumSelectEvts(howMany);
    MFilterList flist;
    flist.AddToList(&selector);

    //
    // --- 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(Hillas.fSize)");
    matrix.AddColumn("HillasSrc.fDist");
    matrix.AddColumn("Hillas.fWidth");
    matrix.AddColumn("Hillas.fLength");
    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
    matrix.AddColumn("abs(Hillas.fM3Long)");
    matrix.AddColumn("Hillas.fConc");
    matrix.AddColumn("Hillas.fConc1");

    pliston.AddToList(pmatrix);

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


    if (WON1)
    {    
      TString outNameImage = typeData;
      outNameImage += "1.root";
      MWriteRootFile write(outNameImage);
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("Hillas",        "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("HillasSrc",     "Events");
      write.AddContainer("MRawRunHeader", "RunHeaders");
      //write.AddContainer("MMcRunHeader","RunHeaders");
      write.AddContainer("MSrcPosCam",    "RunHeaders");
    }

    //+++++++++++++++++++++++++++++++++++++++++++++++++++


    //*****************************
    // tasks to be executed
    
    tliston.AddToList(&read);

    tliston.AddToList(&selbasic);
    tliston.AddToList(&blind);
    tliston.AddToList(&fillsigtheta);
    tliston.AddToList(&clean);

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

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

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

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


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

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

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

    tliston.PrintStatistics(0, kTRUE);

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

    c = pliston.FindObject("MHHillasSrc")->DrawClone();

    //pliston.FindObject("MHHillasExt")->DrawClone();
    c = pliston.FindObject(nxt)->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 1);
    // set refcolumn negative if distribution is not to be changed
    Int_t   refcolumn = 1;
    Int_t   nbins   =  10;
    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 AnalyseCT1 : define reference sample for hadrons" << endl; 
    cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
         << refcolumn << ",  " << nmaxevts << endl;    

    if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) 
    {
      cout << "Macro AnalyseCT1 : 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 << "Macro AnalyseCT1 : 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 = "SigmaTheta_";
      outNameSigTh += typeData;
      outNameSigTh += ".root";

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


    cout << "Macro AnalyseCT1 : End of Job 2" << endl;
    cout << "===================================================" << endl;
 }


  //---------------------------------------------------------------------
   // Job 3
   //======

    // read MC gamma data  

    //    - to pad them
    //      (using the 2D-histogram "sigmabar versus Theta" of the OFF data)

    //    - to write a file of look-alike events (gammas matrix)
    //      (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)
    //      (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 (Job3)
 {
    cout << "=====================================================" << endl;
    cout << "Macro AnalyseCT1 : Start of Job 3" << endl;

    cout << "" << endl;
    cout << "Macro AnalyseCT1 : Job3, WLookMC, WMC1 = " 
         << Job3  << ",  " << WLookMC << ",  " << WMC1 << endl;


    TString typeMC = "MC";
    cout << "typeMC = " << typeMC << endl;
    
    // use for padding sigmabar vs. Theta from ON or OFF data
    TString typeData = "ON";
    //TString typeData = "OFF";
    cout << "typeData = " << typeData << endl;

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


      TString outNameSigTh = "SigmaTheta_";
      outNameSigTh += typeData;
      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;

    plist.AddToList(&tlist);


    //::::::::::::::::::::::::::::::::::::::::::

    MBinning binssize("BinningSize");
    binssize.SetEdgesLog(50, 10, 1.0e5);
    plist.AddToList(&binssize);

    MBinning binsdistc("BinningDist");
    binsdistc.SetEdges(50, 0, 1.4);
    plist.AddToList(&binsdistc);

    MBinning binswidth("BinningWidth");
    binswidth.SetEdges(50, 0, 1.0);
    plist.AddToList(&binswidth);

    MBinning binslength("BinningLength");
    binslength.SetEdges(50, 0, 1.0);
    plist.AddToList(&binslength);

    MBinning binsalpha("BinningAlpha");
    binsalpha.SetEdges(100, -100, 100);
    plist.AddToList(&binsalpha);

    MBinning binsasym("BinningAsym");
    binsasym.SetEdges(50, -1.5, 1.5);
    plist.AddToList(&binsasym);

    MBinning binsht("BinningHeadTail");
    binsht.SetEdges(50, -1.5, 1.5);
    plist.AddToList(&binsht);

    MBinning binsm3l("BinningM3Long");
    binsm3l.SetEdges(50, -1.5, 1.5);
    plist.AddToList(&binsm3l);

    MBinning binsm3t("BinningM3Trans");
    binsm3t.SetEdges(50, -1.5, 1.5);
    plist.AddToList(&binsm3t);

   
    //.....
    MBinning binsb("BinningSigmabar");
    binsb.SetEdges( 100,  0.0,  5.0);
    plist.AddToList(&binsb);

    MBinning binth("BinningTheta");
    binth.SetEdges( 70, -0.5, 69.5);    
    plist.AddToList(&binth);

    MBinning binsdiff("BinningDiffsigma2");
    binsdiff.SetEdges(100, -5.0, 20.0);
    plist.AddToList(&binsdiff);

    //::::::::::::::::::::::::::::::::::::::::::

    //-------------------------------------------
    // create the tasks which should be executed 
    //
    filename = mcfile;
    readname = "ReadCT1MC";
    MCT1ReadPreProc read(filename, readname);

    MSelBasic selbasic;
   
    MBlindPixelCalc blind;
    blind.SetUseInterpolation();
    blind.SetUseBlindPixels();
    // blind.SetUseCetralPixel();
    
    // 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);

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

    MImgCleanStd    clean; //(0, 0);

    // calculate also extended image parameters
    TString fHilName = "Hillas";
    MHillasExt hext(fHilName);
    plist.AddToList(&hext);
    MHillasExt *fHillas = &hext;

    // name of output container for MHillasCalc
    //      = name of MHillas object
    MHillasCalc    hcalc(fHilName);

    // name of output container for MHillasSrcCalc
    //      = name of MHillasSrc object
    TString fHilSrcName = "HillasSrc";
    MHillasSrc hilsrc(fHilSrcName);
    MHillasSrc *fHillasSrc = &hilsrc;
    plist.AddToList(&hilsrc);
    MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);

    MSelStandard selstand(fHilName, fHilSrcName);

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

    TString nxt = "HillasExt";
    MHHillasExt hhilext(nxt, "", fHilName);
    plist.AddToList(&hhilext); 
    TString namext = nxt;
    namext += "[MHHillasExt]";
    MFillH hfill3(namext, fHilSrcName);

    MFillH hfill4("MHHillasSrc", fHilSrcName);

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

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

    TString outName = "matrix_gammas_";
    outName += typeMC;
    outName += "_";
    outName += typeData;
    outName += ".root";


    cout << "" << endl;
    cout << "========================================================" << endl;
    cout << "Matrix for (gammas)" << endl;
    cout << "input file                    = " << filename<< 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("", "", readname);
    selector.SetNumSelectEvts(howMany);
    MFilterList flist;
    flist.AddToList(&selector);

    //
    // --- 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(Hillas.fSize)");
    matrix.AddColumn("HillasSrc.fDist");
    matrix.AddColumn("Hillas.fWidth");
    matrix.AddColumn("Hillas.fLength");
    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
    matrix.AddColumn("abs(Hillas.fM3Long)");
    matrix.AddColumn("Hillas.fConc");
    matrix.AddColumn("Hillas.fConc1");

    plist.AddToList(pmatrix);

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


    if (WMC1)
    {    
      TString outNameImage = typeMC;
      outNameImage += "_";
      outNameImage += typeData;
      outNameImage += "1.root";
      MWriteRootFile write(outNameImage);
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("Hillas",        "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("HillasSrc",     "Events");
      write.AddContainer("MRawRunHeader", "RunHeaders");
      //write.AddContainer("MMcRunHeader","RunHeaders");
      write.AddContainer("MSrcPosCam",    "RunHeaders");
    }


    //+++++++++++++++++++++++++++++++++++++++++++++++++++


    //*****************************
    // tasks to be executed

    tlist.AddToList(&read);

    tlist.AddToList(&selbasic);
    tlist.AddToList(&blind);
    tlist.AddToList(&padthomas);
    tlist.AddToList(&clean);

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

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

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

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


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


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

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

    tlist.PrintStatistics(0, kTRUE);

    //-------------------------------------------
    // Display the histograms
    //
    plist.FindObject("MHHillas")->DrawClone();
    plist.FindObject("MHHillasSrc")->DrawClone();

    //plist.FindObject("MHHillasExt")->DrawClone();
    plist.FindObject(nxt)->DrawClone();

    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;    
    // set refcolumn negative if distribution is not to be changed
    Int_t   refcolumn = 1;
    Int_t   nbins   =  10;
    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 AnalyseCT1 : define reference sample for gammas" << endl; 
    cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
         << refcolumn << ",  " << nmaxevts << endl;    

    if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) 
    {
      cout << "Macro AnalyseCT1 : 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 AnalyseCT1 : matrix of look-alike events for gammas written onto file "
           << outName << endl;

      writejens->Close();
      delete writejens;
    }

    cout << "Macro AnalyseCT1 : End of Job 3" 
         << endl;
    cout << "=========================================================" 
         << endl;
 }



  //---------------------------------------------------------------------
  // Job 4
  //======

    // read OFF1 and MC1 data files  
    //
    //  - produce Neyman-Pearson plot
 
 if (Job4)
 {
    cout << "=====================================================" << endl;
    cout << "Macro AnalyseCT1 : Start of Job 4" << endl;

    cout << "" << endl;
    cout << "Macro AnalyseCT1 : Job4 = " << Job4  << endl;


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

    // use hadron matrix from ON or OFF data
    TString typeData = "ON";
    //TString typeData = "OFF";
    cout << "typeData = " << typeData << endl;

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

    const char* mtxName = "MatrixGammas";


    TString outName = "matrix_gammas_";
    outName += typeMC;
    outName += "_";
    outName += typeData;
    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);
    pmatrixg = &matrixg;

    delete fileg;


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

    const char* mtxName = "MatrixHadrons";

    TString outName = "matrix_hadrons_";
    outName += typeData;
    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);
    pmatrixh = &matrixh;

    delete fileh;

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

    MTaskList tliston;
    MParList pliston;
    pliston.AddToList(&tliston);

    pliston.AddToList(pmatrixg);
    pliston.AddToList(pmatrixh);

    //::::::::::::::::::::::::::::::::::::::::::

    MBinning binssize("BinningSize");
    binssize.SetEdgesLog(50, 10, 1.0e5);
    pliston.AddToList(&binssize);

    MBinning binsdistc("BinningDist");
    binsdistc.SetEdges(50, 0, 1.4);
    pliston.AddToList(&binsdistc);

    MBinning binswidth("BinningWidth");
    binswidth.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binswidth);

    MBinning binslength("BinningLength");
    binslength.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binslength);

    MBinning binsalpha("BinningAlpha");
    binsalpha.SetEdges(100, -100, 100);
    pliston.AddToList(&binsalpha);

    MBinning binsasym("BinningAsym");
    binsasym.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsasym);

    MBinning binsht("BinningHeadTail");
    binsht.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsht);

    MBinning binsm3l("BinningM3Long");
    binsm3l.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3l);

    MBinning binsm3t("BinningM3Trans");
    binsm3t.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3t);

   
    //.....
    MBinning binsb("BinningSigmabar");
    binsb.SetEdges( 100,  0.0,  5.0);
    pliston.AddToList(&binsb);

    MBinning binth("BinningTheta");
    binth.SetEdges( 70, -0.5, 69.5);    
    pliston.AddToList(&binth);

    MBinning binsdiff("BinningDiffsigma2");
    binsdiff.SetEdges(100, -5.0, 20.0);
    pliston.AddToList(&binsdiff);
    //::::::::::::::::::::::::::::::::::::::::::


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

    TString filenameData = typeData;
    filenameData += "1.root";

    TString filenameMC = typeMC;
    filenameMC += "_";
    filenameMC += typeData;
    filenameMC += "1.root";
   

    cout << "filenameData = " << filenameData << endl;
    cout << "filenameMC   = " << filenameMC   << endl;

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

    //.......................................................................
    // g/h separation

    MMultiDimDistCalc ghsep;
    ghsep.SetUseNumRows(25);
    ghsep.SetUseKernelMethod(kFALSE);

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

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

    TString fHilName    = "Hillas"; 
    TString fHilSrcName = "HillasSrc"; 

    Float_t hadcut = 0.2;
    MSelFinal selfinalgh(fHilName, fHilSrcName);
    selfinalgh.SetHadronnessCut(hadcut);
    selfinalgh.SetAlphaCut(100.0);

    MFillH fillh("MHHadronness");

    MSelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetHadronnessCut(hadcut);
    selfinal.SetAlphaCut(20.0);

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

    TString nxt = "HillasExt";
    MHHillasExt hhilext(nxt, "", fHilName);
    pliston.AddToList(&hhilext); 
    TString namext = nxt;
    namext += "[MHHillasExt]";
    MFillH hfill3(namext, fHilSrcName);

    MFillH hfill4("MHHillasSrc", fHilSrcName);


    //*****************************
    // tasks to be executed
    
    tliston.AddToList(&read);

    tliston.AddToList(&ghsep);
    tliston.AddToList(&fillh);
   
    tliston.AddToList(&selfinalgh);
    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);

    tliston.AddToList(&selfinal);

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

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

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

    tliston.PrintStatistics(0, kTRUE);


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

    c = pliston.FindObject("MHHadronness")->DrawClone();
    c->Print();

    //c = pliston.FindObject("MHHillas")->DrawClone();

    c = pliston.FindObject("MHHillasSrc")->DrawClone();

    //pliston.FindObject("MHHillasExt")->DrawClone();
    //c = pliston.FindObject(nxt)->DrawClone();

    c = pliston.FindObject("MHStarMap")->DrawClone();



    cout << "Macro AnalyseCT1 : End of Job 4" << endl;
    cout << "===================================================" << endl;
 }


  //---------------------------------------------------------------------
  // Job 5
  //======

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

    // then read MC1 data file 
    //
    //  - perform the g/h separation
    //  - apply the final cuts
    //
    //  - write a file of MC gamma events (MC2)
    //    (after the g/h separation and after the final cuts)

 if (Job5)
 {
    cout << "=====================================================" << endl;
    cout << "Macro AnalyseCT1 : Start of Job 5" << endl;

    cout << "" << endl;
    cout << "Macro AnalyseCT1 : Job5, WMC2 = " 
         << Job5  << ",  " << WMC2 << endl;

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

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

    // use hadron matrix from ON or OFF data
    TString typeData = "ON";
    //TString typeData = "OFF";
    cout << "typeData = " << typeData << endl;

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

    const char* mtxName = "MatrixGammas";

    TString outName = "matrix_gammas_";
    outName += typeMC;
    outName += "_";
    outName += typeData;
    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);
    pmatrixg = &matrixg;

    delete fileg;

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

    const char* mtxName = "MatrixHadrons";

    TString outName = "matrix_hadrons_";
    outName += typeData;
    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);
    pmatrixh = &matrixh;

    delete fileh;

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


    MTaskList tliston;
    MParList pliston;
    pliston.AddToList(&tliston);

    pliston.AddToList(pmatrixg);
    pliston.AddToList(pmatrixh);

    //::::::::::::::::::::::::::::::::::::::::::

    MBinning binssize("BinningSize");
    binssize.SetEdgesLog(50, 10, 1.0e5);
    pliston.AddToList(&binssize);

    MBinning binsdistc("BinningDist");
    binsdistc.SetEdges(50, 0, 1.4);
    pliston.AddToList(&binsdistc);

    MBinning binswidth("BinningWidth");
    binswidth.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binswidth);

    MBinning binslength("BinningLength");
    binslength.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binslength);

    MBinning binsalpha("BinningAlpha");
    binsalpha.SetEdges(100, -100, 100);
    pliston.AddToList(&binsalpha);

    MBinning binsasym("BinningAsym");
    binsasym.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsasym);

    MBinning binsht("BinningHeadTail");
    binsht.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsht);

    MBinning binsm3l("BinningM3Long");
    binsm3l.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3l);

    MBinning binsm3t("BinningM3Trans");
    binsm3t.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3t);

   
    //.....
    MBinning binsb("BinningSigmabar");
    binsb.SetEdges( 100,  0.0,  5.0);
    pliston.AddToList(&binsb);

    MBinning binth("BinningTheta");
    binth.SetEdges( 70, -0.5, 69.5);    
    pliston.AddToList(&binth);

    MBinning binsdiff("BinningDiffsigma2");
    binsdiff.SetEdges(100, -5.0, 20.0);
    pliston.AddToList(&binsdiff);
    //::::::::::::::::::::::::::::::::::::::::::


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

    TString filenameMC = typeInput;
    filenameMC += "_";
    filenameMC += typeData;
    filenameMC += "1.root";

    readname = "ReadCT1MCdata";
    MReadMarsFile read("Events", filenameMC);
    read.DisableAutoScheme();


    //.......................................................................
    // g/h separation

    MMultiDimDistCalc ghsep;
    ghsep.SetUseNumRows(25);
    ghsep.SetUseKernelMethod(kFALSE);



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

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

    TString fHilName    = "Hillas"; 
    TString fHilSrcName = "HillasSrc"; 

    Float_t hadcut = 0.2;
    MSelFinal selfinalgh(fHilName, fHilSrcName);
    selfinalgh.SetHadronnessCut(hadcut);
    selfinalgh.SetAlphaCut(100.0);

    MFillH fillh("MHHadronness");

    MSelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetHadronnessCut(hadcut);
    selfinal.SetAlphaCut(20.0);

    if (WMC2)
    {    
      TString outNameImage = typeInput;
      outNameImage += "_";
      outNameImage += typeData;
      outNameImage += "2.root";
      MWriteRootFile write(outNameImage);
      write.AddContainer("MHadronness",   "Events");
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("Hillas",        "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("HillasSrc",     "Events");
      write.AddContainer("MRawRunHeader", "RunHeaders");
      //write.AddContainer("MMcRunHeader","RunHeaders");
      write.AddContainer("MSrcPosCam",    "RunHeaders");
    }


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

    TString nxt = "HillasExt";
    MHHillasExt hhilext(nxt, "", fHilName);
    pliston.AddToList(&hhilext); 
    TString namext = nxt;
    namext += "[MHHillasExt]";
    MFillH hfill3(namext, fHilSrcName);

    MFillH hfill4("MHHillasSrc", fHilSrcName);



    //*****************************
    // tasks to be executed
    
    tliston.AddToList(&read);

    tliston.AddToList(&ghsep);
    tliston.AddToList(&fillh);

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

    tliston.AddToList(&selfinal);
    if (WMC2)
      tliston.AddToList(&write);

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

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

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

    tliston.PrintStatistics(0, kTRUE);


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

    c = pliston.FindObject("MHHadronness")->DrawClone();
    c->Print();

    c = pliston.FindObject("MHHillas")->DrawClone();

    c = pliston.FindObject("MHHillasSrc")->DrawClone();

    //pliston.FindObject("MHHillasExt")->DrawClone();
    c = pliston.FindObject(nxt)->DrawClone();

    c = pliston.FindObject("MHStarMap")->DrawClone();



    cout << "Macro AnalyseCT1 : End of Job 5" << endl;
    cout << "===================================================" << endl;
 }


  //---------------------------------------------------------------------
  // Job 6
  //======

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

    // then read ON1 data file 
    //
    //  - perform the g/h separation
    //  - apply the final cuts
    //
    //  - write a file of ON events (ON2)
    //    (after the g/h separation and after the final cuts)
    //    (to be used for the flux calculation)
    //
    //  - do the flux calculation

 if (Job6)
 {
    cout << "=====================================================" << endl;
    cout << "Macro AnalyseCT1 : Start of Job 6" << endl;

    cout << "" << endl;
    cout << "Macro AnalyseCT1 : Job6, WON2 = " 
         << Job6  << ",  " << WON2 << endl;

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

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

    // use hadron matrix from ON or OFF data
    TString typeData = "ON";
    //TString typeData = "OFF";
    cout << "typeData = " << typeData << endl;


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

    const char* mtxName = "MatrixGammas";

    TString outName = "matrix_gammas_";
    outName += typeMC;
    outName += "_";
    outName += typeData;
    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);
    pmatrixg = &matrixg;

    delete fileg;

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

    const char* mtxName = "MatrixHadrons";

    TString outName = "matrix_hadrons_";
    outName += typeData;
    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);
    pmatrixh = &matrixh;

    delete fileh;

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


    MTaskList tliston;
    MParList pliston;
    pliston.AddToList(&tliston);

    pliston.AddToList(pmatrixg);
    pliston.AddToList(pmatrixh);

    //::::::::::::::::::::::::::::::::::::::::::

    MBinning binssize("BinningSize");
    binssize.SetEdgesLog(50, 10, 1.0e5);
    pliston.AddToList(&binssize);

    MBinning binsdistc("BinningDist");
    binsdistc.SetEdges(50, 0, 1.4);
    pliston.AddToList(&binsdistc);

    MBinning binswidth("BinningWidth");
    binswidth.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binswidth);

    MBinning binslength("BinningLength");
    binslength.SetEdges(50, 0, 1.0);
    pliston.AddToList(&binslength);

    MBinning binsalpha("BinningAlpha");
    binsalpha.SetEdges(100, -100, 100);
    pliston.AddToList(&binsalpha);

    MBinning binsasym("BinningAsym");
    binsasym.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsasym);

    MBinning binsht("BinningHeadTail");
    binsht.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsht);

    MBinning binsm3l("BinningM3Long");
    binsm3l.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3l);

    MBinning binsm3t("BinningM3Trans");
    binsm3t.SetEdges(50, -1.5, 1.5);
    pliston.AddToList(&binsm3t);

   
    //.....
    MBinning binsb("BinningSigmabar");
    binsb.SetEdges( 100,  0.0,  5.0);
    pliston.AddToList(&binsb);

    MBinning binth("BinningTheta");
    binth.SetEdges( 70, -0.5, 69.5);    
    pliston.AddToList(&binth);

    MBinning binsdiff("BinningDiffsigma2");
    binsdiff.SetEdges(100, -5.0, 20.0);
    pliston.AddToList(&binsdiff);
    //::::::::::::::::::::::::::::::::::::::::::


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

    TString filenameData = typeInput;
    filenameData += "1.root";

    readname = "ReadCT1ONdata";
    MReadMarsFile read("Events", filenameData);
    read.DisableAutoScheme();


    //.......................................................................
    // g/h separation

    MMultiDimDistCalc ghsep;
    ghsep.SetUseNumRows(25);
    ghsep.SetUseKernelMethod(kFALSE);



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

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

    TString fHilName    = "Hillas"; 
    TString fHilSrcName = "HillasSrc"; 

    Float_t hadcut = 0.2;
    MSelFinal selfinalgh(fHilName, fHilSrcName);
    selfinalgh.SetHadronnessCut(hadcut);
    selfinalgh.SetAlphaCut(100.0);

    MFillH fillh("MHHadronness");

    MSelFinal selfinal(fHilName, fHilSrcName);
    selfinal.SetHadronnessCut(hadcut);
    selfinal.SetAlphaCut(20.0);

    if (WON2)
    {    
      TString outNameImage = typeInput;
      outNameImage += "_";
      outNameImage += typeData;
      outNameImage += "2.root";
      MWriteRootFile write(outNameImage);
      write.AddContainer("MHadronness",   "Events");
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("Hillas",        "Events");
      write.AddContainer("MMcEvt",        "Events");
      write.AddContainer("HillasSrc",     "Events");
      write.AddContainer("MRawRunHeader", "RunHeaders");
      //write.AddContainer("MMcRunHeader","RunHeaders");
      write.AddContainer("MSrcPosCam",    "RunHeaders");
    }


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

    TString nxt = "HillasExt";
    MHHillasExt hhilext(nxt, "", fHilName);
    pliston.AddToList(&hhilext); 
    TString namext = nxt;
    namext += "[MHHillasExt]";
    MFillH hfill3(namext, fHilSrcName);

    MFillH hfill4("MHHillasSrc", fHilSrcName);



    //*****************************
    // tasks to be executed
    
    tliston.AddToList(&read);

    tliston.AddToList(&ghsep);
    tliston.AddToList(&fillh);

    tliston.AddToList(&selfinalgh);

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

    tliston.AddToList(&selfinal);
    if (WON2)
      tliston.AddToList(&write);



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

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

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

    tliston.PrintStatistics(0, kTRUE);


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

    c = pliston.FindObject("MHHadronness")->DrawClone();
    c->Print();

    //c = pliston.FindObject("MHHillas")->DrawClone();

    c = pliston.FindObject("MHHillasSrc")->DrawClone();

    ////pliston.FindObject("MHHillasExt")->DrawClone();
    //c = pliston.FindObject(nxt)->DrawClone();

    c = pliston.FindObject("MHStarMap")->DrawClone();



    cout << "Macro AnalyseCT1 : End of Job 6" << endl;
    cout << "=======================================================" << endl;
 }
  //---------------------------------------------------------------------
}










