Index: trunk/MagicSoft/Mars/macros/AnalyseCT1.C
===================================================================
--- trunk/MagicSoft/Mars/macros/AnalyseCT1.C	(revision 1800)
+++ trunk/MagicSoft/Mars/macros/AnalyseCT1.C	(revision 1809)
@@ -1,51 +1,112 @@
 void AnalyseCT1()
 {
-    const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; 
-    const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; 
+  //-----------------------------------------------
+    //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 *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;
 
-    Bool_t Data  = kFALSE;   // loop over data ?
-    Bool_t WLook = kFALSE;   // write out look-alike events for hadrons ?
-    Bool_t WHist = kFALSE;   // write out histograms for padding ?
-
-    Bool_t MC      = kTRUE;  // loop over MC gamma data ?
-    Bool_t RHist   = kTRUE;  // read in histograms for padding ?
-    Bool_t WLookMC = kTRUE;  // write out look-alike events for gammas ?
-
-
-    cout << "" << endl;
-    cout << "Macro AnalyseCT1 : Data, WHist = " << Data << ",  "
-         << WHist << endl;
-
-    cout << "                         RHist, MC   = " << RHist << ",  "
-         << MC << endl;
-    cout << "" << endl;
+    //************************************************************************
+    // 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 ?
+    //************************************************************************
+
+
+
 
   //---------------------------------------------------------------------
-  // start of section for ON data
-  // (in the CT1 analysis the ON data are also used as a sample representing
-  //  the hadrons for the g/h separation)  
-
-    // read ON data file 
-
-    //  - to produce the 2D-histogram "sigmabar versus Theta"
+    // 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)
 
-    //  - to write a file of look-alike events (for hadrons)
+    //  - 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)
 
-    //  - to write a file of ON events (after the standard cuts, 
-    //                                  before the g/h separation)
+    //  - 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 (Data)
+ if (Job1)
  {
     cout << "=====================================================" << endl;
-    cout << "Macro AnalyseCT1 : Start of section for ON data" << 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;
@@ -112,7 +173,7 @@
     //
 
-    filename = onfile;
-    readname = "ReadCT1data";
-    MCT1ReadPreProc read(filename, readname);
+    filename = offfile;
+    RName = "ReadCT1data";
+    MCT1ReadPreProc read(filename, RName);
 
     MSelBasic selbasic;
@@ -150,5 +211,5 @@
     MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
 
-    MSelStandard selstand(fHillas, fHillasSrc);
+    MSelStandard selstand(fHilName, fHilSrcName);
 
     MFillH hfill1("MHHillas",    fHilName);
@@ -168,6 +229,8 @@
 
     const char* mtxName = "MatrixHadrons";
-    Int_t howMany = 30000;
-    char* outName = "matrix_hadrons.root";
+    Int_t howMany = 50000;
+    TString outName = "matrix_hadrons_";
+    outName += typeData;
+    outName += ".root";
 
     cout << "" << endl;
@@ -182,5 +245,5 @@
 
     // matrix limitation for look-alike events (approximate number)
-    MFEventSelector selector("", "", readname);
+    MFEventSelector selector("", "", RName);
     selector.SetNumSelectEvts(howMany);
     MFilterList flist;
@@ -192,18 +255,16 @@
     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("Hillas.fWidth*Hillas.fLength/Hillas.fSize");
-    matrix.AddColumn("abs(Hillas.fAsym)");
+    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
     matrix.AddColumn("abs(Hillas.fM3Long)");
-    matrix.AddColumn("abs(Hillas.fM3Trans)");
-    matrix.AddColumn("abs(HillasSrc.fHeadTail)");
     matrix.AddColumn("Hillas.fConc");
     matrix.AddColumn("Hillas.fConc1");
-    matrix.AddColumn("HillasSrc.fDist");
-    matrix.AddColumn("log10(Hillas.fSize)");
-    matrix.AddColumn("HillasSrc.fAlpha");
-    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
-    //matrix.AddColumn("MMcEvt.fTheta");
+
     pliston.AddToList(pmatrix);
 
@@ -211,7 +272,21 @@
     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");
+    }
+
     //+++++++++++++++++++++++++++++++++++++++++++++++++++
 
-    MSelFinal selfinal(fHillas, fHillasSrc);
 
     //*****************************
@@ -229,4 +304,10 @@
 
     tliston.AddToList(&selstand);
+    if (WOFF1)
+      tliston.AddToList(&write);
+
+    tliston.AddToList(&flist);
+    tliston.AddToList(&fillmatg);    
+
     tliston.AddToList(&hfill1);
     tliston.AddToList(&hfill2);
@@ -234,8 +315,4 @@
     tliston.AddToList(&hfill4);
 
-    tliston.AddToList(&flist);
-    tliston.AddToList(&fillmatg);
-
-    //tliston.AddToList(&selfinal);
     //*****************************
 
@@ -243,8 +320,11 @@
     // 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;
@@ -255,11 +335,13 @@
     // Display the histograms
     //
-    pliston.FindObject("MHHillas")->DrawClone();
-    pliston.FindObject("MHHillasSrc")->DrawClone();
+    TObject *c;
+    c = pliston.FindObject("MHHillas")->DrawClone();
+
+    c = pliston.FindObject("MHHillasSrc")->DrawClone();
 
     //pliston.FindObject("MHHillasExt")->DrawClone();
-    pliston.FindObject(nxt)->DrawClone();
-
-    pliston.FindObject("MHStarMap")->DrawClone();
+    c = pliston.FindObject(nxt)->DrawClone();
+
+    c = pliston.FindObject("MHStarMap")->DrawClone();
 
 
@@ -277,7 +359,7 @@
     Int_t nmaxevts  = 2000;
 
-    // target distribution for the variable in column refcolumn;    
+    // target distribution for the variable in column refcolumn (start at 1);
     // set refcolumn negative if distribution is not to be changed
-    Int_t   refcolumn = 12;
+    Int_t   refcolumn = 1;
     Int_t   nbins   =  10;
     Float_t frombin = 0.5;
@@ -311,5 +393,5 @@
     // write out look-alike events (for hadrons)
     //
-    if (WLook)
+    if (WLookOFF)
     {
       TFile *writejens = new TFile(outName, "RECREATE", "");
@@ -324,5 +406,5 @@
     //-------------------------------------------
     // Write histograms onto a file
-  if (WHist)
+  if (WHistOFF)
   {
       MHSigmaTheta *sigtheta = 
@@ -330,5 +412,5 @@
       if (!sigtheta)
 	{
-          *fLog << "Object 'SigmaTheta' not found" << endl;
+          cout << "Object 'SigmaTheta' not found" << endl;
           return;
 	}
@@ -337,5 +419,9 @@
       TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
 
-      TFile outfile("SigmaTheta.root", "RECREATE");
+      TString outNameSigTh = "SigmaTheta_";
+      outNameSigTh += typeData;
+      outNameSigTh += ".root";
+
+      TFile outfile(outNameSigTh, "RECREATE");
       fHSigTh->Write();
       fHSigPixTh->Write();
@@ -343,9 +429,9 @@
       outfile.Close();
      
-      cout << "File 'SigmaTheta.root' was written out" << endl;
+      cout << "File " << outNameSigTh << " was written out" << endl;
   }
 
 
-    cout << "Macro AnalyseCT1 : End of section for ON data" << endl;
+    cout << "Macro AnalyseCT1 : End of Job 1" << endl;
     cout << "===================================================" << endl;
  }
@@ -353,5 +439,366 @@
 
   //---------------------------------------------------------------------
-  // start of section for MC gamma data
+  // 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  
@@ -360,28 +807,31 @@
     //      (using the 2D-histogram "sigmabar versus Theta" of the OFF data)
 
-    //    - to write a file of look-alike events (for gammas)
+    //    - 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 
+    //    - 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)
 
-    //    - to write a file of padded MC gamma events (after the final cuts) 
-    //      (to be used for the optimization of the energy estimator
-    //       and for calculating the effective collection areas)
-
- if (MC)
+
+ if (Job3)
  {
-    cout << "=============================================================" 
-         << endl;
-    cout << "Macro : AnalyseCT1 : Start of section for MC gamma data" 
-         << endl;
-
-    MTaskList tlist;
-    MParList plist;
-
-    plist.AddToList(&tlist);
+    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;
 
     //------------------------------------
@@ -389,9 +839,14 @@
     // and                "3D-ThetaPixSigma"
     // and                "3D-ThetaPixDiff"
-  if (RHist)
-  {
-      cout << "Reading in file 'SigmaTheta.root' " << endl;
-
-      TFile *infile = new TFile("SigmaTheta.root");
+
+
+      TString outNameSigTh = "SigmaTheta_";
+      outNameSigTh += typeData;
+      outNameSigTh += ".root";
+
+
+      cout << "Reading in file " << outNameSigTh << endl;
+
+      TFile *infile = new TFile(outNameSigTh);
       infile->ls();
 
@@ -400,5 +855,5 @@
       if (!fHSigmaTheta)
 	{
-          *fLog << "Object '2D-ThetaSigmabar' not found on root file" << endl;
+          cout << "Object '2D-ThetaSigmabar' not found on root file" << endl;
           return;
 	}
@@ -409,5 +864,5 @@
       if (!fHSigmaPixTheta)
 	{
-          *fLog << "Object '3D-ThetaPixSigma' not found on root file" << endl;
+          cout << "Object '3D-ThetaPixSigma' not found on root file" << endl;
           return;
 	}
@@ -418,24 +873,15 @@
       if (!fHSigmaPixTheta)
 	{
-          *fLog << "Object '3D-ThetaPixDiff' not found on root file" << endl;
+          cout << "Object '3D-ThetaPixDiff' not found on root file" << endl;
           return;
 	}
       cout << "Object '3D-ThetaPixDiff' was read in" << endl;
-  }
-  else
-  {
-      MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
-      if (!sigtheta)
-      {
-        cout << "Object with name 'SigmaTheta' not found" << endl;
-        return;
-      }
-
-      TH2D *fHSigmaTheta    = sigtheta->GetSigmaTheta();
-      TH3D *fHSigmaPixTheta = sigtheta->GetSigmaPixTheta();
-      TH3D *fHDiffPixTheta  = sigtheta->GetDiffPixTheta();
-  }
+
     //------------------------------------
 
+    MTaskList tlist;
+    MParList plist;
+
+    plist.AddToList(&tlist);
 
 
@@ -539,5 +985,5 @@
     MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
 
-    MSelStandard selstand(fHillas, fHillasSrc);
+    MSelStandard selstand(fHilName, fHilSrcName);
 
     MFillH hfill1("MHHillas",    fHilName);
@@ -557,6 +1003,12 @@
 
     const char* mtxName = "MatrixGammas";
-    Int_t howMany = 30000;
-    char* outName = "matrix_gammas.root";
+    Int_t howMany = 50000;
+
+    TString outName = "matrix_gammas_";
+    outName += typeMC;
+    outName += "_";
+    outName += typeData;
+    outName += ".root";
+
 
     cout << "" << endl;
@@ -581,18 +1033,16 @@
     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("Hillas.fWidth*Hillas.fLength/Hillas.fSize");
-    matrix.AddColumn("abs(Hillas.fAsym)");
+    matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
     matrix.AddColumn("abs(Hillas.fM3Long)");
-    matrix.AddColumn("abs(Hillas.fM3Trans)");
-    matrix.AddColumn("abs(HillasSrc.fHeadTail)");
     matrix.AddColumn("Hillas.fConc");
     matrix.AddColumn("Hillas.fConc1");
-    matrix.AddColumn("HillasSrc.fDist");
-    matrix.AddColumn("log10(Hillas.fSize)");
-    matrix.AddColumn("HillasSrc.fAlpha");
-    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
-    //matrix.AddColumn("MMcEvt.fTheta");
+
     plist.AddToList(pmatrix);
 
@@ -600,8 +1050,24 @@
     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");
+    }
+
+
     //+++++++++++++++++++++++++++++++++++++++++++++++++++
 
-
-    MSelFinal selfinal(fHillas, fHillasSrc);
 
     //*****************************
@@ -619,4 +1085,10 @@
 
     tlist.AddToList(&selstand);
+    if (WMC1)
+      tlist.AddToList(&write);
+
+    tlist.AddToList(&flist);
+    tlist.AddToList(&fillmatg);
+
     tlist.AddToList(&hfill1);
     tlist.AddToList(&hfill2);
@@ -624,8 +1096,5 @@
     tlist.AddToList(&hfill4);
 
-    tlist.AddToList(&flist);
-    tlist.AddToList(&fillmatg);
-
-    //tlist.AddToList(&selfinal);
+
     //*****************************
 
@@ -634,8 +1103,11 @@
     // 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;
@@ -670,5 +1142,5 @@
     // target distribution for the variable in column refcolumn;    
     // set refcolumn negative if distribution is not to be changed
-    Int_t   refcolumn = 12;
+    Int_t   refcolumn = 1;
     Int_t   nbins   =  10;
     Float_t frombin = 0.5;
@@ -711,9 +1183,859 @@
     }
 
-    cout << "Macro AnalyseCT1 : End of section for MC gamma data" 
+    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;
+ }
+  //---------------------------------------------------------------------
 }
 
@@ -721,2 +2043,8 @@
 
 
+
+
+
+
+
+
