Changeset 2475 for trunk/MagicSoft/Mars/macros

Timestamp:

11/05/03 15:11:25 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/macros

Files:

• CT1Analysis.C (modified) (7 diffs)
• ONOFFCT1Analysis.C (modified) (35 diffs)

trunk/MagicSoft/Mars/macros/CT1Analysis.C

@@ -241 +241 @@
 
     Bool_t JobB_RF_UP  = kFALSE;  
+    Bool_t CTrainRF    = kFALSE;  // create matrices of training events
     Bool_t RTrainRF    = kFALSE;  // read in matrices of training events
-    Bool_t CTrainRF    = kFALSE;  // create matrices of training events
-    Bool_t RTree       = kFALSE;  // read in trees
+    Bool_t RTree       = kFALSE;  // read in trees (otherwise grow them)
     Bool_t WRF         = kFALSE;  // update input root file ?
 
@@ -255 +255 @@
 
     Bool_t JobB_SC_UP  = kFALSE;
-    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file  
+    Bool_t CMatrix     = kFALSE;  // create training and test matrices 
+    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file
     Bool_t WOptimize   = kFALSE;  // do optimization using the training sample
                                   // and write supercuts parameter values 
@@ -1293 +1294 @@
     rfwrite2.AddContainer("MRanTree", "TREE");
 
+    MFillH fillh2("MHRanForestGini");
+
     // list of tasks for the loop over the trees
     
     tlist2.AddToList(&rfgrow2);
     tlist2.AddToList(&rfwrite2);
+    tlist2.AddToList(&fillh2);
 
     //-------------------
@@ -1309 +1313 @@
 
     tlist2.PrintStatistics(0, kTRUE);
+
+    plist2.FindObject("MHRanForestGini")->DrawClone();
+
 
     // get adresses of objects which are used in the next eventloop
@@ -1557 +1564 @@
 
     gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobB_SC_UP, RMatrix, WOptimize, RTest, WSC = "
+    gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = "
          << (JobB_SC_UP ? "kTRUE" : "kFALSE")  << ",  " 
+         << (CMatrix    ? "kTRUE" : "kFALSE")  << ",  "
          << (RMatrix    ? "kTRUE" : "kFALSE")  << ",  "
          << (WOptimize  ? "kTRUE" : "kFALSE")  << ",  "
@@ -1660 +1668 @@
     //--------------------------
     // create matrices and write them onto files 
-    if (!RMatrix)
+    if (CMatrix)
     {
       MH3 &mh3 = *(new MH3("MHillas.fSize"));
@@ -2111 +2119 @@
     //  - read ON1 and MC1 data files  
     //    which should have been updated to contain the hadronnesses
-    //    for the method of NEAREST NEIGHBORS and for the SUOERCUTS
+    //    for the method of Random Forest and for the SUPERCUTS
     //  - produce Neyman-Pearson plots
  

trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C

@@ +1 @@
+
 
 #include "CT1EgyEst.C"
@@ -7 +8 @@
         gLog << "InitBinnings" << endl;
 
+        //--------------------------------------------
         MBinning *binse = new MBinning("BinningE");
         //binse->SetEdgesLog(30, 1.0e2, 1.0e5);
-        binse->SetEdges(30, 2, 5);
+
+        //This is Daniel's binning in energy:
+        binse->SetEdgesLog(14, 296.296, 86497.6);
         plist->AddToList(binse);
+
+        //--------------------------------------------
 
         MBinning *binssize = new MBinning("BinningSize");
@@ -51 +57 @@
 
         MBinning *binth = new MBinning("BinningTheta");
+        // this is Daniel's binning in theta
+        //Double_t yedge[8] = 
+        //  {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99};
+        // this is our binning
         Double_t yedge[9] = 
                        {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
-        TArrayD yed(9,yedge);
+        TArrayD yed;
+        yed.Set(9,yedge);
         binth->SetEdges(yed);
         plist->AddToList(binth);
@@ -63 +74 @@
           zedge[8-i] = cos(yedge[i]/kRad2Deg);
         }
-        TArrayD zed(9,zedge);
+        TArrayD zed;
+        zed.Set(9,zedge);
         bincosth->SetEdges(zed);
         plist->AddToList(bincosth);
@@ -70 +82 @@
         binsdiff->SetEdges(100, -5.0, 20.0);
         plist->AddToList(binsdiff);
+
+        // robert ----------------------------------------------
+        MBinning *binsalphaf = new MBinning("BinningAlphaFlux");
+        binsalphaf->SetEdges(100, -100, 100);
+        plist->AddToList(binsalphaf);
+
+        MBinning *binsdifftime = new MBinning("BinningTimeDiff");
+        binsdifftime->SetEdges(50, 0, 10);
+        plist->AddToList(binsdifftime);
+
+        MBinning *binstime = new MBinning("BinningTime");
+        binstime->SetEdges(50, 44500, 61000);
+        plist->AddToList(binstime);
 }
+
 
 void DeleteBinnings(MParList *plist)
@@ -76 +102 @@
         gLog << "DeleteBinnings" << endl;
 
-        if (!plist)
-        {
-          gLog << "Deletebinnins : MParlist no longer existing" << endl;
-          return;
-        }
-
         TObject *bin;
+
+        //--------------------------------------------
+        bin = plist->FindObject("BinningE");
+        if (bin) delete bin;
+
+        //--------------------------------------------
 
         bin = plist->FindObject("BinningSize");
@@ -121 +147 @@
         if (bin) delete bin;
 
-        gLog << "exit DeleteBinnings" << endl;
+
+        // robert ----------------------------------------------
+        bin = plist->FindObject("BinningAlphaFlux");
+        if (bin) delete bin;
+
+        bin = plist->FindObject("BinningTimeDiff");
+        if (bin) delete bin;
+
+        bin = plist->FindObject("BinningTime");
+        if (bin) delete bin;
 }
+
+
 
 //************************************************************************
@@ -162 +199 @@
     //  - write root file for ON (or OFF or MC) data (ON1.root, ...);
 
-    Bool_t JobA    = kTRUE;  
-    Bool_t WPad    = kTRUE;   // write out padding histograms ?
+    Bool_t JobA    = kFALSE;  
+    Bool_t WPad    = kFALSE;   // write out padding histograms ?
     Bool_t RPad    = kFALSE;   // read in padding histograms ?
     Bool_t Wout    = kFALSE;   // write out root file ON1.root 
@@ -178 +215 @@
 
     Bool_t JobB_RF_UP  = kFALSE;  
-    Bool_t RTrainRF     = kFALSE;  // read in matrices of training events
-    Bool_t CTrainRF     = kFALSE;  // create  matrices of training events
-    Bool_t RTree       = kFALSE;  // read in trees
+    Bool_t CTrainRF    = kFALSE;  // create  matrices of training events
+    Bool_t RTrainRF    = kFALSE;  // read in matrices of training events
+    Bool_t RTree       = kFALSE;  // read in trees (otherwise grow them)
     Bool_t WRF         = kFALSE;  // update input root file ?
 
 
+    // Job B_SC_UP : read ON2.root (or MC2.root) file 
+    //  - depending on WParSC : create (or read in) supercuts parameter values
+    //  - calculate hadroness for the SUPERCUTS
+    //  - update the input files with the hadroness (==>ON3.root or MC3.root)
+
+    Bool_t JobB_SC_UP  = kFALSE;
+    Bool_t CMatrix     = kFALSE;  // create training and test matrices 
+    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file
+    Bool_t WOptimize   = kFALSE;  // do optimization using the training sample
+                                  // and write supercuts parameter values 
+                                  // onto the file parSCfile
+    Bool_t RTest       = kFALSE;  // test the supercuts using the test matrix
+    Bool_t WSC         = kFALSE;  // update input root file ?
+
 
     // Job C: 
-    //  - read ON1.root and MC1.root files
+    //  - read ON3.root and MC3.root files
     //    which should have been updated to contain the hadronnesses  
     //    for the method of 
-    //              NEAREST NEIGHBORS  
+    //              RF
     //              SUPERCUTS and
-    //              RF
     //  - produce Neyman-Pearson plots
 
@@ -199 +249 @@
     // Job D :  
     //  - select g/h separation method XX
-    //  - read ON2 (or MC2) root file
+    //  - read ON3 (or MC3) root file
     //  - apply cuts in hadronness
     //  - make plots
@@ -207 +257 @@
 
 
-    // Job E_EST_UP : 
-    //  - read MC1.root file 
+    // Job E_XX : extended version of E_XX (including flux plots)  
     //  - select g/h separation method XX
-    //  - optimize energy estimation for events passing the final cuts
-    //  - write parameters of energy estimator onto file
-    //  - update ON1.root, OFF1.root and MC1.root files with estimated energy
-    //    (ON_XX1.root, OFF_XX1.root and MC_XX1.root)
-
-    Bool_t JobE_EST_UP  = kFALSE;  
-    Bool_t WESTUP       = kFALSE;  // update root files ?
-
-
-
-    // Job F_XX :  
-    //  - select g/h separation method XX
-    //  - read MC_XX2.root file 
+    //  - read MC root file 
     //  - calculate eff. collection area
-    //  - read ON_XX2.root file 
+    //  - optimize energy estimator
+    //  - read ON root file 
     //  - apply final cuts
     //  - calculate flux
-    //  - write root file for ON data after final cuts (ON3.root))
-
-    Bool_t JobF_XX  = kFALSE;  
-    Bool_t WXX      = kFALSE;  // write out root file ON3.root ?
+    //  - write root file for ON data after final cuts 
+
+
+    Bool_t JobE_XX  = kTRUE;  
+    Bool_t OEEst    = kFALSE;  // optimize energy estimator
+    Bool_t WEX      = kTRUE;  // update root file  ?
+    Bool_t WRobert  = kTRUE;  // write out Robert's file  ?
+
 
 
@@ -277 +319 @@
     //--------------------------------------------------
     // type of data to be padded 
-    TString typeInput = "ON";
+    //TString typeInput = "ON";
     //TString typeInput = "OFF";
-    //TString typeInput = "MC";
+    TString typeInput = "MC";
     gLog << "typeInput = " << typeInput << endl;
 
@@ -740 +782 @@
     // file to be updated (ON, OFF or MC)
 
-    TString typeInput = "ON";
+    //TString typeInput = "ON";
     //TString typeInput = "OFF";
-    //TString typeInput = "MC";
+    TString typeInput = "MC";
     gLog << "typeInput = " << typeInput << endl;
 
@@ -937 +979 @@
     MProgressBar matrixbar;
     MEvtLoop evtloopg;
+    evtloopg.SetName("FillGammaMatrix");
     evtloopg.SetParList(&plistg);
     //evtloopg.ReadEnv(env, "", printEnv);
@@ -997 +1040 @@
     MProgressBar matrixbar;
     MEvtLoop evtlooph;
+    evtlooph.SetName("FillHadronMatrix");
     evtlooph.SetParList(&plisth);
     //evtlooph.ReadEnv(env, "", printEnv);
@@ -1009 +1053 @@
 
 
-    // write out matrices of training events events
+    // write out matrices of training events 
 
     gLog << "" << endl;
@@ -1119 +1163 @@
     rfwrite2.AddContainer("MRanTree", "TREE");
 
+    MFillH fillh2("MHRanForestGini");
+
     // list of tasks for the loop over the trees
     
     tlist2.AddToList(&rfgrow2);
     tlist2.AddToList(&rfwrite2);
+    tlist2.AddToList(&fillh2);
 
     //-------------------
@@ -1135 +1182 @@
 
     tlist2.PrintStatistics(0, kTRUE);
+
+    plist2.FindObject("MHRanForestGini")->DrawClone();
+
 
     // get adresses of objects which are used in the next eventloop
@@ -1212 +1262 @@
 
     //-----------------------------------------------------------------
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
 
 
@@ -1304 +1350 @@
     MProgressBar bar;
     MEvtLoop evtloop;
+    evtloop.SetName("UpdateRootFile");
     evtloop.SetParList(&pliston);
     evtloop.SetProgressBar(&bar);
@@ -1373 +1420 @@
 
 
-
   //---------------------------------------------------------------------
-  // Job C  
-  //======
-
-    //  - read ON1 and MC1 data files  
-    //    which should have been updated to contain the hadronnesses
-    //    for the method of NEAREST NEIGHBORS and for the SUOERCUTS
-    //  - produce Neyman-Pearson plots
- 
- if (JobC)
+  // Job B_SC_UP
+  //============
+
+    //  - create (or read in) optimum supercuts parameter values
+    //
+    //  - calculate the hadroness for the supercuts
+    //
+    //  - update input root file, including the hadroness
+
+
+ if (JobB_SC_UP)
  {
     gLog << "=====================================================" << endl;
-    gLog << "Macro CT1Analysis : Start of Job C" << endl;
+    gLog << "Macro CT1Analysis : Start of Job B_SC_UP" << endl;
 
     gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobC = " << JobC  << endl;
-
-
-    // name of input data file
+    gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = "
+         << (JobB_SC_UP ? "kTRUE" : "kFALSE")  << ",  " 
+         << (CMatrix    ? "kTRUE" : "kFALSE")  << ",  "
+         << (RMatrix    ? "kTRUE" : "kFALSE")  << ",  "
+         << (WOptimize  ? "kTRUE" : "kFALSE")  << ",  "
+         << (RTest      ? "kTRUE" : "kFALSE")  << ",  "
+         << (WSC        ? "kTRUE" : "kFALSE")  << endl;
+
+
+    //--------------------------------------------
+    // file which contains the initial parameter values for the supercuts 
+    // if parSCinit ="" the initial values are taken from the constructor of
+    //                  MCT1Supercuts
+
+    TString parSCinit = outPath;
+    //parSCinit += "parSC_1709d";
+    parSCinit = "";
+
+    gLog << "parSCinit = " << parSCinit << endl;
+
+    //---------------
+    // file onto which the optimal parameter values for the supercuts 
+    // are written
+
+    TString parSCfile = outPath;
+    parSCfile += "parSC_2310a";
+
+    gLog << "parSCfile = " << parSCfile << endl;
+
+    //--------------------------------------------
+    // file to be updated (either ON or MC)
+
+    //TString typeInput = "ON";
+    //TString typeInput = "OFF";
+    TString typeInput = "MC";
+    gLog << "typeInput = " << typeInput << endl;
+
+    // name of input root file
     TString filenameData = outPath;
-    filenameData += "OFF";
+    filenameData += typeInput;
     filenameData += "2.root";
-    gLog << "filenameData = " << filenameData << endl;
-
-    // name of input MC file
-    TString filenameMC = outPath;
-    filenameMC += "MC";
-    filenameMC += "2.root";
-    gLog << "filenameMC   = " << filenameMC   << endl;
+    gLog << "filenameData = " << filenameData << endl; 
+
+    // name of output root file
+    TString outNameImage = outPath;
+    outNameImage += typeInput;
+    outNameImage += "3.root";
+    
+
+    //TString outNameImage = filenameData;
+
+    gLog << "outNameImage = " << outNameImage << endl; 
+
+    //--------------------------------------------
+    // files to be read for optimizing the supercuts
+    // 
+    // for the training
+    TString filenameTrain = outPath;
+    filenameTrain += "ON";
+    filenameTrain += "1.root";
+    Int_t howManyTrain = 800000;
+    gLog << "filenameTrain = " << filenameTrain << ",   howManyTrain = "
+         << howManyTrain  << endl; 
+
+    // for testing
+    TString filenameTest = outPath;
+    filenameTest += "ON";
+    filenameTest += "1.root";
+    Int_t howManyTest = 800000;
+
+    gLog << "filenameTest = " << filenameTest << ",   howManyTest = "
+         << howManyTest  << endl; 
+    
+
+    //--------------------------------------------
+    // files to contain the matrices (generated from filenameTrain and
+    //                                               filenameTest)
+    // 
+    // for the training
+    TString fileMatrixTrain = outPath;
+    fileMatrixTrain += "MatrixTrainSC";
+    fileMatrixTrain += ".root";
+    gLog << "fileMatrixTrain = " << fileMatrixTrain << endl; 
+
+    // for testing
+    TString fileMatrixTest = outPath;
+    fileMatrixTest += "MatrixTestSC";
+    fileMatrixTest += ".root";
+    gLog << "fileMatrixTest = " << fileMatrixTest << endl; 
+
+    
+
+    //---------------------------------------------------------------------
+    // Training and test matrices :
+    // - either create them and write them onto a file
+    // - or read them from a file
+
+
+    MCT1FindSupercuts findsuper;
+    findsuper.SetFilenameParam(parSCfile);
+    findsuper.SetHadronnessName("HadSC");
+
+
+    //--------------------------
+    // create matrices and write them onto files 
+    if (CMatrix)
+    {
+      MH3 &mh3 = *(new MH3("MHillas.fSize"));
+      mh3.SetName("Target distribution for SIZE");
+
+      if (filenameTrain == filenameTest)
+      {
+        if ( !findsuper.DefineTrainTestMatrix(filenameTrain, 
+                              howManyTrain, mh3, howManyTest,  mh3,
+                              fileMatrixTrain, fileMatrixTest)     )
+        {
+          *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl;
+          return;
+        }
+
+      }
+      else
+      {
+        if ( !findsuper.DefineTrainMatrix(filenameTrain, 
+                              howManyTrain, mh3, fileMatrixTrain) )
+        {
+          *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl;
+          return;
+        }
+
+        if ( !findsuper.DefineTestMatrix( filenameTest,  
+                              howManyTest,  mh3, fileMatrixTest)  )
+        {
+          *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl;
+          return;
+        }
+      }
+     }
+
+    //--------------------------
+    // read matrices from files
+    //                              
+
+    if (RMatrix)
+      findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
+    //--------------------------
+
+
+
+    //---------------------------------------------------------------------
+    // optimize supercuts using the training sample
+    // 
+    // the initial values are taken 
+    //     - from the file parSCinit (if != "")
+    //     - or from the arrays params and steps (if their sizes are != 0)
+    //     - or from the MCT1Supercuts constructor
+
+
+if (WOptimize)
+  {
+    gLog << "CT1Analysis.C : optimize the supercuts using the training matrix" 
+         << endl;
+
+    TArrayD params(0);
+    TArrayD steps(0);
+  
+    if (parSCinit == "")
+    {
+      Double_t vparams[104] = {
+      // LengthUp
+        0.315585,  0.001455, 0.203198, 0.005532, -0.001670, -0.020362,
+        0.007388, -0.013463,
+      // LengthLo
+        0.151530,  0.028323, 0.510707, 0.053089,  0.013708,  2.357993,
+        0.000080, -0.007157,
+      // WidthUp
+        0.145412, -0.001771, 0.054462, 0.022280, -0.009893,  0.056353,
+        0.020711, -0.016703,
+      // WidthLo
+        0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101,
+        0.006126, -0.002849,
+      // DistUp
+        1.787943,  0.0,      2.942310, 0.199815,  0.0,       0.249909,
+        0.189697,  0.0,
+      // DistLo
+        0.589406,  0.0,     -0.083964,-0.007975,  0.0,       0.045374,
+       -0.001750,  0.0,
+      // AsymUp
+        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // AsymLo
+       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // ConcUp
+        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // ConcLo
+       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // Leakage1Up
+        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // Leakage1Lo
+       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // AlphaUp
+        13.12344,  0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0                                                 };
+
+      Double_t vsteps[104] = {
+      // LengthUp
+        0.03,      0.0002,   0.02,     0.0006,    0.0002,    0.002,
+        0.0008,    0.002,
+      // LengthLo
+        0.02,      0.003,    0.05,     0.006,     0.002,     0.3,
+        0.0001,    0.0008,
+      // WidthUp
+        0.02,      0.0002,   0.006,    0.003,     0.002,     0.006,
+        0.002,     0.002,
+      // WidthLo
+        0.009,     0.0007,   0.008,    0.0004,    0.0005,    0.002,
+        0.0007,    0.003,
+      // DistUp
+        0.2,       0.0,      0.3,      0.02,      0.0,       0.03,
+        0.02,      0.0
+      // DistLo
+        0.06,      0.0,      0.009,    0.0008,    0.0,       0.005,
+        0.0002,    0.0
+      // AsymUp  
+        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // AsymLo  
+        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // ConcUp  
+        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // ConcLo  
+        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // Leakage1Up  
+        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // Leakage1Lo  
+        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0,
+      // AlphaUp  
+        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
+        0.0,       0.0                                                 };
+
+      params.Set(104, vparams);
+      steps.Set (104, vsteps );
+    }
+
+    Bool_t rf;
+    rf = findsuper.FindParams(parSCinit, params, steps);
+
+    if (!rf) 
+    {
+       gLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
+       return;
+    }
+  }
+
+    //--------------------------------------
+    // test the supercuts on the test sample
+    //    
+
+ if (RTest)
+ {
+    gLog << "CT1Analysis.C : test the supercuts on the test matrix" << endl;
+    Bool_t rt = findsuper.TestParams();
+    if (!rt) 
+    {
+       gLog << "CT1Analysis.C : test of supercuts on the test matrix failed" 
+            << endl;
+    }
+
+ }
 
 
     //-----------------------------------------------------------------
-
+    // Update the input files with the SC hadronness
+    //
+
+ if (WSC)
+ {
+    gLog << "" << endl;
+    gLog << "========================================================" << endl;
+    gLog << "Update input file '" <<  filenameData 
+         << "' with the SC hadronness" << endl;
+
+
+    //----------------------------------------------------
+    // read in optimum parameter values for the supercuts
+
+    TFile inparam(parSCfile);
+    MCT1Supercuts scin;
+    scin.Read("MCT1Supercuts");
+    inparam.Close();
+
+    gLog << "Parameter values for supercuts were read in from file '"
+         << parSCfile << "'" << endl;
+
+    TArrayD supercutsPar;
+    supercutsPar =  scin.GetParameters();
+
+    TArrayD supercutsStep;
+    supercutsStep =  scin.GetStepsizes();
+
+    gLog << "Parameter values for supercuts : " << endl;
+    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
+    {
+      gLog << supercutsPar[i] << ",  ";
+    }
+    gLog << endl;
+
+    gLog << "Step values for supercuts : " << endl;
+    for (Int_t i=0; i<supercutsStep.GetSize(); i++)
+    {
+      gLog << supercutsStep[i] << ",  ";
+    }
+    gLog << endl;
+
+
+    //----------------------------------------------------
     MTaskList tliston;
     MParList pliston;
+
+    // set the parameters of the supercuts
+    MCT1Supercuts supercuts;
+    supercuts.SetParameters(supercutsPar);
+    gLog << "parameter values for the supercuts used for updating the input file ' " 
+         << filenameData << "'" << endl;
+    supercutsPar = supercuts.GetParameters();
+    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
+    {
+      gLog << supercutsPar[i] << ",  ";
+    }
+    gLog << endl;
 
 
@@ -1419 +1787 @@
     //
 
-    MReadMarsFile read("Events", filenameMC);
-    read.AddFile(filenameData);
+    MReadMarsFile read("Events", filenameData);
     read.DisableAutoScheme();
-
-
-    //.......................................................................
-    // names of hadronness containers
-
-    TString hadNNName = "HadNN";
-    TString hadSCName = "HadSC";
-    TString hadRFName = "HadRF";
-
-    //.......................................................................
-
 
     TString fHilName    = "MHillas"; 
@@ -1439 +1795 @@
     TString fImgParName = "MNewImagePar"; 
 
-    Float_t maxhadronness =  0.40;
-    Float_t maxalpha      =  20.0;
-    Float_t maxdist       =  10.0;
-
-    MFCT1SelFinal selfinalgh(fHilNameSrc);
-    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
-    selfinalgh.SetHadronnessName(hadRFName);
-    selfinalgh.SetName("SelFinalgh");
-    MContinue contfinalgh(&selfinalgh);
-    contfinalgh.SetName("ContSelFinalgh");
-
-    //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
-    //fillhadnn.SetName("HhadNN");
-    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
-    fillhadsc.SetName("HhadSC");
-    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
-    fillhadrf.SetName("HhadRF");
-
-    MFCT1SelFinal selfinal(fHilNameSrc);
-    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
-    selfinal.SetHadronnessName(hadRFName);
-    selfinal.SetName("SelFinal");
-    MContinue contfinal(&selfinal);
-    contfinal.SetName("ContSelFinal");
-
-
-    MFillH hfill1("MHHillas",    fHilName);
-    hfill1.SetName("HHillas");
-
-    MFillH hfill2("MHStarMap",   fHilName);
-    hfill2.SetName("HStarMap");
-
-    MFillH hfill3("MHHillasExt",    fHilNameSrc);
-    hfill3.SetName("HHillasExt");
+
+    //.......................................................................
+    // calculation of hadroness for the supercuts
+    // (=0.25 if fullfilled, =0.75 otherwise)
+
+    TString hadSCName = "HadSC";
+    MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc);
+    sccalc.SetHadronnessName(hadSCName);
+
+
+    //.......................................................................
+
+
+      //MWriteRootFile write(outNameImage, "UPDATE");
+      //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
+
     
-    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
-    hfill4.SetName("HHillasSrc");    
-
-    MFillH hfill5("MHNewImagePar", fImgParName);
-    hfill5.SetName("HNewImagePar");
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-
-    //tliston.AddToList(&fillhadnn);
-    tliston.AddToList(&fillhadsc);
-    tliston.AddToList(&fillhadrf);
-   
-    tliston.AddToList(&contfinalgh);
-    tliston.AddToList(&hfill1);
-    tliston.AddToList(&hfill2);
-    tliston.AddToList(&hfill3);
-    tliston.AddToList(&hfill4);
-    tliston.AddToList(&hfill5);
-
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-
-    Int_t maxevents = -1;
-    //Int_t maxevents = 35000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-
-    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
-    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
-    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
-
-    pliston.FindObject("MHHillas")->DrawClone();
-    pliston.FindObject("MHHillasExt")->DrawClone();
-    pliston.FindObject("MHHillasSrc")->DrawClone();
-    pliston.FindObject("MHNewImagePar")->DrawClone();
-    pliston.FindObject("MHStarMap")->DrawClone();
-
-    DeleteBinnings(&pliston);
-
-    gLog << "Macro CT1Analysis : End of Job C" << endl;
-    gLog << "===================================================" << endl;
- }
-
-
-  //---------------------------------------------------------------------
-  // Job D
-  //======
-
-    //  - select g/h separation method XX
-    //  - read ON2 (or MC2) root file 
-    //  - apply cuts in hadronness
-    //  - make plots
-
-
- if (JobD)
- {
-    gLog << "=====================================================" << endl;
-    gLog << "Macro CT1Analysis : Start of Job D" << endl;
-
-    gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobD = " 
-         << JobD  << endl;
-
-    // type of data to be analysed
-    TString typeData = "ON";
-    //TString typeData = "OFF";
-    //TString typeData = "MC";
-    gLog << "typeData = " << typeData << endl;
-
-    TString ext      = "2.root";
-
-
-    //------------------------------
-    // selection of g/h separation method
-    // and definition of final selections
-
-    //TString XX("NN");
-    //TString XX("SC");
-    TString XX("RF");
-    TString fhadronnessName("Had");
-    fhadronnessName += XX;
-    gLog << "fhadronnessName = " << fhadronnessName << endl;
-
-    // maximum values of the hadronness, |ALPHA| and DIST
-    Float_t maxhadronness   = 0.30;
-    Float_t maxalpha        = 20.0;
-    Float_t maxdist         = 10.0;
-    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
-         << maxhadronness << ",  " << maxalpha << ",  " 
-         << maxdist << endl;
-
-
-    //------------------------------
-    // name of data file to be analysed
-    TString filenameData(outPath);
-    filenameData += typeData;
-    filenameData += ext;
-    gLog << "filenameData = " << filenameData << endl;
-
-
-
-    //*************************************************************************
-    //
-    // Analyse the data
-    //
-
-    MTaskList tliston;
-    MParList pliston;
-
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-
-    TString fHilName    = "MHillas"; 
-    TString fHilNameExt = "MHillasExt"; 
-    TString fHilNameSrc = "MHillasSrc"; 
-    TString fImgParName = "MNewImagePar"; 
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameData);
-    read.DisableAutoScheme();
-
-
-    //-----------------------------------------------------------------
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-    MFCT1SelFinal selfinalgh(fHilNameSrc);
-    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
-    selfinalgh.SetHadronnessName(fhadronnessName);
-    selfinalgh.SetName("SelFinalgh");
-    MContinue contfinalgh(&selfinalgh);
-    contfinalgh.SetName("ContSelFinalgh");
-
-    //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
-    //fillhadnn.SetName("HhadNN");
-    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
-    fillhadsc.SetName("HhadSC");
-    MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
-    fillhadrf.SetName("HhadRF");
-
-
-    MFillH hfill1("MHHillas",    fHilName);
-    hfill1.SetName("HHillas");
-
-    MFillH hfill2("MHStarMap",   fHilName);
-    hfill2.SetName("HStarMap");
-
-    MFillH hfill3("MHHillasExt",   fHilNameSrc);
-    hfill3.SetName("HHillasExt");    
-
-    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
-    hfill4.SetName("HHillasSrc");    
-
-    MFillH hfill5("MHNewImagePar", fImgParName);
-    hfill5.SetName("HNewImagePar");    
-
-    MFCT1SelFinal selfinal(fHilNameSrc);
-    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
-    selfinal.SetHadronnessName(fhadronnessName);
-    selfinal.SetName("SelFinal");
-    MContinue contfinal(&selfinal);
-    contfinal.SetName("ContSelFinal");
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-
-    tliston.AddToList(&contfinalgh);
-
-    //tliston.AddToList(&fillhadnn);
-    tliston.AddToList(&fillhadsc);
-    tliston.AddToList(&fillhadrf);
-
-    tliston.AddToList(&hfill1);
-    tliston.AddToList(&hfill2);
-    tliston.AddToList(&hfill3);
-    tliston.AddToList(&hfill4);
-    tliston.AddToList(&hfill5);
-
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-
-    Int_t maxevents = -1;
-    //Int_t maxevents = 10000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-
-
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-
-    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
-    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
-    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
-
-    pliston.FindObject("MHHillas")->DrawClone();
-    pliston.FindObject("MHHillasExt")->DrawClone();
-    pliston.FindObject("MHHillasSrc")->DrawClone();
-    pliston.FindObject("MHNewImagePar")->DrawClone();
-    pliston.FindObject("MHStarMap")->DrawClone();
-
-    //-------------------------------------------
-    // fit alpha distribution to get the number of excess events
-    //
-
-    MHHillasSrc* hillasSrc = 
-      (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
-    TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
-  
-    MHOnSubtraction onsub;
-    onsub.Calc(alphaHist, &pliston, kTRUE, 13.1);
-    //-------------------------------------------
-
-    DeleteBinnings(&pliston);
-
-    gLog << "Macro CT1Analysis : End of Job D" << endl;
-    gLog << "=======================================================" << endl;
- }
-  //---------------------------------------------------------------------
-
-
-  //---------------------------------------------------------------------
-  // Job E_EST_UP
-  //============
-
-    //  - read MC2.root file 
-    //  - select g/h separation method XX
-    //  - optimize energy estimator for events passing final cuts
-    //  - write parameters of energy estimator onto file "energyest_XX.root"
-    //
-    //  - read ON2.root and MC2.root files 
-    //  - update input root file with the estimated energy
-    //    (ON_XX2.root, MC_XX2.root)
-
-
- if (JobE_EST_UP)
- {
-    gLog << "=====================================================" << endl;
-    gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl;
-
-    gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = " 
-         << JobE_EST_UP  << ",  " << WESTUP << endl;
-
-
-    TString typeON  = "ON";
-    TString typeOFF = "OFF";
-    TString typeMC  = "MC";
-    TString ext    = "2.root";
-    TString extout = "3.root";
-
-    //------------------------------
-    // name of MC file to be used for optimizing the energy estimator
-    TString filenameOpt(outPath);
-    filenameOpt += typeMC;
-    filenameOpt += ext; 
-    gLog << "filenameOpt = " << filenameOpt << endl;
-
-    //------------------------------
-    // selection of g/h separation method
-    // and definition of final selections
-
-    //TString XX("NN");
-    //TString XX("SC");
-    TString XX("RF");
-    TString fhadronnessName("Had");
-    fhadronnessName += XX;
-    gLog << "fhadronnessName = " << fhadronnessName << endl;
-
-    // maximum values of the hadronness, |alpha| and dist
-    Float_t maxhadronness   = 0.40;
-    Float_t maxalpha        = 20.0;
-    Float_t maxdist         = 10.0;
-    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
-         << maxhadronness << ",  " << maxalpha << ",  " 
-         << maxdist << endl;
-
-    // name of file containing the parameters of the energy estimator
-    TString energyParName(outPath);
-    energyParName += "energyest_";
-    energyParName += XX;
-    energyParName += ".root";
-    gLog << "energyParName = " << energyParName << endl;
-
-
-    //------------------------------
-    // name of ON file to be updated
-    TString filenameON(outPath);
-    filenameON += typeON;
-    filenameON += ext;
-    gLog << "filenameON = " << filenameON << endl;
-
-    // name of OFF file to be updated
-    TString filenameOFF(outPath);
-    filenameOFF += typeOFF;
-    filenameOFF += ext;
-    gLog << "filenameOFF = " << filenameOFF << endl;
-
-    // name of MC file to be updated
-    TString filenameMC(outPath);
-    filenameMC += typeMC;
-    filenameMC += ext;
-    gLog << "filenameMC = " << filenameMC << endl;
-
-    //------------------------------
-    // name of updated ON file 
-    TString filenameONup(outPath);
-    filenameONup += typeON;
-    filenameONup += "_";
-    filenameONup += XX;
-    filenameONup += extout;
-    gLog << "filenameONup = " << filenameONup << endl;
-
-    // name of updated OFF file 
-    TString filenameOFFup(outPath);
-    filenameOFFup += typeOFF;
-    filenameOFFup += "_";
-    filenameOFFup += XX;
-    filenameOFFup += extout;
-    gLog << "filenameOFFup = " << filenameOFFup << endl;
-
-    // name of updated MC file 
-    TString filenameMCup(outPath);
-    filenameMCup += typeMC;
-    filenameMCup += "_";
-    filenameMCup += XX;
-    filenameMCup += extout;
-    gLog << "filenameMCup = " << filenameMCup << endl;
-
-    //-----------------------------------------------------------
-
-    TString fHilName    = "MHillas"; 
-    TString fHilNameExt = "MHillasExt"; 
-    TString fHilNameSrc = "MHillasSrc"; 
-    TString fImgParName = "MNewImagePar"; 
-
-    //===========================================================
-    //
-    // Optimization of energy estimator
-    //
-
-    TString inpath("");
-    TString outpath("");
-    Int_t howMany = 2000;
-    CT1EEst(inpath,   filenameOpt,   outpath, energyParName, 
-            fHilName, fHilNameSrc,   fhadronnessName,
-            howMany,  maxhadronness, maxalpha, maxdist);
-
-    //-----------------------------------------------------------
-    //
-    // Read in parameters of energy estimator
-    //
-    gLog << "================================================================"
-         << endl;
-    gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '"
-         << energyParName << "'" << endl;
-    TFile enparam(energyParName);
-    MMcEnergyEst mcest("MMcEnergyEst"); 
-    mcest.Read("MMcEnergyEst");
-    enparam.Close();
-
-    TArrayD parA(5);
-    TArrayD parB(7);
-    for (Int_t i=0; i<parA.GetSize(); i++)
-      parA[i] = mcest.GetCoeff(i);
-    for (Int_t i=0; i<parB.GetSize(); i++)
-      parB[i] = mcest.GetCoeff( i+parA.GetSize() );
-
-
-   if (WESTUP)
-   {
-    //==========   start update   ============================================
-    //
-    // Update ON, OFF and MC root files with the estimated energy
-
-    //---------------------------------------------------
-    // Update OFF data 
-    //
-    gLog << "============================================================"
-         << endl;
-    gLog << "Macro CT1Analysis.C : update file '" << filenameOFF
-         << "'" << endl;
-
-    MTaskList tlistoff;
-    MParList plistoff;
-
-
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)plistoff->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameOFF);
-    read.DisableAutoScheme();
-
-    //---------------------------
-    // calculate estimated energy
-
-    MEnergyEstParam eest2(fHilName);
-    eest2.Add(fHilNameSrc);
-
-    eest2.SetCoeffA(parA);
-    eest2.SetCoeffB(parB);
-
-    //.......................................................................
-
-      MWriteRootFile write(filenameOFFup);
+      MWriteRootFile write(outNameImage, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -1951 +1825 @@
       write.AddContainer("MNewImagePar",  "Events");
 
-      //write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
       write.AddContainer("HadRF",         "Events");
-
-      write.AddContainer("MEnergyEst",    "Events");
+      write.AddContainer(hadSCName,       "Events");
+    
 
     //-----------------------------------------------------------------
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    Float_t maxhadronness =  0.40;
+    Float_t maxalpha      =  20.0;
+    Float_t maxdist       =  10.0;
+
+    MFCT1SelFinal selfinalgh(fHilNameSrc);
+    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
+    selfinalgh.SetHadronnessName(hadSCName);
+    selfinalgh.SetName("SelFinalgh");
+    MContinue contfinalgh(&selfinalgh);
+    contfinalgh.SetName("ContSelFinalgh");
+
+    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
+    fillhadsc.SetName("HhadSC");
+
+    MFCT1SelFinal selfinal(fHilNameSrc);
+    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
+    selfinal.SetHadronnessName(hadSCName);
+    selfinal.SetName("SelFinal");
+    MContinue contfinal(&selfinal);
+    contfinal.SetName("ContSelFinal");
+
+    TString mh3name = "abs(Alpha)";
+    MBinning binsalphaabs("Binning"+mh3name);
+    binsalphaabs.SetEdges(50, -2.0, 98.0);
+
+    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
+    alphaabs.SetName(mh3name);
+
+    TH1  &alphahist = alphaabs->GetHist();
+
+    MFillH alpha(&alphaabs);
+    alpha.SetName("FillAlphaAbs");
+
+    MFillH hfill1("MHHillas",    fHilName);
+    hfill1.SetName("HHillas");
+
+    MFillH hfill2("MHStarMap",   fHilName);
+    hfill2.SetName("HStarMap");
+
+    MFillH hfill3("MHHillasExt",    fHilNameSrc);
+    hfill3.SetName("HHillasExt");
+    
+    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
+    hfill4.SetName("HHillasSrc");    
+
+    MFillH hfill5("MHNewImagePar", fImgParName);
+    hfill5.SetName("HNewImagePar");
+
+    //*****************************
+    // entries in MParList
+
+    pliston.AddToList(&tliston);
+    InitBinnings(&pliston);
+
+    pliston.AddToList(&supercuts);
+
+    pliston.AddToList(&binsalphaabs);
+    pliston.AddToList(&alphaabs);
+
+    //*****************************
+    // entries in MTaskList
+    
+    tliston.AddToList(&read);
+
+    tliston.AddToList(&sccalc);
+    tliston.AddToList(&fillhadsc);
+
+    tliston.AddToList(&write);
+    tliston.AddToList(&contfinalgh);
+
+    tliston.AddToList(&alpha);
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+    tliston.AddToList(&hfill5);
+
+    tliston.AddToList(&contfinal);
+
+    //*****************************
+
+    //-------------------------------------------
+    // Execute event loop
+    //
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxevents = -1;
+    if ( !evtloop.Eventloop(maxevents) )
+        return;
+
+    tliston.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
+
+    pliston.FindObject("MHHillas")->DrawClone();
+    pliston.FindObject("MHHillasExt")->DrawClone();
+    pliston.FindObject("MHHillasSrc")->DrawClone();
+    pliston.FindObject("MHNewImagePar")->DrawClone();
+    pliston.FindObject("MHStarMap")->DrawClone();
+
+     //-------------------------------------------
+    // fit alpha distribution to get the number of excess events and
+    // calculate significance of gamma signal in the alpha plot
+  
+    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
+    TH1  &alphaHist = absalpha->GetHist();
+    alphaHist.SetXTitle("|alpha|  [\\circ]");
+    alphaHist.SetName("alpha-macro");
+
+    Double_t alphasig = 13.1;
+    Double_t alphamin = 30.0;
+    Double_t alphamax = 90.0;
+    Int_t    degree   =    2;
+    Double_t significance = -99.0;
+    Bool_t   drawpoly  = kTRUE;
+    Bool_t   fitgauss  = kTRUE;
+    Bool_t   print     = kTRUE;
+
+    MHFindSignificance findsig;
+    findsig.SetRebin(kTRUE);
+    findsig.SetReduceDegree(kFALSE);
+
+    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
+                        alphasig, drawpoly, fitgauss, print);
+    significance = findsig.GetSignificance();
+    Float_t alphasi = findsig.GetAlphasi();
+
+    gLog << "For file '" << filenameData << "' : " << endl;
+    gLog << "Significance of gamma signal after supercuts : "
+         << significance << " (for |alpha| < " << alphasi << " degrees)" 
+         << endl;
+
+    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
+
+    //-------------------------------------------
+
+    DeleteBinnings(&pliston);
+ }
+
+
+    gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl;
+    gLog << "=======================================================" << endl;
+ }
+  //---------------------------------------------------------------------
+
+
+
+  //---------------------------------------------------------------------
+  // Job C  
+  //======
+
+    //  - read ON1 and MC1 data files  
+    //    which should have been updated to contain the hadronnesses
+    //    for the method of Random Forest and for the SUPERCUTS
+    //  - produce Neyman-Pearson plots
+ 
+ if (JobC)
+ {
+    gLog << "=====================================================" << endl;
+    gLog << "Macro CT1Analysis : Start of Job C" << endl;
+
+    gLog << "" << endl;
+    gLog << "Macro CT1Analysis : JobC = " 
+         << (JobC ? "kTRUE" : "kFALSE")  << endl;
+
+
+    // name of input data file
+    TString filenameData = outPath;
+    filenameData += "OFF";
+    filenameData += "3.root";
+    gLog << "filenameData = " << filenameData << endl;
+
+    // name of input MC file
+    TString filenameMC = outPath;
+    filenameMC += "MC";
+    filenameMC += "3.root";
+    gLog << "filenameMC   = " << filenameMC   << endl;
+
+
+    //-----------------------------------------------------------------
+
+    MTaskList tliston;
+    MParList pliston;
+
+
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    //-------------------------------------------
+    // create the tasks which should be executed 
+    //
+
+    MReadMarsFile read("Events", filenameMC);
+    read.AddFile(filenameData);
+    read.DisableAutoScheme();
+
+
+    //.......................................................................
+    // names of hadronness containers
+
+    TString hadSCName = "HadSC";
+    TString hadRFName = "HadRF";
+
+    //.......................................................................
+
+
+    TString fHilName    = "MHillas"; 
+    TString fHilNameExt = "MHillasExt"; 
+    TString fHilNameSrc = "MHillasSrc"; 
+    TString fImgParName = "MNewImagePar"; 
+
+    Float_t maxhadronness =  0.40;
+    Float_t maxalpha      =  20.0;
+    Float_t maxdist       =  10.0;
+
+    MFCT1SelFinal selfinalgh(fHilNameSrc);
+    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
+    selfinalgh.SetHadronnessName(hadRFName);
+    selfinalgh.SetName("SelFinalgh");
+    MContinue contfinalgh(&selfinalgh);
+    contfinalgh.SetName("ContSelFinalgh");
+
+
+    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
+    fillhadsc.SetName("HhadSC");
+    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
+    fillhadrf.SetName("HhadRF");
+
+    MFCT1SelFinal selfinal(fHilNameSrc);
+    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
+    selfinal.SetHadronnessName(hadRFName);
+    selfinal.SetName("SelFinal");
+    MContinue contfinal(&selfinal);
+    contfinal.SetName("ContSelFinal");
+
+
+    MFillH hfill1("MHHillas",    fHilName);
+    hfill1.SetName("HHillas");
+
+    MFillH hfill2("MHStarMap",   fHilName);
+    hfill2.SetName("HStarMap");
+
+    MFillH hfill3("MHHillasExt",    fHilNameSrc);
+    hfill3.SetName("HHillasExt");
+    
+    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
+    hfill4.SetName("HHillasSrc");    
+
+    MFillH hfill5("MHNewImagePar", fImgParName);
+    hfill5.SetName("HNewImagePar");
+
+
+    //*****************************
+    // entries in MParList
+
+    pliston.AddToList(&tliston);
+    InitBinnings(&pliston);
+
+
+    //*****************************
+    // entries in MTaskList
+    
+    tliston.AddToList(&read);
+
+    tliston.AddToList(&fillhadsc);
+    tliston.AddToList(&fillhadrf);
+   
+    tliston.AddToList(&contfinalgh);
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+    tliston.AddToList(&hfill5);
+
+    tliston.AddToList(&contfinal);
+
+    //*****************************
+
+    //-------------------------------------------
+    // Execute event loop
+    //
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    evtloop.SetParList(&pliston);
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxevents = -1;
+    //Int_t maxevents = 35000;
+    if ( !evtloop.Eventloop(maxevents) )
+        return;
+
+    tliston.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+
+    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
+    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
+
+    pliston.FindObject("MHHillas")->DrawClone();
+    pliston.FindObject("MHHillasExt")->DrawClone();
+    pliston.FindObject("MHHillasSrc")->DrawClone();
+    pliston.FindObject("MHNewImagePar")->DrawClone();
+    pliston.FindObject("MHStarMap")->DrawClone();
+
+    DeleteBinnings(&pliston);
+
+    gLog << "Macro CT1Analysis : End of Job C" << endl;
+    gLog << "===================================================" << endl;
+ }
+
+
+
+  //---------------------------------------------------------------------
+  // Job D
+  //======
+
+    //  - select g/h separation method XX
+    //  - read ON2 (or MC2) root file 
+    //  - apply cuts in hadronness
+    //  - make plots
+
+
+ if (JobD)
+ {
+    gLog << "=====================================================" << endl;
+    gLog << "Macro CT1Analysis : Start of Job D" << endl;
+
+    gLog << "" << endl;
+    gLog << "Macro CT1Analysis : JobD = " 
+         << (JobD        ? "kTRUE" : "kFALSE")  << endl;
+
+
+    // type of data to be analysed
+    TString typeData = "ON";
+    //TString typeData = "OFF";
+    //TString typeData = "MC";
+    gLog << "typeData = " << typeData << endl;
+
+    TString ext      = "3.root";
+
+
+    //------------------------------
+    // selection of g/h separation method
+    // and definition of final selections
+
+    //TString XX("SC");
+    TString XX("RF");
+    TString fhadronnessName("Had");
+    fhadronnessName += XX;
+    gLog << "fhadronnessName = " << fhadronnessName << endl;
+
+    // maximum values of the hadronness, |ALPHA| and DIST
+    Float_t maxhadronness   = 0.233;
+    Float_t maxalpha        = 20.0;
+    Float_t maxdist         = 10.0;
+    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
+         << maxhadronness << ",  " << maxalpha << ",  " 
+         << maxdist << endl;
+
+
+    //------------------------------
+    // name of data file to be analysed
+    TString filenameData(outPath);
+    filenameData += typeData;
+    filenameData += ext;
+    gLog << "filenameData = " << filenameData << endl;
+
+
+
+    //*************************************************************************
+    //
+    // Analyse the data
+    //
+
+    MTaskList tliston;
+    MParList pliston;
+
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+
+    TString fHilName    = "MHillas"; 
+    TString fHilNameExt = "MHillasExt"; 
+    TString fHilNameSrc = "MHillasSrc"; 
+    TString fImgParName = "MNewImagePar"; 
+
+    //-------------------------------------------
+    // create the tasks which should be executed 
+    //
+
+    MReadMarsFile read("Events", filenameData);
+    read.DisableAutoScheme();
+
+
+    //-----------------------------------------------------------------
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
+    MFCT1SelFinal selfinalgh(fHilNameSrc);
+    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
+    selfinalgh.SetHadronnessName(fhadronnessName);
+    selfinalgh.SetName("SelFinalgh");
+    MContinue contfinalgh(&selfinalgh);
+    contfinalgh.SetName("ContSelFinalgh");
+
+    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
+    fillhadsc.SetName("HhadSC");
+    MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
+    fillhadrf.SetName("HhadRF");
+
+    TString mh3name = "abs(Alpha)";
+    MBinning binsalphaabs("Binning"+mh3name);
+    binsalphaabs.SetEdges(50, -2.0, 98.0);
+
+    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
+    alphaabs.SetName(mh3name);
+
+    TH1  &alphahist = alphaabs->GetHist();
+
+    MFillH alpha(&alphaabs);
+    alpha.SetName("FillAlphaAbs");
+
+    MFillH hfill1("MHHillas",    fHilName);
+    hfill1.SetName("HHillas");
+
+    MFillH hfill2("MHStarMap",   fHilName);
+    hfill2.SetName("HStarMap");
+
+    MFillH hfill3("MHHillasExt",   fHilNameSrc);
+    hfill3.SetName("HHillasExt");    
+
+    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
+    hfill4.SetName("HHillasSrc");    
+
+    MFillH hfill5("MHNewImagePar", fImgParName);
+    hfill5.SetName("HNewImagePar");    
 
     MFCT1SelFinal selfinal(fHilNameSrc);
     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
     selfinal.SetHadronnessName(fhadronnessName);
+    selfinal.SetName("SelFinal");
     MContinue contfinal(&selfinal);
+    contfinal.SetName("ContSelFinal");
 
 
@@ -1968 +2295 @@
     // entries in MParList
 
-    plistoff.AddToList(&tlistoff);
-    InitBinnings(&plistoff);
-
+    pliston.AddToList(&tliston);
+    InitBinnings(&pliston);
+    pliston.AddToList(&binsalphaabs);
+    pliston.AddToList(&alphaabs);
 
     //*****************************
     // entries in MTaskList
     
-    tlistoff.AddToList(&read);
-    tlistoff.AddToList(&eest2);
-    tlistoff.AddToList(&write);
-    tlistoff.AddToList(&contfinal);
+    tliston.AddToList(&read);
+
+    tliston.AddToList(&contfinalgh);
+
+    tliston.AddToList(&fillhadsc);
+    tliston.AddToList(&fillhadrf);
+
+    tliston.AddToList(&alpha);
+    tliston.AddToList(&hfill1);
+    tliston.AddToList(&hfill2);
+    tliston.AddToList(&hfill3);
+    tliston.AddToList(&hfill4);
+    tliston.AddToList(&hfill5);
+
+    tliston.AddToList(&contfinal);
 
     //*****************************
@@ -1987 +2326 @@
     MProgressBar bar;
     MEvtLoop evtloop;
-    evtloop.SetParList(&plistoff);
+    evtloop.SetParList(&pliston);
     evtloop.SetProgressBar(&bar);
 
     Int_t maxevents = -1;
-    //Int_t maxevents = 1000;
+    //Int_t maxevents = 10000;
     if ( !evtloop.Eventloop(maxevents) )
         return;
 
-    tlistoff.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&plistoff);
-
-    //---------------------------------------------------
-
-    //---------------------------------------------------
-    // Update ON data
+    tliston.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // Display the histograms
+    //
+
+    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
+    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
+
+    pliston.FindObject("MHHillas")->DrawClone();
+    pliston.FindObject("MHHillasExt")->DrawClone();
+    pliston.FindObject("MHHillasSrc")->DrawClone();
+    pliston.FindObject("MHNewImagePar")->DrawClone();
+    pliston.FindObject("MHStarMap")->DrawClone();
+
+
+    //-------------------------------------------
+
+    // fit alpha distribution to get the number of excess events and
+    // calculate significance of gamma signal in the alpha plot
+  
+    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
+    TH1  &alphaHist = absalpha->GetHist();
+    alphaHist.SetXTitle("|alpha|  [\\circ]");
+    alphaHist.SetName("alpha-JobD");
+
+    Double_t alphasig = 13.1;
+    Double_t alphamin = 30.0;
+    Double_t alphamax = 90.0;
+    Int_t    degree   =    2;
+    Double_t significance = -99.0;
+    Bool_t   drawpoly  = kTRUE;
+    Bool_t   fitgauss  = kTRUE;
+    Bool_t   print     = kTRUE;
+
+    MHFindSignificance findsig;
+    findsig.SetRebin(kTRUE);
+    findsig.SetReduceDegree(kFALSE);
+
+    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
+                        alphasig, drawpoly, fitgauss, print);
+    significance = findsig.GetSignificance();
+    Float_t alphasi = findsig.GetAlphasi();
+
+    gLog << "For file '" << filenameData << "' : " << endl;
+    gLog << "Significance of gamma signal after supercuts : "
+         << significance << " (for |alpha| < " << alphasi << " degrees)" 
+         << endl;
+
+    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
+
+    //-------------------------------------------
+
+
+    DeleteBinnings(&pliston);
+
+    gLog << "Macro CT1Analysis : End of Job D" << endl;
+    gLog << "=======================================================" << endl;
+ }
+  //---------------------------------------------------------------------
+
+
+
+
+
+  //---------------------------------------------------------------------
+  // Job E_XX
+  //=========
+
+    //  - select g/h separation method XX
+    //  - read MC_XX2.root file 
+    //  - calculate eff. collection area
+    //  - read ON_XX2.root file 
+    //  - apply final cuts
+    //  - calculate flux
+    //  - write root file for ON data after final cuts (ON_XX3.root))
+
+
+ if (JobE_XX)
+ {
+    gLog << "=====================================================" << endl;
+    gLog << "Macro CT1Analysis : Start of Job E_XX" << endl;
+
+    gLog << "" << endl;
+    gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = " 
+         << (JobE_XX ? "kTRUE" : "kFALSE")  << ",  " 
+         << (OEEst ?   "kTRUE" : "kFALSE")  << ",  " 
+         << (WEX     ? "kTRUE" : "kFALSE")  << endl;
+
+
+    // type of data to be analysed
+    //TString typeData = "ON";
+    //TString typeData = "OFF";
+    TString typeData = "MC";
+    gLog << "typeData = " << typeData << endl;
+
+    TString typeMC   = "MC";
+    TString ext      = "3.root";
+    TString extout   = "4.root";
+
+    //------------------------------
+    // selection of g/h separation method
+    // and definition of final selections
+
+    //TString XX("SC");
+    TString XX("RF");
+    TString fhadronnessName("Had");
+    fhadronnessName += XX;
+    gLog << "fhadronnessName = " << fhadronnessName << endl;
+
+    // maximum values of the hadronness, |ALPHA| and DIST
+    Float_t maxhadronness   = 0.23;
+    Float_t maxalpha        = 20.0;
+    Float_t maxdist         = 10.0;
+    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
+         << maxhadronness << ",  " << maxalpha << ",  " 
+         << maxdist << endl;
+
+    //------------------------------
+    // name of MC file to be used for optimizing the energy estimator
+    TString filenameOpt(outPath);
+    filenameOpt += typeMC;
+    filenameOpt += ext; 
+    gLog << "filenameOpt = " << filenameOpt << endl;
+
+    //------------------------------
+    // name of file containing the parameters of the energy estimator
+    TString energyParName(outPath);
+    energyParName += "energyest_";
+    energyParName += XX;
+    energyParName += ".root";
+    gLog << "energyParName = " << energyParName << endl;
+
+    //------------------------------
+    // name of MC file to be used for calculating the eff. collection areas
+    TString filenameArea(outPath);
+    filenameArea += typeMC;
+    filenameArea += ext; 
+    gLog << "filenameArea = " << filenameArea << endl;
+
+    //------------------------------
+    // name of file containing the eff. collection areas
+    TString collareaName(outPath);
+    collareaName += "area_";
+    collareaName += XX;
+    collareaName += ".root";
+    gLog << "collareaName = " << collareaName << endl;
+
+    //------------------------------
+    // name of data file to be analysed
+    TString filenameData(outPath);
+    filenameData += typeData;
+    filenameData += ext;
+    gLog << "filenameData = " << filenameData << endl;
+
+    //------------------------------
+    // name of output data file (after the final cuts)
+    TString filenameDataout(outPath);
+    filenameDataout += typeData;
+    filenameDataout += "_";
+    filenameDataout += XX;
+    filenameDataout += extout;
+    gLog << "filenameDataout = " << filenameDataout << endl;
+
+    //------------------------------
+    // name of file containing histograms for flux calculastion
+    TString filenameResults(outPath);
+    filenameResults += typeData;
+    filenameResults += "Results_";
+    filenameResults += XX;
+    filenameResults += extout;
+    gLog << "filenameResults = " << filenameResults << endl;
+
+
+    //====================================================================
+    gLog << "-----------------------------------------------" << endl;
+    gLog << "Start calculation of effective collection areas" << endl;
+    MParList  parlist;
+    MTaskList tasklist;
+
+    //---------------------------------------
+    // Setup the tasks to be executed
+    //
+    MReadMarsFile reader("Events", filenameArea);
+    reader.DisableAutoScheme();
+
+    MFCT1SelFinal cuthadrons;
+    cuthadrons.SetHadronnessName(fhadronnessName);
+    cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
+
+    MContinue conthadrons(&cuthadrons);
+
+    MHMcCT1CollectionArea collarea;
+    collarea.SetEaxis(MHMcCT1CollectionArea::kLinear);
+
+    MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
+    filler.SetName("CollectionArea");
+
+    //********************************
+    // entries in MParList
+
+    parlist.AddToList(&tasklist);
+    InitBinnings(&parlist);
+    parlist.AddToList(&collarea);
+
+    //********************************
+    // entries in MTaskList
+
+    tasklist.AddToList(&reader);   
+    tasklist.AddToList(&conthadrons);
+    tasklist.AddToList(&filler);
+
+    //********************************
+
+    //-----------------------------------------
+    // Execute event loop
+    //
+    MEvtLoop magic;
+    magic.SetParList(&parlist);
+
+    MProgressBar bar;
+    magic.SetProgressBar(&bar);
+    if (!magic.Eventloop())
+        return;
+
+    tasklist.PrintStatistics(0, kTRUE);
+
+    // Calculate effective collection areas 
+    // and display the histograms
+    //
+    //MHMcCT1CollectionArea *collarea = 
+    //     (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
+    collarea.CalcEfficiency();
+    collarea.DrawClone();
+
+    // save binnings for call to CT1EEst
+    MBinning *binsE     = (MBinning*)parlist.FindObject("BinningE");
+    if (!binsE)
+        {
+          gLog << "Object 'BinningE' not found in MParList" << endl;
+          return;
+        }
+    MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
+    if (!binsTheta)
+        {
+          gLog << "Object 'BinningTheta' not found in MParList" << endl;
+          return;
+        }
+
+
+    //---------------------------------------------
+    // Write histograms to a file 
+    //
+
+    TFile f(collareaName, "RECREATE");
+    collarea.GetHist()->Write();
+    collarea.GetHAll()->Write();
+    collarea.GetHSel()->Write();
+    f.Close();
+
+    gLog << "Collection area plots written onto file " << collareaName << endl;
+
+    gLog << "Calculation of effective collection areas done" << endl;
+    gLog << "-----------------------------------------------" << endl;    
+    //------------------------------------------------------------------
+
+
+    TString fHilName    = "MHillas"; 
+    TString fHilNameExt = "MHillasExt"; 
+    TString fHilNameSrc = "MHillasSrc"; 
+    TString fImgParName = "MNewImagePar"; 
+
+
+ if (OEEst)
+ { 
+   //===========================================================
+    //
+    // Optimization of energy estimator
+    //
+    gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl;
+
+    TString inpath("");
+    TString outpath("");
+    Int_t howMany = 2000;
+    CT1EEst(inpath,   filenameOpt,   outpath, energyParName, 
+            fHilName, fHilNameSrc,   fhadronnessName,
+            howMany,  maxhadronness, maxalpha, maxdist,
+            binsE, binsTheta);
+    gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
+ }
+
+ if (WEX)
+ {
+    //-----------------------------------------------------------
+    //
+    // Read in parameters of energy estimator ("MMcEnergyEst")
+    //                   and migration matrix ("MHMcEnergyMigration")
+    //
+    gLog << "================================================================"
+         << endl;
+    gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
+         << energyParName << "'" << endl;
+    TFile enparam(energyParName);
+    enparam.ls();
+    MMcEnergyEst mcest("MMcEnergyEst"); 
+    mcest.Read("MMcEnergyEst");
+
+    //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
+    gLog << "Parameters of energy estimator were read in" << endl;
+
+
+    gLog << "Read in Migration matrix" << endl;   
+
+    MHMcEnergyMigration mighiston("MHMcEnergyMigration");
+    mighiston.Read("MHMcEnergyMigration");
+    //MHMcEnergyMigration &mighiston = 
+    //      *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
+
+    gLog << "Migration matrix was read in" << endl;
+
+
+    TArrayD parA(mcest.GetNumCoeffA());
+    TArrayD parB(mcest.GetNumCoeffB());
+    for (Int_t i=0; i<parA.GetSize(); i++)
+      parA[i] = mcest.GetCoeff(i);
+    for (Int_t i=0; i<parB.GetSize(); i++)
+      parB[i] = mcest.GetCoeff( i+parA.GetSize() );
+
+    //*************************************************************************
+    //
+    // Analyse the data
     //
     gLog << "============================================================"
          << endl;
-    gLog << "Macro CT1Analysis.C : update file '" << filenameON
-         << "'" << endl;
+    gLog << "Analyse the data" << endl;
 
     MTaskList tliston;
     MParList pliston;
-
 
     // geometry is needed in  MHHillas... classes 
@@ -2016 +2678 @@
              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
 
+
     //-------------------------------------------
     // create the tasks which should be executed 
     //
 
-    MReadMarsFile read("Events", filenameON);
+    MReadMarsFile read("Events", filenameData);
     read.DisableAutoScheme();
 
-    //---------------------------
-    // calculate estimated energy
-
-    MEnergyEstParam eest2(fHilName);
-    eest2.Add(fHilNameSrc);
-
-    eest2.SetCoeffA(parA);
-    eest2.SetCoeffB(parB);
-
     //.......................................................................
 
-      MWriteRootFile write(filenameONup);
+      gLog << "CT1Analysis.C : write root file '" << filenameDataout 
+           << "'" << endl;
+   
+      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+
+
+      MWriteRootFile write(filenameDataout, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -2053 +2713 @@
       write.AddContainer("MEnergyEst",    "Events");
 
-    //-----------------------------------------------------------------
-
-    MFCT1SelFinal selfinal(fHilNameSrc);
-    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
-    selfinal.SetHadronnessName(fhadronnessName);
-    MContinue contfinal(&selfinal);
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-    tliston.AddToList(&eest2);
-    tliston.AddToList(&write);
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-
-    Int_t maxevents = -1;
-    //Int_t maxevents = 1000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&pliston);
-
-    //---------------------------------------------------
-
-    //---------------------------------------------------
-    // Update MC data
-    //
-    gLog << "============================================================"
-         << endl;
-    gLog << "Macro CT1Analysis.C : update file '" << filenameMC
-         << "'" << endl;
-
-    MTaskList tlistmc;
-    MParList plistmc;
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameMC);
-    read.DisableAutoScheme();
-
-    //---------------------------
-    // calculate estimated energy
-
-    MEnergyEstParam eest2(fHilName);
-    eest2.Add(fHilNameSrc);
-
-    eest2.SetCoeffA(parA);
-    eest2.SetCoeffB(parB);
-
-    //.......................................................................
-
-      MWriteRootFile write(filenameMCup);
-
-      write.AddContainer("MRawRunHeader", "RunHeaders");
-      write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      write.AddContainer("ThetaOrig",     "Events");
-      write.AddContainer("MSrcPosCam",    "Events");
-      write.AddContainer("MSigmabar",     "Events");
-      write.AddContainer("MHillas",       "Events");
-      write.AddContainer("MHillasExt",    "Events");
-      write.AddContainer("MHillasSrc",    "Events");
-      write.AddContainer("MNewImagePar",  "Events");
-
-      //write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
-      write.AddContainer("HadRF",         "Events");
-
-      write.AddContainer("MEnergyEst",    "Events");
-
-    //-----------------------------------------------------------------
-
-    MFCT1SelFinal selfinal(fHilNameSrc);
-    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
-    selfinal.SetHadronnessName(fhadronnessName);
-    MContinue contfinal(&selfinal);
-
-
-    //*****************************
-    // entries in MParList
-
-    plistmc.AddToList(&tlistmc);
-    InitBinnings(&plistmc);
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tlistmc.AddToList(&read);
-    tlistmc.AddToList(&eest2);
-    tlistmc.AddToList(&write);
-    tlistmc.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&plistmc);
-    evtloop.SetProgressBar(&bar);
-
-    Int_t maxevents = -1;
-    //Int_t maxevents = 1000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tlistmc.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&plistmc);
-
-
-    //==========   end update   ============================================
-   }
-    
-    enparam.Close();
-
-    gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl;
-    gLog << "=======================================================" << endl;
- }
-  //---------------------------------------------------------------------
-
-
-  //---------------------------------------------------------------------
-  // Job F_XX
-  //=========
-
-    //  - select g/h separation method XX
-    //  - read MC_XX2.root file 
-    //  - calculate eff. collection area
-    //  - read ON_XX2.root file 
-    //  - apply final cuts
-    //  - calculate flux
-    //  - write root file for ON data after final cuts (ON_XX3.root))
-
-
- if (JobF_XX)
- {
-    gLog << "=====================================================" << endl;
-    gLog << "Macro CT1Analysis : Start of Job F_XX" << endl;
-
-    gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobF_XX, WXX = " 
-         << JobF_XX  << ",  " << WXX << endl;
-
-    // type of data to be analysed
-    //TString typeData = "ON";
-    //TString typeData = "OFF";
-    TString typeData = "MC";
-    gLog << "typeData = " << typeData << endl;
-
-    TString typeMC   = "MC";
-    TString ext      = "3.root";
-    TString extout   = "4.root";
-
-    //------------------------------
-    // selection of g/h separation method
-    // and definition of final selections
-
-    //TString XX("NN");
-    //TString XX("SC");
-    TString XX("RF");
-    TString fhadronnessName("Had");
-    fhadronnessName += XX;
-    gLog << "fhadronnessName = " << fhadronnessName << endl;
-
-    // maximum values of the hadronness, |ALPHA| and DIST
-    Float_t maxhadronness   = 0.40;
-    Float_t maxalpha        = 20.0;
-    Float_t maxdist         = 10.0;
-    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
-         << maxhadronness << ",  " << maxalpha << ",  " 
-         << maxdist << endl;
-
-
-    //------------------------------
-    // name of MC file to be used for calculating the eff. collection areas
-    TString filenameArea(outPath);
-    filenameArea += typeMC;
-    filenameArea += "_";
-    filenameArea += XX;
-    filenameArea += ext; 
-    gLog << "filenameArea = " << filenameArea << endl;
-
-    //------------------------------
-    // name of file containing the eff. collection areas
-    TString collareaName(outPath);
-    collareaName += "area_";
-    collareaName += XX;
-    collareaName += ".root";
-    gLog << "collareaName = " << collareaName << endl;
-
-    //------------------------------
-    // name of data file to be analysed
-    TString filenameData(outPath);
-    filenameData += typeData;
-    filenameData += "_";
-    filenameData += XX;
-    filenameData += ext;
-    gLog << "filenameData = " << filenameData << endl;
-
-    //------------------------------
-    // name of output data file (after the final cuts)
-    TString filenameDataout(outPath);
-    filenameDataout += typeData;
-    filenameDataout += "_";
-    filenameDataout += XX;
-    filenameDataout += extout;
-    gLog << "filenameDataout = " << filenameDataout << endl;
-
-
-    //====================================================================
-    gLog << "-----------------------------------------------" << endl;
-    gLog << "Start calculation of effective collection areas" << endl;
-    MParList  parlist;
-    MTaskList tasklist;
-
-    //---------------------------------------
-    // Setup the tasks to be executed
-    //
-    MReadMarsFile reader("Events", filenameArea);
-    reader.DisableAutoScheme();
-
-    MFCT1SelFinal cuthadrons;
-    cuthadrons.SetHadronnessName(fhadronnessName);
-    cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
-
-    MContinue conthadrons(&cuthadrons);
-
-    //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
-    //MHMcCT1CollectionArea* collarea;
-
-    MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
-    filler.SetName("CollectionArea");
-
-    //********************************
-    // entries in MParList
-
-    parlist.AddToList(&tasklist);
-    InitBinnings(&parlist);
-    //parlist.AddToList(collarea);
-
-    //********************************
-    // entries in MTaskList
-
-    tasklist.AddToList(&reader);   
-    tasklist.AddToList(&conthadrons);
-    tasklist.AddToList(&filler);
-
-    //********************************
-
-    //-----------------------------------------
-    // Execute event loop
-    //
-    MEvtLoop magic;
-    magic.SetParList(&parlist);
-
-    MProgressBar bar;
-    magic.SetProgressBar(&bar);
-    if (!magic.Eventloop())
-        return;
-
-    tasklist.PrintStatistics(0, kTRUE);
-
-    // Calculate effective collection areas 
-    // and display the histograms
-    //
-
-    MHMcCT1CollectionArea *collarea = 
-        (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
-
-    collarea->CalcEfficiency();
-    collarea->DrawClone("lego");
-
-    //---------------------------------------------
-    // Write histograms to a file 
-    //
-
-    TFile f(collareaName, "RECREATE");
-    collarea->GetHist()->Write();
-    collarea->GetHAll()->Write();
-    collarea->GetHSel()->Write();
-    f.Close();
-
-    //delete collarea;
-
-    gLog << "Calculation of effective collection areas done" << endl;
-    gLog << "-----------------------------------------------" << endl;    
-    //------------------------------------------------------------------
-
-
-    //*************************************************************************
-    //
-    // Analyse the data
-    //
-
-    MTaskList tliston;
-    MParList pliston;
-
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-    TString fHilName    = "MHillas"; 
-    TString fHilNameExt = "MHillasExt"; 
-    TString fHilNameSrc = "MHillasSrc"; 
-    TString fImgParName = "MNewImagePar"; 
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameData);
-    read.DisableAutoScheme();
-
-    //.......................................................................
-
-
-      MWriteRootFile write(filenameDataout);
-
-      write.AddContainer("MRawRunHeader", "RunHeaders");
-      write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      write.AddContainer("ThetaOrig",     "Events");
-      write.AddContainer("MSrcPosCam",    "Events");
-      write.AddContainer("MSigmabar",     "Events");
-      write.AddContainer("MHillas",       "Events");
-      write.AddContainer("MHillasExt",    "Events");
-      write.AddContainer("MHillasSrc",    "Events");
-      write.AddContainer("MNewImagePar",  "Events");
-
-      //write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
-      write.AddContainer("HadRF",         "Events");
-
-      write.AddContainer("MEnergyEst",    "Events");
-
 
     //-----------------------------------------------------------------
@@ -2431 +2733 @@
     fillhadrf.SetName("HhadRF");
 
+    //---------------------------
+    // calculate estimated energy
+
+    MEnergyEstParam eeston(fHilName);
+    eeston.Add(fHilNameSrc);
+
+    eeston.SetCoeffA(parA);
+    eeston.SetCoeffB(parB);
+
+    //---------------------------
+    // calculate estimated energy using Daniel's parameters
+
+    //MEnergyEstParamDanielMkn421 eeston(fHilName);
+    //eeston.Add(fHilNameSrc);
+    //eeston.SetCoeffA(parA);
+    //eeston.SetCoeffB(parB);
+
+
+    //---------------------------
+
 
     MFillH hfill1("MHHillas",    fHilName);
@@ -2446 +2768 @@
     MFillH hfill5("MHNewImagePar", fImgParName);
     hfill5.SetName("HNewImagePar");    
+
+    //---------------------------
+    // new from Robert
+
+    MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
+    hfill6.SetName("HTimeDiffTheta");
+
+    MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
+    hfill6a.SetName("HTimeDiffTime");
+
+    MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
+    hfill7.SetName("HAlphaEnergyTheta");
+
+    MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
+    hfill7a.SetName("HAlphaEnergyTime");
+
+    MFillH hfill7b("MHThetabarTime", fHilNameSrc);
+    hfill7b.SetName("HThetabarTime");
+
+    MFillH hfill7c("MHEnergyTime", "MMcEvt");
+    hfill7c.SetName("HEnergyTime");
+
+
+    //---------------------------
 
     MFCT1SelFinal selfinal(fHilNameSrc);
@@ -2467 +2813 @@
     tliston.AddToList(&read);
 
+    // robert      
+    tliston.AddToList(&hfill6);   //timediff
+    tliston.AddToList(&hfill6a);   //timediff
+
     tliston.AddToList(&contfinalgh);
-
-    if (WXX)
-      tliston.AddToList(&write);
+    tliston.AddToList(&eeston);
+
+    tliston.AddToList(&write);
 
     //tliston.AddToList(&fillhadnn);
@@ -2482 +2832 @@
     tliston.AddToList(&hfill5);
 
+    //robert
+    tliston.AddToList(&hfill7);
+    tliston.AddToList(&hfill7a);
+    tliston.AddToList(&hfill7b);
+    tliston.AddToList(&hfill7c);
+
     tliston.AddToList(&contfinal);
 
@@ -2495 +2851 @@
 
     Int_t maxevents = -1;
-    //Int_t maxevents = 10000;
     if ( !evtloop.Eventloop(maxevents) )
         return;
@@ -2507 +2862 @@
 
     //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
+
+    gLog << "before hadRF" << endl;
     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
+
+    gLog << "before hadSC" << endl;
     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
 
+    gLog << "before MHHillas" << endl;
     pliston.FindObject("MHHillas")->DrawClone();
+
+    gLog << "before MHHillasExt" << endl;
     pliston.FindObject("MHHillasExt")->DrawClone();
+
+    gLog << "before MHHillasSrc" << endl;
     pliston.FindObject("MHHillasSrc")->DrawClone();
+
+    gLog << "before MHNewImagePar" << endl;
     pliston.FindObject("MHNewImagePar")->DrawClone();
+
+    gLog << "before MHStarMap" << endl;
     pliston.FindObject("MHStarMap")->DrawClone();
 
+    gLog << "before DeleteBinnings" << endl;
+
     DeleteBinnings(&pliston);
 
-    gLog << "Macro CT1Analysis : End of Job F_XX" << endl;
+    gLog << "before Robert's code" << endl;
+
+
+//rwagner write all relevant histograms onto a file
+
+  if (WRobert)
+  {
+    gLog << "=======================================================" << endl;
+    gLog << "Write results onto file '" << filenameResults << "'" << endl;
+
+    TFile outfile(filenameResults,"recreate");
+
+    MHHillasSrc* hillasSrc = 
+      (MHHillasSrc*)(pliston->FindObject("MHHillasSrc"));
+        TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
+    alphaHist->Write();
+    gLog << "Alpha plot has been written out" << endl;    
+
+
+    MHAlphaEnergyTheta* aetH = 
+      (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
+    TH3D* aetHist = (TH3D*)(aetH->GetHist());
+    aetHist->SetName("aetHist");
+    aetHist->Write();
+    gLog << "AlphaEnergyTheta plot has been written out" << endl;    
+
+    MHAlphaEnergyTime* aetH2 = 
+      (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime"));
+    TH3D* aetHist2 = (TH3D*)(aetH2->GetHist());
+    aetHist2->SetName("aetimeHist");
+//     aetHist2->DrawClone();
+    aetHist2->Write();
+    gLog << "AlphaEnergyTime plot has been written out" << endl;    
+
+    MHThetabarTime* tbt = 
+      (MHThetabarTime*)(pliston->FindObject("MHThetabarTime"));
+    TProfile* tbtHist = (TProfile*)(tbt->GetHist());
+    tbtHist->SetName("tbtHist");
+    tbtHist->Write();
+    gLog << "ThetabarTime plot has been written out" << endl;    
+
+    MHEnergyTime* ent = 
+      (MHEnergyTime*)(pliston->FindObject("MHEnergyTime"));
+    TH2D* entHist = (TH2D*)(ent->GetHist());
+    entHist->SetName("entHist");
+    entHist->Write();
+    gLog << "EnergyTime plot has been written out" << endl;    
+    
+    MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
+    TH2D* timeHist = (TH2D*)(time->GetHist());
+    timeHist->SetName("MHTimeDiffTheta");
+    timeHist->SetTitle("Time diffs");
+    timeHist->Write();
+    gLog << "TimeDiffTheta plot has been written out" << endl;    
+
+
+    MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
+    TH2D* timeHist2 = (TH2D*)(time2->GetHist());
+    timeHist2->SetName("MHTimeDiffTime");
+    timeHist2->SetTitle("Time diffs");
+    timeHist2->Write();
+    gLog << "TimeDiffTime plot has been written out" << endl;    
+
+//rwagner write also collareas to same file
+    collarea->GetHist()->Write();
+    collarea->GetHAll()->Write();
+    collarea->GetHSel()->Write();
+    gLog << "Effective collection areas have been written out" << endl;        
+
+//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data
+//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file
+
+    gLog << "before closing outfile" << endl;
+
+    //outfile.Close();
+    gLog << "Results were written onto file '" << filenameResults 
+         << "'" << endl;
+    gLog << "=======================================================" << endl;
+  }
+
+  }
+
+    gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
     gLog << "=======================================================" << endl;
  }
@@ -2527 +2979 @@
 
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+