Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C
===================================================================
--- trunk/MagicSoft/Mars/macros/CT1Analysis.C	(revision 2317)
+++ trunk/MagicSoft/Mars/macros/CT1Analysis.C	(revision 2357)
@@ -22,4 +22,5 @@
 #include "MWriteRootFile.h"
 
+//#include "TH3D.h"
 
 
@@ -77,12 +78,12 @@
 
         MBinning *binth = new MBinning("BinningTheta");
-        //Double_t yedge[9] = 
-        //               {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
         // 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};
-        //TArrayD yed(9,yedge);
+        //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;
-        yed.Set(8,yedge);
+        yed.Set(9,yedge);
         binth->SetEdges(yed);
         plist->AddToList(binth);
@@ -94,5 +95,6 @@
           zedge[8-i] = cos(yedge[i]/kRad2Deg);
 	}
-        TArrayD zed(9,zedge);
+        TArrayD zed;
+        zed.Set(9,zedge);
         bincosth->SetEdges(zed);
         plist->AddToList(bincosth);
@@ -116,4 +118,5 @@
 }
 
+
 void DeleteBinnings(MParList *plist)
 {
@@ -122,4 +125,10 @@
         TObject *bin;
 
+        //--------------------------------------------
+        bin = plist->FindObject("BinningE");
+        if (bin) delete bin;
+
+        //--------------------------------------------
+
         bin = plist->FindObject("BinningSize");
         if (bin) delete bin;
@@ -159,5 +168,6 @@
         if (bin) delete bin;
 
-        //robert
+
+        // robert ----------------------------------------------
         bin = plist->FindObject("BinningAlphaFlux");
         if (bin) delete bin;
@@ -170,7 +180,11 @@
 }
 
+
 //************************************************************************
 void CT1Analysis()
 {
+
+  gLog << "Entry CT1Analysis()" << endl;
+
       gLog.SetNoColors();
 
@@ -189,8 +203,8 @@
 
       // path for input for Mars
-      TString inPath = "~magican/ct1test/wittek/marsoutput/optionB/";
+      TString inPath = "~magican/ct1test/wittek/marsoutput/optionC/";
 
       // path for output from Mars
-      TString outPath = "~magican/ct1test/wittek/marsoutput/optionB/";
+      TString outPath = "~magican/ct1test/wittek/marsoutput/optionC/";
 
       //-----------------------------------------------
@@ -217,27 +231,4 @@
     Bool_t JobA_MC  = kFALSE;  
     Bool_t WMC1     = kFALSE;  // write out root file MC1.root ?
-
-
-    // Job B_NN_UP : read ON1.root (or MC1.root) file 
-    //  - depending on RlookNN : create (or read in) hadron and gamma matrix
-    //  - calculate hadroness for method of NEAREST NEIGHBORS
-    //    and for the SUPERCUTS
-    //  - update the input files with the hadronesses (ON1.root or MC1.root)
-
-    Bool_t JobB_NN_UP  = kFALSE;  
-    Bool_t RLookNN     = kFALSE;  // read in look-alike events
-    Bool_t WNN         = kFALSE;  // update input root file ?
-
-
-    // Job B_SC_UP : read ON1.root (or MC1.root) file 
-    //  - depending on WParSC : create (or read in) supercuts parameter values
-    //  - calculate hadroness for the SUPERCUTS
-    //  - update the input files with the hadroness (ON1.root or MC1.root)
-
-    Bool_t JobB_SC_UP  = kTRUE;
-    Bool_t RMatrix     = kFALSE;  // read matrices from file  
-    Bool_t WParSC      = kTRUE;  // do optimization and write supercuts 
-                                  // parameter values onto a file
-    Bool_t WSC         = kFALSE;  // update input root file ?
 
 
@@ -252,4 +243,18 @@
     Bool_t RTree       = kFALSE;  // read in trees
     Bool_t WRF         = kFALSE;  // update input root file ?
+
+
+
+
+    // Job B_SC_UP : read ON1.root (or MC1.root) file 
+    //  - depending on WParSC : create (or read in) supercuts parameter values
+    //  - calculate hadroness for the SUPERCUTS
+    //  - update the input files with the hadroness (ON1.root or MC1.root)
+
+    Bool_t JobB_SC_UP  = kFALSE;
+    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file  
+    Bool_t WParSC      = kFALSE;  // do optimization and write supercuts 
+                                 // parameter values onto a file
+    Bool_t WSC         = kFALSE; // update input root file ?
 
 
@@ -299,6 +304,9 @@
     //  - write root file for ON data after final cuts 
 
-    Bool_t JobF_XX  = kFALSE;  
-    Bool_t WFX      = kFALSE;  // write out root file  ?
+    gLog << "before setting switches for JobF_XX" << endl;
+
+    Bool_t JobF_XX  = kTRUE;  
+    Bool_t WFX      = kTRUE;  // write out root file  ?
+    Bool_t WRobert  = kTRUE;  // write out Robert's file  ?
 
 
@@ -327,5 +335,7 @@
     gLog << "" << endl;
     gLog << "Macro CT1Analysis : JobA_ON, WHistON, WON1 = " 
-         << JobA_ON  << ",  " << WHistON << ",  " << WON1 << endl;
+         << (JobA_ON ? "kTRUE" : "kFALSE")  << ",  " 
+         << (WHistON ? "kTRUE" : "kFALSE")  << ",  " 
+         << (WON1    ? "kTRUE" : "kFALSE")  << endl;
 
 
@@ -424,7 +434,6 @@
     contstandard.SetName("SelStandard");
 
-    if (WON1)
-    {
-      MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
+      //MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
+      MWriteRootFile write(outNameImage, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -438,5 +447,5 @@
       write.AddContainer("MHillasSrc",    "Events");
       write.AddContainer("MNewImagePar",  "Events");
-    }
+
 
     //$$$$$$$$$$$$$$$$$$$$$$$$$$$$
@@ -500,5 +509,4 @@
 
     Int_t maxevents = -1;
-    //Int_t maxevents = 1000;
     if ( !evtloop.Eventloop(maxevents) )
         return;
@@ -587,5 +595,6 @@
     gLog << "" << endl;
     gLog << "Macro CT1Analysis : JobA_MC, WMC1 = " 
-         << JobA_MC  << ",  " << WMC1 << endl;
+         << (JobA_MC ? "kTRUE" : "kFALSE")  << ",  " 
+         << (WMC1    ? "kTRUE" : "kFALSE")  << endl;
 
 
@@ -688,6 +697,11 @@
     MCT1ReadPreProc read(filenamein);
 
+    MBlindPixelCalc blindbeforepad;
+    blindbeforepad.SetUseBlindPixels();
+    blindbeforepad.SetName("BlindBeforePadding");
+
     MBlindPixelCalc blind;
     blind.SetUseBlindPixels();
+    blind.SetName("BlindAfterPadding");
 
     MFCT1SelBasic selbasic;
@@ -758,7 +772,6 @@
 
 
-    if (WMC1)
-    {
-      MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
+      //MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
+      MWriteRootFile write(outNameImage, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -772,5 +785,5 @@
       write.AddContainer("MHillasSrc",    "Events");
       write.AddContainer("MNewImagePar",  "Events");
-    }
+
 
 
@@ -788,4 +801,5 @@
 
     tlist.AddToList(&read);
+    tlist.AddToList(&blindbeforepad);
     tlist.AddToList(&padthomas);
     tlist.AddToList(&blind);
@@ -855,28 +869,35 @@
  }
 
-
-
   //---------------------------------------------------------------------
-  // Job B_NN_UP
+  // Job B_RF_UP
   //============
 
-    //  - read in the matrices of look-alike events for gammas and hadrons
-
-    // then read ON1.root (or MC1.root) file 
-    //
-    //  - calculate the hadroness for the method of nearest neighbors
-    //
-    //  - update input root file, including the hadroness
-
-
- if (JobB_NN_UP)
+
+    //  - create (or read in) the matrices of look-alike events for gammas 
+    //    and hadrons
+    //  - create (or read in) the trees
+    //  - then read ON1.root (or MC1.root) file 
+    //  - calculate the hadroness for the method of RANDOM FOREST
+    //  - update input root file with the hadroness
+
+
+ if (JobB_RF_UP)
  {
     gLog << "=====================================================" << endl;
-    gLog << "Macro CT1Analysis : Start of Job B_NN_UP" << endl;
+    gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
 
     gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobB_NN_UP, RLookNN, WNN = " 
-         << JobB_NN_UP  << ",  " << RLookNN << ",  " << WNN << endl;
-
+    gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " 
+         << (JobB_RF_UP ? "kTRUE" : "kFALSE")  << ",  " 
+         << (RLookRF    ? "kTRUE" : "kFALSE")  << ",  "
+         << (RTree      ? "kTRUE" : "kFALSE")  << ",  "
+         << (WRF        ? "kTRUE" : "kFALSE")  << endl;
+
+
+    //--------------------------------------------
+    // parameters for the random forest
+    Int_t NumTrees = 100;
+    Int_t NumTry   =   3;
+    Int_t NdSize   =   1;
 
 
@@ -898,5 +919,4 @@
     outNameImage += typeInput;
     outNameImage += "2.root";
-
     //TString outNameImage = filenameData;
 
@@ -909,5 +929,5 @@
     filenameON += "ON";
     filenameON += "1.root";
-    Int_t howManyHadrons = 500000;
+    Int_t howManyHadrons = 35000;
     gLog << "filenameON = " << filenameON << ",   howManyHadrons = "
          << howManyHadrons  << endl; 
@@ -918,5 +938,5 @@
     filenameMC += "MC";
     filenameMC += "1.root";
-    Int_t howManyGammas = 50000;
+    Int_t howManyGammas = 35000;
     gLog << "filenameMC = " << filenameMC << ",   howManyGammas = "
          << howManyGammas   << endl; 
@@ -926,5 +946,5 @@
 
     TString outNameGammas = outPath;
-    outNameGammas += "matrix_gammas_";
+    outNameGammas += "RFmatrix_gammas_";
     outNameGammas += "MC";
     outNameGammas += ".root";
@@ -934,5 +954,5 @@
 
     TString outNameHadrons = outPath;
-    outNameHadrons += "matrix_hadrons_";
+    outNameHadrons += "RFmatrix_hadrons_";
     outNameHadrons += typeMatrixHadrons;
     outNameHadrons += ".root";
@@ -940,9 +960,19 @@
 
     MHMatrix matrixg("MatrixGammas");
+    matrixg.EnableGraphicalOutput();
+
     MHMatrix matrixh("MatrixHadrons");
+    matrixh.EnableGraphicalOutput();
+
+    //--------------------------------------------
+    // file of trees of the random forest 
+
+    TString outRF = outPath;
+    outRF += "RF.root";
+
 
    //*************************************************************************
    // read in matrices of look-alike events
-if (RLookNN)
+if (RLookRF)
   {
     const char* mtxName = "MatrixGammas";
@@ -961,5 +991,5 @@
 
     matrixg.Read(mtxName);
-    matrixg.Print("cols");
+    matrixg.Print("SizeCols");
 
 
@@ -981,5 +1011,5 @@
 
     matrixh.Read(mtxName);
-    matrixh.Print("cols");
+    matrixh.Print("SizeCols");
   }
 
@@ -987,5 +1017,5 @@
    //*************************************************************************
    // create matrices of look-alike events
-if (!RLookNN)
+if (!RLookRF)
   {
     gLog << "" << endl;
@@ -1011,13 +1041,23 @@
     MFParticleId fgamma("MMcEvt", '=', kGAMMA);
     fgamma.SetName("gammaID");
+
     MFParticleId fhadrons("MMcEvt", '!', kGAMMA);
     fhadrons.SetName("hadronID)");
 
-    MFEventSelector selectorg;
-    selectorg.SetNumSelectEvts(howManyGammas);
+    TString mgname("costhg");
+    MBinning bing("Binning"+mgname);
+    bing.SetEdges(10, 0., 1.0);
+
+    MH3 gref("cos(MMcEvt.fTelescopeTheta)");
+    gref.SetName(mgname);
+    MH::SetBinning(&gref.GetHist(), &bing);
+    for (Int_t i=1; i<=gref.GetNbins(); i++)
+      gref.GetHist().SetBinContent(i, 1.0);
+
+
+    MFEventSelector2 selectorg(gref);
+    selectorg.SetNumMax(howManyGammas);
     selectorg.SetName("selectGammas");
-    MFEventSelector selectorh;
-    selectorh.SetNumSelectEvts(howManyHadrons);
-    selectorh.SetName("selectHadrons");
+
 
     MFillH fillmatg("MatrixGammas");
@@ -1029,5 +1069,5 @@
     fillmath.SetName("fillHadrons");
 
-    
+
     //*****************************   fill gammas   ***  
     // entries in MFilterList
@@ -1065,8 +1105,26 @@
     tlistg.PrintStatistics(0, kTRUE);
 
-
+    
     //*****************************   fill hadrons   ***  
 
     gLog << " Hadrons :" << endl;
+
+    TString mhname("costhh");
+    MBinning binh("Binning"+mhname);
+    binh.SetEdges(10, 0., 1.0);
+
+    MH3 href("cos(MMcEvt.fTelescopeTheta)");
+    href.SetName(mhname);
+    MH::SetBinning(&href.GetHist(), &binh);
+    for (Int_t i=1; i<=href.GetNbins(); i++)
+      href.GetHist().SetBinContent(i, 1.0);
+
+    MFEventSelector2 selectorh(href);
+    //selectorh.SetNumMax(howManyHadrons);
+    // select as many hadrons as gammas
+    selectorh.SetNumMax(matrixg.GetM().GetNrows());
+    selectorh.SetName("selectHadrons");
+
+
     // entries in MFilterList
 
@@ -1102,80 +1160,9 @@
 
     tlisth.PrintStatistics(0, kTRUE);
+
+
+
     //*****************************************************  
 
-    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-    //
-    // ----------------------------------------------------------
-    //  Definition of the reference sample of look-alike events
-    //  (this is a subsample of the original sample)
-    // ----------------------------------------------------------
-    //
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    // Select a maximum of nmaxevts events from the sample of look-alike 
-    // events. They will form the reference sample.
-    Int_t nmaxevts  = 2000;
-
-    // target distribution for the variable in column refcolumn (start at 0);
-    //Int_t   refcolumn = 0;
-    //Int_t   nbins   =   5;
-    //Float_t frombin = 0.5;
-    //Float_t tobin   = 1.0;
-    //TH1F *thsh = new TH1F("thsh","target distribution", 
-    //                       nbins, frombin, tobin);
-    //thsh->SetDirectory(NULL);
-    //thsh->SetXTitle("cos( \\Theta)");
-    //thsh->SetYTitle("Counts");
-    //Float_t dbin = (tobin-frombin)/nbins;
-    //Float_t lbin = frombin +dbin*0.5;
-    //for (Int_t j=1; j<=nbins; j++) 
-    //{
-    //  thsh->Fill(lbin,1.0);
-    //  lbin += dbin;
-    //}
-
-    Int_t   refcolumn = 0;
-    MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning");
-    TH1F *thsh = new TH1F();
-    thsh->SetNameTitle("thsh","target distribution");
-    MH::SetBinning(thsh, binscostheta);
-    Int_t nbins = thsh->GetNbinsX();
-    Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900};
-    for (Int_t j=1; j<=nbins; j++) 
-    {
-      thsh->Fill(thsh->GetBinCenter(j), cont[j-1]);
-    }
-
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; 
-    gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
-         << refcolumn << ",  " << nmaxevts << endl;    
-
-    matrixg.EnableGraphicalOutput();
-    Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts);
-
-    if ( !matrixok ) 
-    {
-      gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
-      return;
-    }
-
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; 
-    gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
-         << refcolumn << ",  " << nmaxevts << endl;    
-
-    matrixh.EnableGraphicalOutput();
-    matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts);
-    delete thsh;
-
-    if ( !matrixok ) 
-    {
-      gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
-      return;
-    }
-    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
     // write out look-alike events
@@ -1189,5 +1176,5 @@
       // "gammas"
       gLog << "Gammas :" << endl;    
-      matrixg.Print("cols");
+      matrixg.Print("SizeCols");
 
       TFile writeg(outNameGammas, "RECREATE", "");
@@ -1201,5 +1188,5 @@
       // "hadrons"
       gLog << "Hadrons :" << endl;    
-      matrixh.Print("cols");
+      matrixh.Print("SizeCols");
 
       TFile writeh(outNameHadrons, "RECREATE", "");
@@ -1214,14 +1201,126 @@
 
 
+    MRanForest *fRanForest;
+    MRanTree *fRanTree;
     //-----------------------------------------------------------------
-    // Update the input files with the NN and SC hadronness
-    //
-
- if (WNN)
+    // read in the trees of the random forest 
+    if (RTree)
+    {
+      MParList plisttr;
+      MTaskList tlisttr;
+      plisttr.AddToList(&tlisttr);
+
+      MReadTree readtr("TREE", outRF);
+      readtr.DisableAutoScheme();
+
+      MRanForestFill rffill;
+      rffill.SetNumTrees(NumTrees);
+
+      // list of tasks for the loop over the trees
+
+      tlisttr.AddToList(&readtr);
+      tlisttr.AddToList(&rffill);
+
+      //-------------------
+      // Execute tree loop
+      //
+      MEvtLoop evtlooptr;
+      evtlooptr.SetParList(&plisttr);
+      if (!evtlooptr.Eventloop())
+        return;
+
+      tlisttr.PrintStatistics(0, kTRUE);
+
+
+    // get adresses of objects which are used in the next eventloop
+    fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
+    if (!fRanForest)
+    {
+        *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
+    if (!fRanTree)                                  
+    {                                                                          
+        *fLog << err << dbginf << "MRanTree not found... aborting." << endl;    
+        return kFALSE;
+    }
+
+    }
+
+    //-----------------------------------------------------------------
+    // grow the trees of the random forest (event loop = tree loop)
+
+    if (!RTree)
+    {
+
+    gLog << "" << endl;
+    gLog << "========================================================" << endl;
+    gLog << "Macro CT1Analysis : start growing trees" << endl;
+
+    MTaskList tlist2;
+    MParList plist2;
+    plist2.AddToList(&tlist2);
+
+    plist2.AddToList(&matrixg);
+    plist2.AddToList(&matrixh);
+
+    MRanForestGrow rfgrow2;
+    rfgrow2.SetNumTrees(NumTrees);
+    rfgrow2.SetNumTry(NumTry);
+    rfgrow2.SetNdSize(NdSize);
+
+    MWriteRootFile rfwrite2(outRF);
+    rfwrite2.AddContainer("MRanTree", "TREE");
+
+    // list of tasks for the loop over the trees
+    
+    tlist2.AddToList(&rfgrow2);
+    tlist2.AddToList(&rfwrite2);
+
+    //-------------------
+    // Execute tree loop
+    //
+    MEvtLoop treeloop;
+    treeloop.SetParList(&plist2);
+
+    if ( !treeloop.Eventloop() )
+        return;
+
+    tlist2.PrintStatistics(0, kTRUE);
+
+    // get adresses of objects which are used in the next eventloop
+    fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
+    if (!fRanForest)
+    {
+        *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
+    if (!fRanTree)                                  
+    {                                                                          
+        *fLog << err << dbginf << "MRanTree not found... aborting." << endl;    
+        return kFALSE;
+    }
+
+    }
+    // end of growing the trees of the random forest
+    //-----------------------------------------------------------------
+
+
+
+
+    //-----------------------------------------------------------------
+    // Update the input files with the RF hadronness
+    //
+
+ if (WRF)
   {
     gLog << "" << endl;
     gLog << "========================================================" << endl;
     gLog << "Update input file '" <<  filenameData 
-         << "' with the NN and SC hadronness" << endl;
+         << "' with the RF hadronness" << endl;
 
     MTaskList tliston;
@@ -1246,19 +1345,9 @@
 
     //.......................................................................
-    // calculation of hadroness for method of Nearest Neighbors
-
-    TString hadNNName = "HadNN";
-    MMultiDimDistCalc nncalc;
-    nncalc.SetUseNumRows(25);
-    nncalc.SetUseKernelMethod(kFALSE);
-    nncalc.SetHadronnessName(hadNNName);
-
-    //.......................................................................
-    // calculation of hadroness for the supercuts
-    // (=0.25 if fullfilled, =0.75 otherwise)
-
-    TString hadSCName = "HadSC";
-    MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc);
-    sccalc.SetHadronnessName(hadSCName);
+    // calculate hadronnes for method of RANDOM FOREST
+
+    TString hadRFName = "HadRF";
+    MRanForestCalc rfcalc;
+    rfcalc.SetHadronnessName(hadRFName);
 
     //.......................................................................
@@ -1278,6 +1367,5 @@
       write.AddContainer("MNewImagePar",  "Events");
 
-      write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
+      write.AddContainer(hadRFName,       "Events");
 
 
@@ -1287,4 +1375,5 @@
              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
 
+
     Float_t maxhadronness =  0.40;
     Float_t maxalpha      =  20.0;
@@ -1293,21 +1382,30 @@
     MFCT1SelFinal selfinalgh(fHilNameSrc);
     selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
-    selfinalgh.SetHadronnessName(hadNNName);
+    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 fillranfor("MHRanForest");
+    fillranfor.SetName("HRanForest");
+
+    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
+    fillhadrf.SetName("HhadRF");
 
     MFCT1SelFinal selfinal(fHilNameSrc);
     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
-    selfinal.SetHadronnessName(hadNNName);
+    selfinal.SetHadronnessName(hadRFName);
     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);
+    MFillH alpha(&alphaabs);
+    alpha.SetName("FillAlphaAbs");
 
     MFillH hfill1("MHHillas",    fHilName);
@@ -1332,7 +1430,9 @@
     InitBinnings(&pliston);
 
-    pliston.AddToList(&matrixg);
-    pliston.AddToList(&matrixh);
-
+    pliston.AddToList(fRanForest);
+    pliston.AddToList(fRanTree);
+
+    pliston.AddToList(&binsalphaabs);
+    pliston.AddToList(&alphaabs);
 
     //*****************************
@@ -1341,9 +1441,12 @@
     tliston.AddToList(&read);
 
-    tliston.AddToList(&nncalc);
-    tliston.AddToList(&sccalc);
-    tliston.AddToList(&fillhadnn);
-    tliston.AddToList(&fillhadsc);
-
+    tliston.AddToList(&rfcalc);
+    tliston.AddToList(&fillranfor);
+    tliston.AddToList(&fillhadrf);
+
+    tliston.AddToList(&write);
+    tliston.AddToList(&contfinalgh);
+
+    tliston.AddToList(&alpha);
     tliston.AddToList(&hfill1);
     tliston.AddToList(&hfill2);
@@ -1352,7 +1455,4 @@
     tliston.AddToList(&hfill5);
 
-    tliston.AddToList(&write);
-
-    tliston.AddToList(&contfinalgh);
     tliston.AddToList(&contfinal);
 
@@ -1373,9 +1473,10 @@
     tliston.PrintStatistics(0, kTRUE);
 
+
     //-------------------------------------------
     // Display the histograms
     //
-    pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
-    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
+    pliston.FindObject("MHRanForest")->DrawClone();
+    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
 
     pliston.FindObject("MHHillas")->DrawClone();
@@ -1385,13 +1486,49 @@
     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]");
+
+    Double_t alphasig = 13.1;
+    Double_t alphamin = 30.0;
+    Double_t alphamax = 90.0;
+    Int_t    degree   =    4;
+    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_NN_UP" << endl;
+    gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
     gLog << "=======================================================" << endl;
  }
   //---------------------------------------------------------------------
+
+
 
 
@@ -1413,13 +1550,36 @@
 
     gLog << "" << endl;
-    //gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = "
-    //     << JobB_SC_UP  << ",  " << WParSC << ",  " << WSC << endl;
-
+    gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = "
+         << (JobB_SC_UP ? "kTRUE" : "kFALSE")  << ",  " 
+         << (WParSC     ? "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_2209a";
+
+    gLog << "parSCfile = " << parSCfile << endl;
 
     //--------------------------------------------
     // file to be updated (either ON or MC)
 
-    TString typeInput = "ON";
-    //TString typeInput = "MC";
+    //TString typeInput = "ON";
+    TString typeInput = "MC";
     gLog << "typeInput = " << typeInput << endl;
 
@@ -1434,5 +1594,4 @@
     outNameImage += typeInput;
     outNameImage += "3.root";
-    outNameImage += "xxxxx";
     
 
@@ -1447,5 +1606,5 @@
     TString filenameTrain = outPath;
     filenameTrain += "ON";
-    filenameTrain += "2.root";
+    filenameTrain += "1.root";
     Int_t howManyTrain = 800000;
     gLog << "filenameTrain = " << filenameTrain << ",   howManyTrain = "
@@ -1455,6 +1614,6 @@
     TString filenameTest = outPath;
     filenameTest += "ON";
-    filenameTest += "2.root";
-    Int_t howManyTest = 100000;
+    filenameTest += "1.root";
+    Int_t howManyTest = 800000;
 
     gLog << "filenameTest = " << filenameTest << ",   howManyTest = "
@@ -1479,19 +1638,11 @@
 
     
-    //--------------------------------------------
-    // file which contains the optimum parameter values for the supercuts 
-
-    TString parSCfile = outPath;
-    parSCfile += "parSC";
-
-    gLog << "parSCfile = " << parSCfile << endl;
-
 
     //---------------------------------------------------------------------
-    // optimize supercuts using the training sample
-    // and test them using the test sample
-
-if (WParSC)
-  {
+    // Training and test matrices :
+    // - either create them and write them onto a file
+    // - or read them from a file
+
+
     MCT1FindSupercuts findsuper;
     findsuper.SetFilenameParam(parSCfile);
@@ -1500,5 +1651,5 @@
 
     //--------------------------
-    // create matrices 
+    // create matrices and write them onto files 
     if (!RMatrix)
     {
@@ -1509,6 +1660,6 @@
       {
         if ( !findsuper.DefineTrainTestMatrix(filenameTrain, 
-                                              howManyTrain, mh3,
-                                              howManyTest,  mh3) )
+                              howManyTrain, mh3, howManyTest,  mh3,
+                              fileMatrixTrain, fileMatrixTest)     )
         {
           *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl;
@@ -1519,5 +1670,6 @@
       else
       {
-        if ( !findsuper.DefineTrainMatrix(filenameTrain, howManyTrain, mh3) )
+        if ( !findsuper.DefineTrainMatrix(filenameTrain, 
+                              howManyTrain, mh3, fileMatrixTrain) )
         {
           *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl;
@@ -1525,5 +1677,6 @@
         }
 
-	if ( !findsuper.DefineTestMatrix( filenameTest,  howManyTest,  mh3) )
+	if ( !findsuper.DefineTestMatrix( filenameTest,  
+                              howManyTest,  mh3, fileMatrixTest)  )
         {
           *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl;
@@ -1531,24 +1684,124 @@
         }
       }
- 
-      // writing doesn't work ???
-      //findsuper.WriteMatrix(fileMatrixTrain, fileMatrixTest);
-    }
+     }
 
     //--------------------------
     // read matrices from files
-    //                              NOT YET TESTED
-    //if (RMatrix)
-    //  findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
+    //                              
+
+    if (RMatrix)
+      findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
     //--------------------------
 
 
-    //--------------------------
-    // optimize the supercuts
-    Bool_t rf = findsuper.FindParams();
+
+    //---------------------------------------------------------------------
+    // 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 (WParSC)
+  {
+    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) 
     {
-      *fLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
+       gLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
        return;
     }
@@ -1565,5 +1818,6 @@
     //
 
-
+ if (WSC)
+ {
     gLog << "" << endl;
     gLog << "========================================================" << endl;
@@ -1580,5 +1834,5 @@
     inparam.Close();
 
-    gLog << "Optimum parameter values for supercuts were read in from file '"
+    gLog << "Parameter values for supercuts were read in from file '"
          << parSCfile << "'" << endl;
 
@@ -1586,8 +1840,18 @@
     supercutsPar =  scin.GetParameters();
 
-    gLog << "Optimum parameter values for supercuts : " << endl;
-    for (Int_t i=0; i<72; i++)
+    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;
@@ -1604,5 +1868,5 @@
          << filenameData << "'" << endl;
     supercutsPar = supercuts.GetParameters();
-    for (Int_t i=0; i<72; i++)
+    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
     {
       gLog << supercutsPar[i] << ",  ";
@@ -1640,8 +1904,7 @@
 
 
- if (WSC)
-  {
       //MWriteRootFile write(outNameImage, "UPDATE");
-      MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
+      //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
+      MWriteRootFile write(outNameImage, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -1656,7 +1919,6 @@
       write.AddContainer("MNewImagePar",  "Events");
 
-      write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
-  }
+      write.AddContainer("HadRF",         "Events");
+      write.AddContainer(hadSCName,       "Events");
 
 
@@ -1693,4 +1955,7 @@
     MH3 alphaabs("abs(MHillasSrc.fAlpha)");
     alphaabs.SetName(mh3name);
+
+    TH1  &alphahist = alphaabs->GetHist();
+
     MFillH alpha(&alphaabs);
     alpha.SetName("FillAlphaAbs");
@@ -1716,7 +1981,9 @@
     pliston.AddToList(&tliston);
     InitBinnings(&pliston);
+
+    pliston.AddToList(&supercuts);
+
     pliston.AddToList(&binsalphaabs);
     pliston.AddToList(&alphaabs);
-    pliston.AddToList(&supercuts);
 
     //*****************************
@@ -1727,4 +1994,6 @@
     tliston.AddToList(&sccalc);
     tliston.AddToList(&fillhadsc);
+
+    tliston.AddToList(&write);
     tliston.AddToList(&contfinalgh);
 
@@ -1735,7 +2004,4 @@
     tliston.AddToList(&hfill4);
     tliston.AddToList(&hfill5);
-
-    if (WSC)
-      tliston.AddToList(&write);
 
     tliston.AddToList(&contfinal);
@@ -1776,9 +2042,10 @@
     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   =    4;
+    Int_t    degree   =    2;
     Double_t significance = -99.0;
     Bool_t   drawpoly  = kTRUE;
@@ -1788,5 +2055,5 @@
     MHFindSignificance findsig;
     findsig.SetRebin(kTRUE);
-    findsig.SetReduceDegree(kTRUE);
+    findsig.SetReduceDegree(kFALSE);
 
     findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
@@ -1805,4 +2072,6 @@
 
     DeleteBinnings(&pliston);
+ }
+
 
     gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl;
@@ -1810,661 +2079,4 @@
  }
   //---------------------------------------------------------------------
-
-
-  //---------------------------------------------------------------------
-  // Job B_RF_UP
-  //============
-
-
-    //  - create (or read in) the matrices of look-alike events for gammas 
-    //    and hadrons
-    //  - create (or read in) the trees
-    //  - then read ON1.root (or MC1.root) file 
-    //  - calculate the hadroness for the method of RANDOM FOREST
-    //  - update input root file with the hadroness
-
-
- if (JobB_RF_UP)
- {
-    gLog << "=====================================================" << endl;
-    gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
-
-    gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " 
-         << JobB_RF_UP  << ",  " << RLookRF << ",  " << RTree << ",  "
-         << WRF << endl;
-
-
-
-    //--------------------------------------------
-    // parameters for the random forest
-    Int_t NumTrees = 100;
-    Int_t NumTry   =   3;
-    Int_t NdSize   =   3;
-
-
-    //--------------------------------------------
-    // file to be updated (either ON or MC)
-
-    TString typeInput = "ON";
-    //TString typeInput = "MC";
-    gLog << "typeInput = " << typeInput << endl;
-
-    // name of input root file
-    TString filenameData = outPath;
-    filenameData += typeInput;
-    filenameData += "2.root";
-    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 generating the look-alike events
-    // "hadrons" :
-    TString filenameON = outPath;
-    filenameON += "ON";
-    filenameON += "1.root";
-    Int_t howManyHadrons = 1000000;
-    gLog << "filenameON = " << filenameON << ",   howManyHadrons = "
-         << howManyHadrons  << endl; 
-    
-
-    // "gammas" :
-    TString filenameMC = outPath;
-    filenameMC += "MC";
-    filenameMC += "1.root";
-    Int_t howManyGammas = 50000;
-    gLog << "filenameMC = " << filenameMC << ",   howManyGammas = "
-         << howManyGammas   << endl; 
-    
-    //--------------------------------------------
-    // files of look-alike events 
-
-    TString outNameGammas = outPath;
-    outNameGammas += "RFmatrix_gammas_";
-    outNameGammas += "MC";
-    outNameGammas += ".root";
-
-    TString typeMatrixHadrons = "ON";
-    gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
-
-    TString outNameHadrons = outPath;
-    outNameHadrons += "RFmatrix_hadrons_";
-    outNameHadrons += typeMatrixHadrons;
-    outNameHadrons += ".root";
-
-
-    MHMatrix matrixg("MatrixGammas");
-    MHMatrix matrixh("MatrixHadrons");
-
-    //--------------------------------------------
-    // file of trees of the random forest 
-
-    TString outRF = outPath;
-    outRF += "RF.root";
-
-
-   //*************************************************************************
-   // read in matrices of look-alike events
-if (RLookRF)
-  {
-    const char* mtxName = "MatrixGammas";
-
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << "Get matrix for (gammas)" << endl;
-    gLog << "matrix name        = " << mtxName << endl;
-    gLog << "name of root file  = " << outNameGammas << endl;
-    gLog << "" << endl;
-
-
-    // read in the object with the name 'mtxName' from file 'outNameGammas'
-    //
-    TFile fileg(outNameGammas); 
-
-    matrixg.Read(mtxName);
-    matrixg.Print("cols");
-
-
-    //***************************************************************** 
-
-    const char* mtxName = "MatrixHadrons";
-
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << " Get matrix for (hadrons)" << endl;
-    gLog << "matrix name        = " << mtxName << endl;
-    gLog << "name of root file  = " << outNameHadrons << endl;
-    gLog << "" << endl;
-
-
-    // read in the object with the name 'mtxName' from file 'outNameHadrons'
-    //
-    TFile fileh(outNameHadrons); 
-
-    matrixh.Read(mtxName);
-    matrixh.Print("cols");
-  }
-
-
-   //*************************************************************************
-   // create matrices of look-alike events
-if (!RLookRF)
-  {
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << " Create matrices of look-alike events" << endl;
-    gLog << " Gammas :" << endl;
-
-
-    MParList  plistg;
-    MTaskList tlistg;
-    MFilterList flistg;
-
-    MParList  plisth;
-    MTaskList tlisth;
-    MFilterList flisth;
-
-    MReadMarsFile  readg("Events", filenameMC);
-    readg.DisableAutoScheme();
-
-    MReadMarsFile  readh("Events", filenameON);
-    readh.DisableAutoScheme();
-
-    MFParticleId fgamma("MMcEvt", '=', kGAMMA);
-    fgamma.SetName("gammaID");
-    MFParticleId fhadrons("MMcEvt", '!', kGAMMA);
-    fhadrons.SetName("hadronID)");
-
-    MFEventSelector selectorg;
-    selectorg.SetNumSelectEvts(howManyGammas);
-    selectorg.SetName("selectGammas");
-    MFEventSelector selectorh;
-    selectorh.SetNumSelectEvts(howManyHadrons);
-    selectorh.SetName("selectHadrons");
-
-    MFillH fillmatg("MatrixGammas");
-    fillmatg.SetFilter(&flistg);
-    fillmatg.SetName("fillGammas");
-
-    MFillH fillmath("MatrixHadrons");
-    fillmath.SetFilter(&flisth);
-    fillmath.SetName("fillHadrons");
-
-    
-    //*****************************   fill gammas   ***  
-    // entries in MFilterList
-
-    flistg.AddToList(&fgamma);
-    flistg.AddToList(&selectorg);
-
-    //*****************************  
-    // entries in MParList
-    
-    plistg.AddToList(&tlistg);
-    InitBinnings(&plistg);
-
-    plistg.AddToList(&matrixg);
-
-    //*****************************
-    // entries in MTaskList
-    
-    tlistg.AddToList(&readg);
-    tlistg.AddToList(&flistg);
-    tlistg.AddToList(&fillmatg);
-
-    //*****************************
-
-    MProgressBar matrixbar;
-    MEvtLoop evtloopg;
-    evtloopg.SetParList(&plistg);
-    evtloopg.ReadEnv(env, "", printEnv);
-    evtloopg.SetProgressBar(&matrixbar);
-
-    Int_t maxevents = -1;
-    if (!evtloopg.Eventloop(maxevents))
-        return;
-
-    tlistg.PrintStatistics(0, kTRUE);
-
-
-    //*****************************   fill hadrons   ***  
-
-    gLog << " Hadrons :" << endl;
-    // entries in MFilterList
-
-    flisth.AddToList(&fhadrons);
-    flisth.AddToList(&selectorh);
-
-    //*****************************  
-    // entries in MParList
-    
-    plisth.AddToList(&tlisth);
-    InitBinnings(&plisth);
-
-    plisth.AddToList(&matrixh);
-
-    //*****************************
-    // entries in MTaskList
-    
-    tlisth.AddToList(&readh);
-    tlisth.AddToList(&flisth);
-    tlisth.AddToList(&fillmath);
-
-    //*****************************
-
-    MProgressBar matrixbar;
-    MEvtLoop evtlooph;
-    evtlooph.SetParList(&plisth);
-    evtlooph.ReadEnv(env, "", printEnv);
-    evtlooph.SetProgressBar(&matrixbar);
-
-    Int_t maxevents = -1;
-    if (!evtlooph.Eventloop(maxevents))
-        return;
-
-    tlisth.PrintStatistics(0, kTRUE);
-    //*****************************************************  
-
-    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-    //
-    // ----------------------------------------------------------
-    //  Definition of the reference sample of look-alike events
-    //  (this is a subsample of the original sample)
-    // ----------------------------------------------------------
-    //
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    // Select a maximum of nmaxevts events from the sample of look-alike 
-    // events. They will form the reference sample.
-    Int_t nmaxevts  = 10000;
-
-    // target distribution for the variable in column refcolumn (start at 0);
-    //Int_t   refcolumn = 0;
-    //Int_t   nbins   =   5;
-    //Float_t frombin = 0.5;
-    //Float_t tobin   = 1.0;
-    //TH1F *thsh = new TH1F("thsh","target distribution", 
-    //                       nbins, frombin, tobin);
-    //thsh->SetDirectory(NULL);
-    //thsh->SetXTitle("cos( \\Theta)");
-    //thsh->SetYTitle("Counts");
-    //Float_t dbin = (tobin-frombin)/nbins;
-    //Float_t lbin = frombin +dbin*0.5;
-    //for (Int_t j=1; j<=nbins; j++) 
-    //{
-    //  thsh->Fill(lbin,1.0);
-    //  lbin += dbin;
-    //}
-
-    Int_t   refcolumn = 0;
-    MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning");
-    TH1F *thsh = new TH1F();
-    thsh->SetNameTitle("thsh","target distribution");
-    MH::SetBinning(thsh, binscostheta);
-    Int_t nbins = thsh->GetNbinsX();
-    Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900};
-    for (Int_t j=1; j<=nbins; j++) 
-    {
-      thsh->Fill(thsh->GetBinCenter(j), cont[j-1]);
-    }
-
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; 
-    gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
-         << refcolumn << ",  " << nmaxevts << endl;    
-
-    matrixg.EnableGraphicalOutput();
-    Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts);
-
-    if ( !matrixok ) 
-    {
-      gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
-      return;
-    }
-
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; 
-    gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
-         << refcolumn << ",  " << nmaxevts << endl;    
-
-    matrixh.EnableGraphicalOutput();
-    matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts);
-    delete thsh;
-
-    if ( !matrixok ) 
-    {
-      gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
-      return;
-    }
-    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-    // write out look-alike events
-
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << "Write out look=alike events" << endl;
-
-
-      //-------------------------------------------
-      // "gammas"
-      gLog << "Gammas :" << endl;    
-      matrixg.Print("cols");
-
-      TFile writeg(outNameGammas, "RECREATE", "");
-      matrixg.Write();
-
-      gLog << "" << endl;
-      gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
-           << outNameGammas << endl;
-
-      //-------------------------------------------
-      // "hadrons"
-      gLog << "Hadrons :" << endl;    
-      matrixh.Print("cols");
-
-      TFile writeh(outNameHadrons, "RECREATE", "");
-      matrixh.Write();
-
-      gLog << "" << endl;
-      gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
-           << outNameHadrons << endl;
-
-  }
-   //**********   end of creating matrices of look-alike events   ***********
-
-
-    MRanForest *fRanForest;
-    MRanTree *fRanTree;
-    //-----------------------------------------------------------------
-    // read in the trees of the random forest 
-    if (RTree)
-    {
-      MParList plisttr;
-      MTaskList tlisttr;
-      plisttr.AddToList(&tlisttr);
-
-      MReadTree readtr("TREE", outRF);
-      readtr.DisableAutoScheme();
-
-      MRanForestFill rffill;
-      rffill.SetNumTrees(NumTrees);
-
-      // list of tasks for the loop over the trees
-
-      tlisttr.AddToList(&readtr);
-      tlisttr.AddToList(&rffill);
-
-      //-------------------
-      // Execute tree loop
-      //
-      MEvtLoop evtlooptr;
-      evtlooptr.SetParList(&plisttr);
-      if (!evtlooptr.Eventloop())
-        return;
-
-      tlisttr.PrintStatistics(0, kTRUE);
-
-
-    // get adresses of objects which are used in the next eventloop
-    fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
-    if (!fRanForest)
-    {
-        *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
-    if (!fRanTree)                                  
-    {                                                                          
-        *fLog << err << dbginf << "MRanTree not found... aborting." << endl;    
-        return kFALSE;
-    }
-
-    }
-
-    //-----------------------------------------------------------------
-    // grow the trees of the random forest (event loop = tree loop)
-
-    if (!RTree)
-    {
-
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << "Macro CT1Analysis : start growing trees" << endl;
-
-    MTaskList tlist2;
-    MParList plist2;
-    plist2.AddToList(&tlist2);
-
-    plist2.AddToList(&matrixg);
-    plist2.AddToList(&matrixh);
-
-    MRanForestGrow rfgrow2;
-    rfgrow2.SetNumTrees(NumTrees);
-    rfgrow2.SetNumTry(NumTry);
-    rfgrow2.SetNdSize(NdSize);
-
-    MWriteRootFile rfwrite2(outRF);
-    rfwrite2.AddContainer("MRanTree", "TREE");
-
-    // list of tasks for the loop over the trees
-    
-    tlist2.AddToList(&rfgrow2);
-    tlist2.AddToList(&rfwrite2);
-
-    //-------------------
-    // Execute tree loop
-    //
-    MEvtLoop treeloop;
-    treeloop.SetParList(&plist2);
-
-    if ( !treeloop.Eventloop() )
-        return;
-
-    tlist2.PrintStatistics(0, kTRUE);
-
-    // get adresses of objects which are used in the next eventloop
-    fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
-    if (!fRanForest)
-    {
-        *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
-    if (!fRanTree)                                  
-    {                                                                          
-        *fLog << err << dbginf << "MRanTree not found... aborting." << endl;    
-        return kFALSE;
-    }
-
-    }
-    // end of growing the trees of the random forest
-    //-----------------------------------------------------------------
-
-
-
-
-    //-----------------------------------------------------------------
-    // Update the input files with the RF hadronness
-    //
-
- if (WRF)
-  {
-    gLog << "" << endl;
-    gLog << "========================================================" << endl;
-    gLog << "Update input file '" <<  filenameData 
-         << "' with the RF hadronness" << endl;
-
-    MTaskList tliston;
-    MParList pliston;
-
-
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameData);
-    read.DisableAutoScheme();
-
-    TString fHilName    = "MHillas"; 
-    TString fHilNameExt = "MHillasExt"; 
-    TString fHilNameSrc = "MHillasSrc"; 
-    TString fImgParName = "MNewImagePar"; 
-
-    //.......................................................................
-    // calculate hadronnes for method of RANDOM FOREST
-
-    TString hadRFName = "HadRF";
-    MRanForestCalc rfcalc;
-    rfcalc.SetHadronnessName(hadRFName);
-
-    //.......................................................................
-
-      //MWriteRootFile write(outNameImage, "UPDATE");
-      MWriteRootFile write(outNameImage, "RECREATE");
-
-      write.AddContainer("MRawRunHeader", "RunHeaders");
-      write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      write.AddContainer("ThetaOrig",     "Events");
-      write.AddContainer("MSrcPosCam",    "Events");
-      write.AddContainer("MSigmabar",     "Events");
-      write.AddContainer("MHillas",       "Events");
-      write.AddContainer("MHillasExt",    "Events");
-      write.AddContainer("MHillasSrc",    "Events");
-      write.AddContainer("MNewImagePar",  "Events");
-
-      write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
-
-      write.AddContainer(hadRFName,       "Events");
-
-
-    //-----------------------------------------------------------------
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-
-    Float_t maxhadronness =  0.40;
-    Float_t maxalpha      =  20.0;
-    Float_t maxdist       =  10.0;
-
-    MFCT1SelFinal selfinalgh(fHilNameSrc);
-    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
-    selfinalgh.SetHadronnessName(hadRFName);
-    selfinalgh.SetName("SelFinalgh");
-    MContinue contfinalgh(&selfinalgh);
-    contfinalgh.SetName("ContSelFinalgh");
-
-    MFillH fillranfor("MHRanForest");
-    fillranfor.SetName("HRanForest");
-
-    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
-    fillhadrf.SetName("HhadRF");
-
-    MFCT1SelFinal selfinal(fHilNameSrc);
-    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
-    selfinal.SetHadronnessName(hadRFName);
-    selfinal.SetName("SelFinal");
-    MContinue contfinal(&selfinal);
-    contfinal.SetName("ContSelFinal");
-
-
-    MFillH hfill1("MHHillas",    fHilName);
-    hfill1.SetName("HHillas");
-
-    MFillH hfill2("MHStarMap",   fHilName);
-    hfill2.SetName("HStarMap");
-
-    MFillH hfill3("MHHillasExt",    fHilNameSrc);
-    hfill3.SetName("HHillasExt");
-    
-    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
-    hfill4.SetName("HHillasSrc");    
-
-    MFillH hfill5("MHNewImagePar", fImgParName);
-    hfill5.SetName("HNewImagePar");
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-    pliston.AddToList(fRanForest);
-    pliston.AddToList(fRanTree);
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-
-    tliston.AddToList(&rfcalc);
-    tliston.AddToList(&fillranfor);
-    tliston.AddToList(&fillhadrf);
-
-    tliston.AddToList(&hfill1);
-    tliston.AddToList(&hfill2);
-    tliston.AddToList(&hfill3);
-    tliston.AddToList(&hfill4);
-    tliston.AddToList(&hfill5);
-
-    tliston.AddToList(&write);
-
-    tliston.AddToList(&contfinalgh);
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-
-    Int_t maxevents = -1;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-    pliston.FindObject("MHRanForest")->DrawClone();
-    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
-
-    pliston.FindObject("MHHillas")->DrawClone();
-    pliston.FindObject("MHHillasExt")->DrawClone();
-    pliston.FindObject("MHHillasSrc")->DrawClone();
-    pliston.FindObject("MHNewImagePar")->DrawClone();
-    pliston.FindObject("MHStarMap")->DrawClone();
-
-    DeleteBinnings(&pliston);
-  }
-
-    gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
-    gLog << "=======================================================" << endl;
- }
-  //---------------------------------------------------------------------
-
 
 
@@ -2484,5 +2096,6 @@
 
     gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobC = " << JobC  << endl;
+    gLog << "Macro CT1Analysis : JobC = " 
+         << (JobC       ? "kTRUE" : "kFALSE")  << endl;
 
 
@@ -2522,5 +2135,5 @@
     // names of hadronness containers
 
-    TString hadNNName = "HadNN";
+    //TString hadNNName = "HadNN";
     TString hadSCName = "HadSC";
     TString hadRFName = "HadRF";
@@ -2540,11 +2153,11 @@
     MFCT1SelFinal selfinalgh(fHilNameSrc);
     selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
-    selfinalgh.SetHadronnessName(hadNNName);
+    selfinalgh.SetHadronnessName(hadSCName);
     selfinalgh.SetName("SelFinalgh");
     MContinue contfinalgh(&selfinalgh);
     contfinalgh.SetName("ContSelFinalgh");
 
-    MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
-    fillhadnn.SetName("HhadNN");
+    //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
+    //fillhadnn.SetName("HhadNN");
     MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
     fillhadsc.SetName("HhadSC");
@@ -2554,5 +2167,5 @@
     MFCT1SelFinal selfinal(fHilNameSrc);
     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
-    selfinal.SetHadronnessName(hadNNName);
+    selfinal.SetHadronnessName(hadSCName);
     selfinal.SetName("SelFinal");
     MContinue contfinal(&selfinal);
@@ -2588,5 +2201,5 @@
     tliston.AddToList(&read);
 
-    tliston.AddToList(&fillhadnn);
+    //tliston.AddToList(&fillhadnn);
     tliston.AddToList(&fillhadsc);
     tliston.AddToList(&fillhadrf);
@@ -2623,5 +2236,5 @@
     //
 
-    pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
+    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
@@ -2657,5 +2270,6 @@
     gLog << "" << endl;
     gLog << "Macro CT1Analysis : JobD = " 
-         << JobD  << endl;
+         << (JobD        ? "kTRUE" : "kFALSE")  << endl;
+
 
     // type of data to be analysed
@@ -2868,5 +2482,7 @@
     gLog << "" << endl;
     gLog << "Macro CT1Analysis : JobE_XX, WXX = " 
-         << JobE_XX  << ",  " << WXX << endl;
+         << (JobE_XX ? "kTRUE" : "kFALSE")  << ",  " 
+         << (WXX     ? "kTRUE" : "kFALSE")  << endl;
+
 
     // type of data to be analysed
@@ -3060,4 +2676,7 @@
     gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
 
+
+  if (WXX)
+  {
     //-----------------------------------------------------------
     //
@@ -3119,7 +2738,7 @@
     //.......................................................................
 
-    if (WXX)
-    {
-      MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+
+      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+      MWriteRootFile write(filenameDataout, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -3139,5 +2758,5 @@
 
       write.AddContainer("MEnergyEst",    "Events");
-    }
+
 
     //-----------------------------------------------------------------
@@ -3265,4 +2884,6 @@
 
     DeleteBinnings(&pliston);
+
+  }
 
     gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
@@ -3292,5 +2913,7 @@
     gLog << "" << endl;
     gLog << "Macro CT1Analysis : JobF_XX, WFX = " 
-         << JobF_XX  << ",  " << WFX << endl;
+         << (JobF_XX ? "kTRUE" : "kFALSE")  << ",  " 
+         << (WFX     ? "kTRUE" : "kFALSE")  << endl;
+
 
     // type of data to be analysed
@@ -3493,4 +3116,7 @@
     gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
 
+
+  if (WFX)
+  {
     //-----------------------------------------------------------
     //
@@ -3552,10 +3178,9 @@
     //.......................................................................
 
-    //if (WFX)
-    //{
       gLog << "CT1Analysis.C : write root file '" << filenameDataout 
            << "'" << endl;
    
       //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+      /*
       MWriteRootFile write(filenameDataout, "RECREATE");
 
@@ -3571,10 +3196,10 @@
       write.AddContainer("MNewImagePar",  "Events");
 
-      write.AddContainer("HadNN",         "Events");
+      //write.AddContainer("HadNN",         "Events");
       write.AddContainer("HadSC",         "Events");
       write.AddContainer("HadRF",         "Events");
 
       write.AddContainer("MEnergyEst",    "Events");
-      //}
+      */
 
     //-----------------------------------------------------------------
@@ -3590,6 +3215,6 @@
     contfinalgh.SetName("ContSelFinalgh");
 
-    MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
-    fillhadnn.SetName("HhadNN");
+    //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
+    //fillhadnn.SetName("HhadNN");
     MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
     fillhadsc.SetName("HhadSC");
@@ -3600,15 +3225,15 @@
     // calculate estimated energy
 
-    //MEnergyEstParam eeston(fHilName);
-    //eeston.Add(fHilNameSrc);
-
-    //eeston.SetCoeffA(parA);
-    //eeston.SetCoeffB(parB);
+    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);
+    //MEnergyEstParamDanielMkn421 eeston(fHilName);
+    //eeston.Add(fHilNameSrc);
     //eeston.SetCoeffA(parA);
     //eeston.SetCoeffB(parB);
@@ -3684,8 +3309,7 @@
     tliston.AddToList(&eeston);
 
-    if (WFX)
-      tliston.AddToList(&write);
-
-    tliston.AddToList(&fillhadnn);
+    //tliston.AddToList(&write);
+
+    //tliston.AddToList(&fillhadnn);
     tliston.AddToList(&fillhadsc);
     tliston.AddToList(&fillhadrf);
@@ -3727,19 +3351,38 @@
     //
 
-    pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
+    //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 << "before Robert's code" << endl;
+
 
 //rwagner write all relevant histograms onto a file
 
+  if (WRobert)
+  {
     gLog << "=======================================================" << endl;
     gLog << "Write results onto file '" << filenameResults << "'" << endl;
@@ -3751,5 +3394,7 @@
         TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
     alphaHist->Write();
-    
+    gLog << "Alpha plot has been written out" << endl;    
+
+
     MHAlphaEnergyTheta* aetH = 
       (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
@@ -3757,4 +3402,5 @@
     aetHist->SetName("aetHist");
     aetHist->Write();
+    gLog << "AlphaEnergyTheta plot has been written out" << endl;    
 
     MHAlphaEnergyTime* aetH2 = 
@@ -3764,4 +3410,5 @@
 //     aetHist2->DrawClone();
     aetHist2->Write();
+    gLog << "AlphaEnergyTime plot has been written out" << endl;    
 
     MHThetabarTime* tbt = 
@@ -3770,4 +3417,5 @@
     tbtHist->SetName("tbtHist");
     tbtHist->Write();
+    gLog << "ThetabarTime plot has been written out" << endl;    
 
     MHEnergyTime* ent = 
@@ -3776,4 +3424,5 @@
     entHist->SetName("entHist");
     entHist->Write();
+    gLog << "EnergyTime plot has been written out" << endl;    
     
     MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
@@ -3782,4 +3431,6 @@
     timeHist->SetTitle("Time diffs");
     timeHist->Write();
+    gLog << "TimeDiffTheta plot has been written out" << endl;    
+
 
     MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
@@ -3788,4 +3439,5 @@
     timeHist2->SetTitle("Time diffs");
     timeHist2->Write();
+    gLog << "TimeDiffTime plot has been written out" << endl;    
 
 //rwagner write also collareas to same file
@@ -3793,12 +3445,18 @@
     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 F_XX" << endl;
Index: trunk/MagicSoft/Mars/macros/CT1EgyEst.C
===================================================================
--- trunk/MagicSoft/Mars/macros/CT1EgyEst.C	(revision 2317)
+++ trunk/MagicSoft/Mars/macros/CT1EgyEst.C	(revision 2357)
@@ -33,5 +33,5 @@
 #include "MHMatrix.h"
 #include "MEnergyEstParam.h"
-#include "MEnergyEstParamDanielMkn421.h"
+//#include "MEnergyEstParamDanielMkn421.h"
 #include "MMatrixLoop.h"
 #include "MChisqEval.h"
@@ -206,13 +206,13 @@
   //
 
-  //MEnergyEstParam eest2(hilName);
+  MEnergyEstParam eest2(hilName);
+  eest2.Add(hilSrcName);
+
+  eest2.SetCoeffA(parA);
+  eest2.SetCoeffB(parB);
+
+  // estimate energy using Daniel's parameters
+  //MEnergyEstParamDanielMkn421 eest2(hilName);
   //eest2.Add(hilSrcName);
-
-  //eest2.SetCoeffA(parA);
-  //eest2.SetCoeffB(parB);
-
-  // estimate energy using Daniel's parameters
-  MEnergyEstParamDanielMkn421 eest2(hilName);
-  eest2.Add(hilSrcName);
 
 
