Changeset 2357 for trunk/MagicSoft/Mars

Timestamp:

09/24/03 13:29:00 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• macros/CT1Analysis.C (modified) (98 diffs)
• macros/CT1EgyEst.C (modified) (2 diffs)
• manalysis/MCT1FindSupercuts.cc (modified) (29 diffs)
• manalysis/MCT1FindSupercuts.h (modified) (1 diff)
• manalysis/MCT1Supercuts.cc (modified) (8 diffs)
• manalysis/MCT1Supercuts.h (modified) (5 diffs)
• manalysis/MCT1SupercutsCalc.cc (modified) (10 diffs)
• manalysis/MCT1SupercutsCalc.h (modified) (4 diffs)
• manalysis/MMinuitInterface.cc (modified) (5 diffs)
• mfilter/MFEventSelector2.cc (modified) (8 diffs)
• mfilter/MFEventSelector2.h (modified) (1 diff)
• mhist/MHCT1Supercuts.cc (modified) (1 diff)
• mhist/MHFindSignificance.cc (modified) (10 diffs)
• mhist/MHMatrix.h (modified) (1 diff)
• mmontecarlo/MMcEnergyEst.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+
+
+ 2003/09/24: Wolfgang Wittek
+
+   * mfilter/MFEventSelector2.[h,cc]
+     - execution statistics added
+
+   * mhist/MHFindSignificance.cc
+     - add fHist->UseCurrentStyle()
+       to get the y-axis + labels drawn
+
+   * mhist/MHMatrix.h
+     - replace   Int_t fNumRow  //!   
+            by   Int_t fNumRow  //
+       because otherwise fNumRow is not defined when MHMatrix object is read in
+       after it had been written out
+
+   * mhist/MHCT1Supercuts.cc
+     - change title of object
+
+   * manalysis/MMinuitInterface.cc
+     - add arguments maxcalls and tolerance to SIMPLEX call
+
+   * manalysis/MCT1SupercutsCalc.[h,cc]
+     - add variables asymmetry, conc, leakage
+
+   * manalysis/MCT1Supercuts.[h,cc]
+     - add variables asymmetry, conc, leakage
+     - add TArrayD fStepsizes (initial step sizes for the parameters)
+
+   * manalysis/MCT1FindSupercuts.cc
+     - replace MGeomCamCT1Daniel by MGeomCamCT1
+     - arguments 'parSCinit', 'params' and 'steps' added in FindParams() ;
+         parSCinit is the name of the file containing the initial
+         values of the parameters
+         'params' and 'steps' are the initial values in case parSCinit == ""
+     - add member functions GetMatrixTrain() and GetMatrixTest()
+     - remove member function WriteMatrix()
+       because it didn't work; now the matrices are written out in
+       DefineTrainMatrix(), DefineTestMatri() and DefineTrainTestMatrix()
+
+   * macros/CT1EgyEst.C
+     - don't use Daniel's energy estimator
+
+   * mmontecarlo/MMcEnergyEst.cc
+     - extend printout of comments
+
+
 
  2003/09/17: Abelardo Moralejo
@@ -51 +99 @@
      - added 'Time Spectra of Cosmics' button
      - added size argument to instatiation of MHFadcCam
-
 
 

trunk/MagicSoft/Mars/macros/CT1Analysis.C

@@ -22 +22 @@
 #include "MWriteRootFile.h"
 
+//#include "TH3D.h"
 
 
@@ -77 +78 @@
 
         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 +95 @@
           zedge[8-i] = cos(yedge[i]/kRad2Deg);
         }
-        TArrayD zed(9,zedge);
+        TArrayD zed;
+        zed.Set(9,zedge);
         bincosth->SetEdges(zed);
         plist->AddToList(bincosth);
@@ -116 +118 @@
 }
 
+
 void DeleteBinnings(MParList *plist)
 {
@@ -122 +125 @@
         TObject *bin;
 
+        //--------------------------------------------
+        bin = plist->FindObject("BinningE");
+        if (bin) delete bin;
+
+        //--------------------------------------------
+
         bin = plist->FindObject("BinningSize");
         if (bin) delete bin;
@@ -159 +168 @@
         if (bin) delete bin;
 
-        //robert
+
+        // robert ----------------------------------------------
         bin = plist->FindObject("BinningAlphaFlux");
         if (bin) delete bin;
@@ -170 +180 @@
 }
 
+
 //************************************************************************
 void CT1Analysis()
 {
+
+  gLog << "Entry CT1Analysis()" << endl;
+
       gLog.SetNoColors();
 
@@ -189 +203 @@
 
       // 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 +231 @@
     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 +243 @@
     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 +304 @@
     //  - 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 +335 @@
     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 +434 @@
     contstandard.SetName("SelStandard");
 
-    if (WON1)
-    {
-      MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
+      //MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
+      MWriteRootFile write(outNameImage, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -438 +447 @@
       write.AddContainer("MHillasSrc",    "Events");
       write.AddContainer("MNewImagePar",  "Events");
-    }
+
 
     //$$$$$$$$$$$$$$$$$$$$$$$$$$$$
@@ -500 +509 @@
 
     Int_t maxevents = -1;
-    //Int_t maxevents = 1000;
     if ( !evtloop.Eventloop(maxevents) )
         return;
@@ -587 +595 @@
     gLog << "" << endl;
     gLog << "Macro CT1Analysis : JobA_MC, WMC1 = " 
-         << JobA_MC  << ",  " << WMC1 << endl;
+         << (JobA_MC ? "kTRUE" : "kFALSE")  << ",  " 
+         << (WMC1    ? "kTRUE" : "kFALSE")  << endl;
 
 
@@ -688 +697 @@
     MCT1ReadPreProc read(filenamein);
 
+    MBlindPixelCalc blindbeforepad;
+    blindbeforepad.SetUseBlindPixels();
+    blindbeforepad.SetName("BlindBeforePadding");
+
     MBlindPixelCalc blind;
     blind.SetUseBlindPixels();
+    blind.SetName("BlindAfterPadding");
 
     MFCT1SelBasic selbasic;
@@ -758 +772 @@
 
 
-    if (WMC1)
-    {
-      MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
+      //MWriteRootFile &write = *(new MWriteRootFile(outNameImage));
+      MWriteRootFile write(outNameImage, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -772 +785 @@
       write.AddContainer("MHillasSrc",    "Events");
       write.AddContainer("MNewImagePar",  "Events");
-    }
+
 
 
@@ -788 +801 @@
 
     tlist.AddToList(&read);
+    tlist.AddToList(&blindbeforepad);
     tlist.AddToList(&padthomas);
     tlist.AddToList(&blind);
@@ -855 +869 @@
  }
 
-
-
   //---------------------------------------------------------------------
-  // 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 +919 @@
     outNameImage += typeInput;
     outNameImage += "2.root";
-
     //TString outNameImage = filenameData;
 
@@ -909 +929 @@
     filenameON += "ON";
     filenameON += "1.root";
-    Int_t howManyHadrons = 500000;
+    Int_t howManyHadrons = 35000;
     gLog << "filenameON = " << filenameON << ",   howManyHadrons = "
          << howManyHadrons  << endl; 
@@ -918 +938 @@
     filenameMC += "MC";
     filenameMC += "1.root";
-    Int_t howManyGammas = 50000;
+    Int_t howManyGammas = 35000;
     gLog << "filenameMC = " << filenameMC << ",   howManyGammas = "
          << howManyGammas   << endl; 
@@ -926 +946 @@
 
     TString outNameGammas = outPath;
-    outNameGammas += "matrix_gammas_";
+    outNameGammas += "RFmatrix_gammas_";
     outNameGammas += "MC";
     outNameGammas += ".root";
@@ -934 +954 @@
 
     TString outNameHadrons = outPath;
-    outNameHadrons += "matrix_hadrons_";
+    outNameHadrons += "RFmatrix_hadrons_";
     outNameHadrons += typeMatrixHadrons;
     outNameHadrons += ".root";
@@ -940 +960 @@
 
     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 +991 @@
 
     matrixg.Read(mtxName);
-    matrixg.Print("cols");
+    matrixg.Print("SizeCols");
 
 
@@ -981 +1011 @@
 
     matrixh.Read(mtxName);
-    matrixh.Print("cols");
+    matrixh.Print("SizeCols");
   }
 
@@ -987 +1017 @@
    //*************************************************************************
    // create matrices of look-alike events
-if (!RLookNN)
+if (!RLookRF)
   {
     gLog << "" << endl;
@@ -1011 +1041 @@
     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 +1069 @@
     fillmath.SetName("fillHadrons");
 
-    
+
     //*****************************   fill gammas   ***  
     // entries in MFilterList
@@ -1065 +1105 @@
     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 +1160 @@
 
     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 +1176 @@
       // "gammas"
       gLog << "Gammas :" << endl;    
-      matrixg.Print("cols");
+      matrixg.Print("SizeCols");
 
       TFile writeg(outNameGammas, "RECREATE", "");
@@ -1201 +1188 @@
       // "hadrons"
       gLog << "Hadrons :" << endl;    
-      matrixh.Print("cols");
+      matrixh.Print("SizeCols");
 
       TFile writeh(outNameHadrons, "RECREATE", "");
@@ -1214 +1201 @@
 
 
+    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 +1345 @@
 
     //.......................................................................
-    // 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 +1367 @@
       write.AddContainer("MNewImagePar",  "Events");
 
-      write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
+      write.AddContainer(hadRFName,       "Events");
 
 
@@ -1287 +1375 @@
              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
 
+
     Float_t maxhadronness =  0.40;
     Float_t maxalpha      =  20.0;
@@ -1293 +1382 @@
     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 +1430 @@
     InitBinnings(&pliston);
 
-    pliston.AddToList(&matrixg);
-    pliston.AddToList(&matrixh);
-
+    pliston.AddToList(fRanForest);
+    pliston.AddToList(fRanTree);
+
+    pliston.AddToList(&binsalphaabs);
+    pliston.AddToList(&alphaabs);
 
     //*****************************
@@ -1341 +1441 @@
     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 +1455 @@
     tliston.AddToList(&hfill5);
 
-    tliston.AddToList(&write);
-
-    tliston.AddToList(&contfinalgh);
     tliston.AddToList(&contfinal);
 
@@ -1373 +1473 @@
     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 +1486 @@
     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 +1550 @@
 
     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 +1594 @@
     outNameImage += typeInput;
     outNameImage += "3.root";
-    outNameImage += "xxxxx";
     
 
@@ -1447 +1606 @@
     TString filenameTrain = outPath;
     filenameTrain += "ON";
-    filenameTrain += "2.root";
+    filenameTrain += "1.root";
     Int_t howManyTrain = 800000;
     gLog << "filenameTrain = " << filenameTrain << ",   howManyTrain = "
@@ -1455 +1614 @@
     TString filenameTest = outPath;
     filenameTest += "ON";
-    filenameTest += "2.root";
-    Int_t howManyTest = 100000;
+    filenameTest += "1.root";
+    Int_t howManyTest = 800000;
 
     gLog << "filenameTest = " << filenameTest << ",   howManyTest = "
@@ -1479 +1638 @@
 
     
-    //--------------------------------------------
-    // 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 +1651 @@
 
     //--------------------------
-    // create matrices 
+    // create matrices and write them onto files 
     if (!RMatrix)
     {
@@ -1509 +1660 @@
       {
         if ( !findsuper.DefineTrainTestMatrix(filenameTrain, 
-                                              howManyTrain, mh3,
-                                              howManyTest,  mh3) )
+                              howManyTrain, mh3, howManyTest,  mh3,
+                              fileMatrixTrain, fileMatrixTest)     )
         {
           *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl;
@@ -1519 +1670 @@
       else
       {
-        if ( !findsuper.DefineTrainMatrix(filenameTrain, howManyTrain, mh3) )
+        if ( !findsuper.DefineTrainMatrix(filenameTrain, 
+                              howManyTrain, mh3, fileMatrixTrain) )
         {
           *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl;
@@ -1525 +1677 @@
         }
 
-        if ( !findsuper.DefineTestMatrix( filenameTest,  howManyTest,  mh3) )
+        if ( !findsuper.DefineTestMatrix( filenameTest,  
+                              howManyTest,  mh3, fileMatrixTest)  )
         {
           *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl;
@@ -1531 +1684 @@
         }
       }
- 
-      // 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 +1818 @@
     //
 
-
+ if (WSC)
+ {
     gLog << "" << endl;
     gLog << "========================================================" << endl;
@@ -1580 +1834 @@
     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 +1840 @@
     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 +1868 @@
          << 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 +1904 @@
 
 
- 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 +1919 @@
       write.AddContainer("MNewImagePar",  "Events");
 
-      write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
-  }
+      write.AddContainer("HadRF",         "Events");
+      write.AddContainer(hadSCName,       "Events");
 
 
@@ -1693 +1955 @@
     MH3 alphaabs("abs(MHillasSrc.fAlpha)");
     alphaabs.SetName(mh3name);
+
+    TH1  &alphahist = alphaabs->GetHist();
+
     MFillH alpha(&alphaabs);
     alpha.SetName("FillAlphaAbs");
@@ -1716 +1981 @@
     pliston.AddToList(&tliston);
     InitBinnings(&pliston);
+
+    pliston.AddToList(&supercuts);
+
     pliston.AddToList(&binsalphaabs);
     pliston.AddToList(&alphaabs);
-    pliston.AddToList(&supercuts);
 
     //*****************************
@@ -1727 +1994 @@
     tliston.AddToList(&sccalc);
     tliston.AddToList(&fillhadsc);
+
+    tliston.AddToList(&write);
     tliston.AddToList(&contfinalgh);
 
@@ -1735 +2004 @@
     tliston.AddToList(&hfill4);
     tliston.AddToList(&hfill5);
-
-    if (WSC)
-      tliston.AddToList(&write);
 
     tliston.AddToList(&contfinal);
@@ -1776 +2042 @@
     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 +2055 @@
     MHFindSignificance findsig;
     findsig.SetRebin(kTRUE);
-    findsig.SetReduceDegree(kTRUE);
+    findsig.SetReduceDegree(kFALSE);
 
     findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
@@ -1805 +2072 @@
 
     DeleteBinnings(&pliston);
+ }
+
 
     gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl;
@@ -1810 +2079 @@
  }
   //---------------------------------------------------------------------
-
-
-  //---------------------------------------------------------------------
-  // 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 +2096 @@
 
     gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobC = " << JobC  << endl;
+    gLog << "Macro CT1Analysis : JobC = " 
+         << (JobC       ? "kTRUE" : "kFALSE")  << endl;
 
 
@@ -2522 +2135 @@
     // names of hadronness containers
 
-    TString hadNNName = "HadNN";
+    //TString hadNNName = "HadNN";
     TString hadSCName = "HadSC";
     TString hadRFName = "HadRF";
@@ -2540 +2153 @@
     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 +2167 @@
     MFCT1SelFinal selfinal(fHilNameSrc);
     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
-    selfinal.SetHadronnessName(hadNNName);
+    selfinal.SetHadronnessName(hadSCName);
     selfinal.SetName("SelFinal");
     MContinue contfinal(&selfinal);
@@ -2588 +2201 @@
     tliston.AddToList(&read);
 
-    tliston.AddToList(&fillhadnn);
+    //tliston.AddToList(&fillhadnn);
     tliston.AddToList(&fillhadsc);
     tliston.AddToList(&fillhadrf);
@@ -2623 +2236 @@
     //
 
-    pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
+    //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
@@ -2657 +2270 @@
     gLog << "" << endl;
     gLog << "Macro CT1Analysis : JobD = " 
-         << JobD  << endl;
+         << (JobD        ? "kTRUE" : "kFALSE")  << endl;
+
 
     // type of data to be analysed
@@ -2868 +2482 @@
     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 +2676 @@
     gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
 
+
+  if (WXX)
+  {
     //-----------------------------------------------------------
     //
@@ -3119 +2738 @@
     //.......................................................................
 
-    if (WXX)
-    {
-      MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+
+      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+      MWriteRootFile write(filenameDataout, "RECREATE");
 
       write.AddContainer("MRawRunHeader", "RunHeaders");
@@ -3139 +2758 @@
 
       write.AddContainer("MEnergyEst",    "Events");
-    }
+
 
     //-----------------------------------------------------------------
@@ -3265 +2884 @@
 
     DeleteBinnings(&pliston);
+
+  }
 
     gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
@@ -3292 +2913 @@
     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 +3116 @@
     gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
 
+
+  if (WFX)
+  {
     //-----------------------------------------------------------
     //
@@ -3552 +3178 @@
     //.......................................................................
 
-    //if (WFX)
-    //{
       gLog << "CT1Analysis.C : write root file '" << filenameDataout 
            << "'" << endl;
    
       //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+      /*
       MWriteRootFile write(filenameDataout, "RECREATE");
 
@@ -3571 +3196 @@
       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 +3215 @@
     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 +3225 @@
     // 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 +3309 @@
     tliston.AddToList(&eeston);
 
-    if (WFX)
-      tliston.AddToList(&write);
-
-    tliston.AddToList(&fillhadnn);
+    //tliston.AddToList(&write);
+
+    //tliston.AddToList(&fillhadnn);
     tliston.AddToList(&fillhadsc);
     tliston.AddToList(&fillhadrf);
@@ -3727 +3351 @@
     //
 
-    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 +3394 @@
         TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
     alphaHist->Write();
-    
+    gLog << "Alpha plot has been written out" << endl;    
+
+
     MHAlphaEnergyTheta* aetH = 
       (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
@@ -3757 +3402 @@
     aetHist->SetName("aetHist");
     aetHist->Write();
+    gLog << "AlphaEnergyTheta plot has been written out" << endl;    
 
     MHAlphaEnergyTime* aetH2 = 
@@ -3764 +3410 @@
 //     aetHist2->DrawClone();
     aetHist2->Write();
+    gLog << "AlphaEnergyTime plot has been written out" << endl;    
 
     MHThetabarTime* tbt = 
@@ -3770 +3417 @@
     tbtHist->SetName("tbtHist");
     tbtHist->Write();
+    gLog << "ThetabarTime plot has been written out" << endl;    
 
     MHEnergyTime* ent = 
@@ -3776 +3424 @@
     entHist->SetName("entHist");
     entHist->Write();
+    gLog << "EnergyTime plot has been written out" << endl;    
     
     MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
@@ -3782 +3431 @@
     timeHist->SetTitle("Time diffs");
     timeHist->Write();
+    gLog << "TimeDiffTheta plot has been written out" << endl;    
+
 
     MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
@@ -3788 +3439 @@
     timeHist2->SetTitle("Time diffs");
     timeHist2->Write();
+    gLog << "TimeDiffTime plot has been written out" << endl;    
 
 //rwagner write also collareas to same file
@@ -3793 +3445 @@
     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;

trunk/MagicSoft/Mars/macros/CT1EgyEst.C

@@ -33 +33 @@
 #include "MHMatrix.h"
 #include "MEnergyEstParam.h"
-#include "MEnergyEstParamDanielMkn421.h"
+//#include "MEnergyEstParamDanielMkn421.h"
 #include "MMatrixLoop.h"
 #include "MChisqEval.h"
@@ -206 +206 @@
   //
 
-  //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);
 
 

trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc

@@ -57 +57 @@
 #include "MFEventSelector2.h"
 #include "MFillH.h"
-#include "MGeomCamCT1Daniel.h"
+//#include "MGeomCamCT1Daniel.h"
+#include "MGeomCamCT1.h"
 #include "MH3.h"
 #include "MHCT1Supercuts.h"
@@ -95 +96 @@
                          Double_t *par, Int_t iflag)
 {
+    //cout <<  "entry fcnSupercuts" << endl;
+
     //-------------------------------------------------------------
     // save pointer to the MINUIT object for optimizing the supercuts
@@ -106 +109 @@
 
     MParList *plistfcn   = (MParList*)evtloopfcn->GetParList();
+    //MTaskList *tasklistfcn   = (MTaskList*)plistfcn->FindObject("MTaskList");
 
     MCT1Supercuts *super = (MCT1Supercuts*)plistfcn->FindObject("MCT1Supercuts");
@@ -130 +134 @@
     // for testing
     //TArrayD checkparameters = super->GetParameters();
+    //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ",  " 
+    //     << fNparx  << endl;
     //gLog << "fcnsupercuts : i, par, checkparameters =" << endl;
     //for (Int_t i=0; i<fNparx; i++)
@@ -142 +148 @@
     evtloopfcn->Eventloop();
 
+    //tasklistfcn->PrintStatistics(0, kTRUE);
+
     MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3");
     if (!alpha)
@@ -147 +155 @@
 
     TH1 &alphaHist = alpha->GetHist();
-    //alphaHist.SetXTitle("|\\alpha|  [\\circ]");
+    alphaHist.SetName("alpha-fcnSupercuts");
 
     //-------------------------------------------
@@ -161 +169 @@
     const Double_t alphamin = 30.0;
     const Double_t alphamax = 90.0;
-    const Int_t    degree   =    4;
+    const Int_t    degree   =    2;
 
     Bool_t drawpoly;
@@ -240 +248 @@
     //---------------------------
     // camera geometry is needed for conversion mm ==> degree
-    fCam = new MGeomCamCT1Daniel; 
-
-    // matrices to contain the the training/test samples
+    //fCam = new MGeomCamCT1Daniel; 
+    fCam = new MGeomCamCT1; 
+
+    // matrices to contain the training/test samples
     fMatrixTrain = new MHMatrix("MatrixTrain");
     fMatrixTest  = new MHMatrix("MatrixTest");
@@ -277 +286 @@
 //
 Bool_t MCT1FindSupercuts::DefineTrainMatrix(const TString &nametrain,
-                          const Int_t howmanytrain, MH3 &hreftrain)
+                          const Int_t howmanytrain, MH3 &hreftrain,
+                          const TString &filetrain)
 {
     if (nametrain.IsNull() || howmanytrain <= 0)
@@ -343 +353 @@
     *fLog << "=============================================" << endl;
 
+    //--------------------------
+    // write out training matrix
+
+    if (filetrain != "")
+    {
+      TFile filetr(filetrain, "RECREATE");
+      fMatrixTrain->Write();
+      filetr.Close();
+
+      *fLog << "MCT1FindSupercuts::DefineTrainMatrix; Training matrix was written onto file '"
+            << filetrain << "'" << endl;
+    }
+
+
     return kTRUE;
 }
@@ -356 +380 @@
 //
 Bool_t MCT1FindSupercuts::DefineTestMatrix(const TString &nametest,
-                          const Int_t howmanytest, MH3 &hreftest)
+                          const Int_t howmanytest, MH3 &hreftest,
+                          const TString &filetest)
 {
     if (nametest.IsNull() || howmanytest<=0)
@@ -422 +447 @@
     *fLog << "=============================================" << endl;
 
+    //--------------------------
+    // write out test matrix
+
+    if (filetest != "")
+    {
+      TFile filete(filetest, "RECREATE", "");
+      fMatrixTest->Write();
+      filete.Close();
+
+      *fLog << "MCT1FindSupercuts::DefineTestMatrix; Test matrix was written onto file '"
+            << filetest << "'" << endl;
+    }
+
+
   return kTRUE;
 }
@@ -434 +473 @@
                           const TString &name,
                           const Int_t howmanytrain, MH3 &hreftrain,
-                          const Int_t howmanytest,  MH3 &hreftest)
+                          const Int_t howmanytest,  MH3 &hreftest,
+                          const TString &filetrain, const TString &filetest)
 {
     *fLog << "=============================================" << endl;
@@ -511 +551 @@
     evtloop.SetProgressBar(&bar);
 
-    if (!evtloop.Eventloop())
+    Int_t maxev = -1;
+    if (!evtloop.Eventloop(maxev))
       return kFALSE;
 
@@ -520 +561 @@
     if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
     {
-      *fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events ("
+      *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
             << generatedtrain 
             << ") is incompatible with the no.of requested events ("
@@ -530 +571 @@
     if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
     {
-      *fLog << "MCT1FindSupercuts::DefineTestMatrix; no.of generated events ("
+      *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
             << generatedtest 
             << ") is incompatible with the no.of requested events ("
@@ -541 +582 @@
 
 
+    //--------------------------
+    // write out training matrix
+
+    if (filetrain != "")
+    {
+      TFile filetr(filetrain, "RECREATE");
+      fMatrixTrain->Write();
+      filetr.Close();
+
+      *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '"
+            << filetrain << "'" << endl;
+    }
+
+    //--------------------------
+    // write out test matrix
+
+    if (filetest != "")
+    {
+      TFile filete(filetest, "RECREATE", "");
+      fMatrixTest->Write();
+      filete.Close();
+
+      *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '"
+            << filetest << "'" << endl;
+    }
+
   return kTRUE;
 }
@@ -548 +615 @@
 // Read training and test matrices from files
 //
-//
-// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$
 //
 
@@ -563 +628 @@
   *fLog << "MCT1FindSupercuts::ReadMatrix; Training matrix was read in from file '"
         << filetrain << "'" << endl;
+  filetr.Close();
 
 
@@ -574 +640 @@
   *fLog << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '"
         << filetest << "'" << endl;
-
+  filete.Close();
 
   return kTRUE;  
 }
 
-// --------------------------------------------------------------------------
-//
-// Write training and test matrices onto files
-//
-//
-// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$
-//
-//
-Bool_t MCT1FindSupercuts::WriteMatrix(const TString &filetrain, const TString &filetest)
-{
-  //--------------------------
-  // write out training matrix
-
-  TFile filetr(filetrain, "RECREATE");
-
-  *fLog << "nach TFile" << endl;
-
-  fMatrixTrain->Write();
-
-  *fLog << "MCT1FindSupercuts::WriteMatrix; Training matrix was written onto file '"
-        << filetrain << "'" << endl;
-
-
-  //--------------------------
-  // write out test matrix
-
-  TFile filete(filetest, "RECREATE");
-  fMatrixTest->Print("SizeCols");
-  fMatrixTest->Write();
-
-  *fLog << "MCT1FindSupercuts::WriteMatrix; Test matrix was written onto file '"
-        << filetest << "'" << endl;
-
-
-  return kTRUE;  
-}
 
 //------------------------------------------------------------------------
@@ -645 +675 @@
 //                       put the supercuts hadronness
 //
-// - for the minimization, the starting values of the parameters are taken as 
-//   MCT1Supercuts::GetParameters(fVinit)
+// - for the minimization, the starting values of the parameters are taken  
+//     - from file parSCinit (if it is != "")
+//     - or from the arrays params and/or steps 
 //
 //----------------------------------------------------------------------
-Bool_t MCT1FindSupercuts::FindParams()
+Bool_t MCT1FindSupercuts::FindParams(TString parSCinit,
+                                     TArrayD &params, TArrayD &steps)
 {
     // Setup the event loop which will be executed in the 
     //                 fcnSupercuts function  of MINUIT
     //
+    // parSCinit is the name of the file containing the initial values 
+    // of the parameters; 
+    // if parSCinit = ""   'params' and 'steps' are taken as initial values
 
     *fLog << "=============================================" << endl;
@@ -685 +720 @@
     MMatrixLoop loop(fMatrixTrain);
 
+    //--------------------------------
     // create container for the supercut parameters
-    // and set them to their default values
+    // and set them to their initial values
     MCT1Supercuts super;
 
+    // take initial values from file parSCinit
+    if (parSCinit != "")
+    {
+      TFile inparam(parSCinit);
+      super.Read("MCT1Supercuts");
+      inparam.Close();
+      *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from file "
+            << parSCinit << endl;
+    }
+
+    // take initial values from 'params' and/or 'steps'
+    else if (params.GetSize() != 0  || steps.GetSize()  != 0 )
+    {
+      if (params.GetSize()  != 0)
+      {
+        *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from 'params'"
+              << endl;
+        super.SetParameters(params);
+      }
+      if (steps.GetSize()  != 0)
+      {
+        *fLog << "MCT1FindSupercuts::FindParams; initial step sizes are taken from 'steps'"
+              << endl;
+        super.SetStepsizes(steps);
+      }
+    }
+    else
+    {
+        *fLog << "MCT1FindSupercuts::FindParams; initial values and step sizes are taken from the MCT1Supercits constructor"
+              << endl;
+    }
+
+
+    //--------------------------------
     // calculate supercuts hadronness
     fCalcHadTrain->SetHadronnessName(fHadronnessName);
@@ -761 +831 @@
     //
 
-    // get initial values of parameters from MCT1Supercuts
+    // get initial values of parameters 
     fVinit = super.GetParameters();
+    fStep  = super.GetStepsizes();
 
     TString name[fVinit.GetSize()];
@@ -776 +847 @@
         name[i]   = "p";
         name[i]  += i+1;
-        fStep[i]  = TMath::Abs(fVinit[i]/10.0);
+        //fStep[i]  = TMath::Abs(fVinit[i]/10.0);
         fLimlo[i] = -100.0;
         fLimup[i] =  100.0;
         fFix[i]   =     0;
-
-        // vary only first 48 parameters
-        if (i >= 48)
-        {
-          fStep[i] = 0.0;
-          fFix[i]  =   1;
-        }
-    }
+    }
+
+    // these parameters make no sense, fix them at 0.0
+    fVinit[33] = 0.0;
+    fStep[33]  = 0.0;
+    fFix[33]   = 1;
+
+    fVinit[36] = 0.0;
+    fStep[36]  = 0.0;
+    fFix[36]   = 1;
+
+    fVinit[39] = 0.0;
+    fStep[39]  = 0.0;
+    fFix[39]   = 1;
+
+    fVinit[41] = 0.0;
+    fStep[41]  = 0.0;
+    fFix[41]   = 1;
+
+    fVinit[44] = 0.0;
+    fStep[44]  = 0.0;
+    fFix[44]   = 1;
+
+    fVinit[47] = 0.0;
+    fStep[47]  = 0.0;
+    fFix[47]   = 1;
+
+    // vary only first 48 parameters
+    //for (UInt_t i=0; i<fNpar; i++)
+    //{
+    //    if (i >= 48)
+    //  {
+    //      fStep[i] = 0.0;
+    //      fFix[i]  =   1;
+    //  }
+    //}
+ 
 
     // -------------------------------------------
@@ -910 +1010 @@
     MH3 alphabefore("MatrixTest[7]");
     alphabefore.SetName(mh3NameB);
+
+    TH1 &alphahistb = alphabefore.GetHist();
+    alphahistb.SetName("alphaBefore-TestParams");
+
     MFillH fillalphabefore(&alphabefore);
     fillalphabefore.SetName("FillAlphaBefore");
@@ -924 +1028 @@
     MH3 alphaafter("MatrixTest[7]");
     alphaafter.SetName(mh3NameA);
+
+    TH1 &alphahista = alphaafter.GetHist();
+    alphahista.SetName("alphaAfter-TestParams");
+
     MFillH fillalphaafter(&alphaafter);
     fillalphaafter.SetName("FillAlphaAfter");
@@ -978 +1086 @@
     TH1  &alphaHist2 = alphaAfter->GetHist();
     alphaHist2.SetXTitle("|\\alpha|  [\\circ]");
+    alphaHist2.SetName("alpha-TestParams)");
 
     TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300);
@@ -999 +1108 @@
     const Double_t alphamin = 30.0;
     const Double_t alphamax = 90.0;
-    const Int_t    degree   =    4;
+    const Int_t    degree   =    2;
     const Bool_t   drawpoly  = kTRUE;
     const Bool_t   fitgauss  = kFALSE;
@@ -1025 +1134 @@
     return kTRUE;
 }
+

trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h

@@ -89 +89 @@
 
   Bool_t DefineTrainMatrix(const TString &name, const Int_t howmany,
-                           MH3 &href);
+                           MH3 &href, const TString &filetrain); 
 
   Bool_t DefineTestMatrix(const TString &name, const Int_t howmany,
-                          MH3 &href);
+                          MH3 &href,  const TString &filetest);
 
   Bool_t DefineTrainTestMatrix(const TString &name, 
-                               const Int_t howmanytrain, MH3 &hreftrain,
-                               const Int_t howmanytest,  MH3 &hreftest);
+                         const Int_t howmanytrain, MH3 &hreftrain,
+                         const Int_t howmanytest,  MH3 &hreftest,
+                         const TString &filetrain, const TString &filetest);
+
+  MHMatrix *GetMatrixTrain() { return fMatrixTrain; }
+  MHMatrix *GetMatrixTest()  { return fMatrixTest;  }
 
   Bool_t ReadMatrix( const TString &filetrain, const TString &filetest);
-  Bool_t WriteMatrix(const TString &filetrain, const TString &filetest);
 
-  Bool_t FindParams();
+  Bool_t FindParams(TString parSCinit, TArrayD &params, TArrayD &steps);
   Bool_t TestParams();
 

trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc

@@ -50 +50 @@
 //
 MCT1Supercuts::MCT1Supercuts(const char *name, const char *title)
-    : fParameters(72),
-    fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+8),
+  : fParameters(104), fStepsizes(104),
+    fLengthUp(fParameters.GetArray()),   fLengthLo(fParameters.GetArray()+8),
     fWidthUp(fParameters.GetArray()+16), fWidthLo(fParameters.GetArray()+24),
-    fDistUp(fParameters.GetArray()+32), fDistLo(fParameters.GetArray()+40),
-    fAsymUp(fParameters.GetArray()+48), fAsymLo(fParameters.GetArray()+56),
-    fAlphaUp(fParameters.GetArray()+64)
+    fDistUp(fParameters.GetArray()+32),  fDistLo(fParameters.GetArray()+40),
+    fAsymUp(fParameters.GetArray()+48),  fAsymLo(fParameters.GetArray()+56),
+    fConcUp(fParameters.GetArray()+64),  fConcLo(fParameters.GetArray()+72),
+    fLeakage1Up(fParameters.GetArray()+80), fLeakage1Lo(fParameters.GetArray()+88),
+    fAlphaUp(fParameters.GetArray()+96)
 {
     fName  = name  ? name  : "MCT1Supercuts";
@@ -71 +73 @@
 void MCT1Supercuts::InitParameters()
 {
+    //---------------------------------------------------
+    //  these are Daniel's original values for Mkn 421
+
     fLengthUp[0] =  0.315585;
     fLengthUp[1] =  0.001455;
@@ -80 +85 @@
     fLengthUp[7] = -0.013463;
 
+    fLengthLo[0] =  0.151530;
+    fLengthLo[1] =  0.028323;
+    fLengthLo[2] =  0.510707;
+    fLengthLo[3] =  0.053089;
+    fLengthLo[4] =  0.013708;
+    fLengthLo[5] =  2.357993;
+    fLengthLo[6] =  0.000080;
+    fLengthLo[7] = -0.007157;
+
     fWidthUp[0] =  0.145412;
     fWidthUp[1] = -0.001771;
@@ -89 +103 @@
     fWidthUp[7] = -0.016703;
 
+    fWidthLo[0] =  0.089187;
+    fWidthLo[1] = -0.006430;
+    fWidthLo[2] =  0.074442;
+    fWidthLo[3] =  0.003738;
+    fWidthLo[4] = -0.004256;
+    fWidthLo[5] = -0.014101;
+    fWidthLo[6] =  0.006126;
+    fWidthLo[7] = -0.002849;
+
     fDistUp[0] =  1.787943;
     fDistUp[1] =  0;
@@ -98 +121 @@
     fDistUp[7] =  0;
 
-    fLengthLo[0] =  0.151530;
-    fLengthLo[1] =  0.028323;
-    fLengthLo[2] =  0.510707;
-    fLengthLo[3] =  0.053089;
-    fLengthLo[4] =  0.013708;
-    fLengthLo[5] =  2.357993;
-    fLengthLo[6] =  0.000080;
-    fLengthLo[7] = -0.007157;
-
-    fWidthLo[0] =  0.089187;
-    fWidthLo[1] = -0.006430;
-    fWidthLo[2] =  0.074442;
-    fWidthLo[3] =  0.003738;
-    fWidthLo[4] = -0.004256;
-    fWidthLo[5] = -0.014101;
-    fWidthLo[6] =  0.006126;
-    fWidthLo[7] = -0.002849;
-
     fDistLo[0] =  0.589406;
     fDistLo[1] =  0;
@@ -124 +129 @@
     fDistLo[6] = -0.001750;
     fDistLo[7] =  0;
-
-    fAsymUp[0] =  0.061267;
-    fAsymUp[1] =  0.014462;
-    fAsymUp[2] =  0.014327;
-    fAsymUp[3] =  0.014540;
-    fAsymUp[4] =  0.013391;
-    fAsymUp[5] =  0.012319;
-    fAsymUp[6] =  0.010444;
-    fAsymUp[7] =  0.008328;
-
-    fAsymLo[0] = -0.012055;
-    fAsymLo[1] =  0.009157;
-    fAsymLo[2] =  0.005441;
-    fAsymLo[3] =  0.000399;
-    fAsymLo[4] =  0.001433;
-    fAsymLo[5] = -0.002050;
-    fAsymLo[6] = -0.000104;
-    fAsymLo[7] = -0.001188;
+    
+
+    // dummy values
+
+    fAsymUp[0] =  1.e10;
+    fAsymUp[1] =  0.0;
+    fAsymUp[2] =  0.0;
+    fAsymUp[3] =  0.0;
+    fAsymUp[4] =  0.0;
+    fAsymUp[5] =  0.0;
+    fAsymUp[6] =  0.0;
+    fAsymUp[7] =  0.0;
+
+    fAsymLo[0] = -1.e10;
+    fAsymLo[1] =  0.0;
+    fAsymLo[2] =  0.0;
+    fAsymLo[3] =  0.0;
+    fAsymLo[4] =  0.0;
+    fAsymLo[5] =  0.0;
+    fAsymLo[6] =  0.0;
+    fAsymLo[7] =  0.0;
+
+    fConcUp[0] =  1.e10;
+    fConcUp[1] =  0.0;
+    fConcUp[2] =  0.0;
+    fConcUp[3] =  0.0;
+    fConcUp[4] =  0.0;
+    fConcUp[5] =  0.0;
+    fConcUp[6] =  0.0;
+    fConcUp[7] =  0.0;
+
+    fConcLo[0] = -1.e10;
+    fConcLo[1] =  0.0;
+    fConcLo[2] =  0.0;
+    fConcLo[3] =  0.0;
+    fConcLo[4] =  0.0;
+    fConcLo[5] =  0.0;
+    fConcLo[6] =  0.0;
+    fConcLo[7] =  0.0;
+
+    fLeakage1Up[0] =  1.e10;
+    fLeakage1Up[1] =  0.0;
+    fLeakage1Up[2] =  0.0;
+    fLeakage1Up[3] =  0.0;
+    fLeakage1Up[4] =  0.0;
+    fLeakage1Up[5] =  0.0;
+    fLeakage1Up[6] =  0.0;
+    fLeakage1Up[7] =  0.0;
+
+    fLeakage1Lo[0] = -1.e10;
+    fLeakage1Lo[1] =  0.0;
+    fLeakage1Lo[2] =  0.0;
+    fLeakage1Lo[3] =  0.0;
+    fLeakage1Lo[4] =  0.0;
+    fLeakage1Lo[5] =  0.0;
+    fLeakage1Lo[6] =  0.0;
+    fLeakage1Lo[7] =  0.0;
 
     fAlphaUp[0] = 13.123440; 
@@ -151 +195 @@
     fAlphaUp[6] = 0;
     fAlphaUp[7] = 0;
+
+    //---------------------------------------------------
+    // fStepsizes 
+    // if == 0.0    the parameter will be fixed in the minimization
+    //    != 0.0    initial step sizes for the parameters
+
+    // LengthUp
+    fStepsizes[0] = 0.03;
+    fStepsizes[1] = 0.0002;
+    fStepsizes[2] = 0.02;
+    fStepsizes[3] = 0.0006;
+    fStepsizes[4] = 0.0002;
+    fStepsizes[5] = 0.002;
+    fStepsizes[6] = 0.0008;
+    fStepsizes[7] = 0.002;
+
+    // LengthLo
+    fStepsizes[8]  = 0.02;
+    fStepsizes[9]  = 0.003;
+    fStepsizes[10] = 0.05;
+    fStepsizes[11] = 0.006;
+    fStepsizes[12] = 0.002;
+    fStepsizes[13] = 0.3;
+    fStepsizes[14] = 0.0001;
+    fStepsizes[15] = 0.0008;
+
+    // WidthUp
+    fStepsizes[16] = 0.02;
+    fStepsizes[17] = 0.0002;
+    fStepsizes[18] = 0.006;
+    fStepsizes[19] = 0.003;
+    fStepsizes[20] = 0.002;
+    fStepsizes[21] = 0.006;
+    fStepsizes[22] = 0.002;
+    fStepsizes[23] = 0.002;
+
+    // WidthLo
+    fStepsizes[24] = 0.009;
+    fStepsizes[25] = 0.0007;
+    fStepsizes[26] = 0.008;
+    fStepsizes[27] = 0.0004;
+    fStepsizes[28] = 0.0005;
+    fStepsizes[29] = 0.002;
+    fStepsizes[30] = 0.0007;
+    fStepsizes[31] = 0.003;
+
+    // DistUp
+    fStepsizes[32] = 0.2;
+    fStepsizes[33] = 0.0;
+    fStepsizes[34] = 0.3;
+    fStepsizes[35] = 0.02;
+    fStepsizes[36] = 0.0;
+    fStepsizes[37] = 0.03;
+    fStepsizes[38] = 0.02;
+    fStepsizes[39] = 0.0;
+
+    // DistLo
+    fStepsizes[40] = 0.06;
+    fStepsizes[41] = 0.0;
+    fStepsizes[42] = 0.009;
+    fStepsizes[43] = 0.0008;
+    fStepsizes[44] = 0.0;
+    fStepsizes[45] = 0.005;
+    fStepsizes[46] = 0.0002;
+    fStepsizes[47] = 0.0;
+
+    // AsymUp
+    fStepsizes[48] = 0.0;
+    fStepsizes[49] = 0.0;
+    fStepsizes[50] = 0.0;
+    fStepsizes[51] = 0.0;
+    fStepsizes[52] = 0.0;
+    fStepsizes[53] = 0.0;
+    fStepsizes[54] = 0.0;
+    fStepsizes[55] = 0.0;
+
+    // AsymLo
+    fStepsizes[56] = 0.0;
+    fStepsizes[57] = 0.0;
+    fStepsizes[58] = 0.0;
+    fStepsizes[59] = 0.0;
+    fStepsizes[60] = 0.0;
+    fStepsizes[61] = 0.0;
+    fStepsizes[62] = 0.0;
+    fStepsizes[63] = 0.0;
+
+    // ConcUp
+    fStepsizes[64] = 0.0;
+    fStepsizes[65] = 0.0;
+    fStepsizes[66] = 0.0;
+    fStepsizes[67] = 0.0;
+    fStepsizes[68] = 0.0;
+    fStepsizes[69] = 0.0;
+    fStepsizes[70] = 0.0;
+    fStepsizes[71] = 0.0;
+
+    // ConcLo
+    fStepsizes[72] = 0.0;
+    fStepsizes[73] = 0.0;
+    fStepsizes[74] = 0.0;
+    fStepsizes[75] = 0.0;
+    fStepsizes[76] = 0.0;
+    fStepsizes[77] = 0.0;
+    fStepsizes[78] = 0.0;
+    fStepsizes[79] = 0.0;
+
+    // Leakage1Up
+    fStepsizes[80] = 0.0;
+    fStepsizes[81] = 0.0;
+    fStepsizes[82] = 0.0;
+    fStepsizes[83] = 0.0;
+    fStepsizes[84] = 0.0;
+    fStepsizes[85] = 0.0;
+    fStepsizes[86] = 0.0;
+    fStepsizes[87] = 0.0;
+
+    // Leakage1Lo
+    fStepsizes[88] = 0.0;
+    fStepsizes[89] = 0.0;
+    fStepsizes[90] = 0.0;
+    fStepsizes[91] = 0.0;
+    fStepsizes[92] = 0.0;
+    fStepsizes[93] = 0.0;
+    fStepsizes[94] = 0.0;
+    fStepsizes[95] = 0.0;
+
+    // AlphaUp
+    fStepsizes[96]  = 0.0;
+    fStepsizes[97]  = 0.0;
+    fStepsizes[98]  = 0.0;
+    fStepsizes[99]  = 0.0;
+    fStepsizes[100] = 0.0;
+    fStepsizes[101] = 0.0;
+    fStepsizes[102] = 0.0;
+    fStepsizes[103] = 0.0;
 }
 
@@ -173 +352 @@
 }
 
-
-
-
-
-
-
-
-
-
-
-
+// --------------------------------------------------------------------------
+//
+// Set the step sizes from the array 'd'
+//
+//
+Bool_t MCT1Supercuts::SetStepsizes(const TArrayD &d)
+{
+    if (d.GetSize() != fStepsizes.GetSize())
+    {
+        *fLog << err << "Sizes of d and of fStepsizes are different : "
+              << d.GetSize() << ",  " << fStepsizes.GetSize() << endl;
+        return kFALSE;
+    }
+
+    fStepsizes = d;
+
+    return kTRUE;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h

@@ -14 +14 @@
 private:
     TArrayD fParameters; // supercut parameters
+    TArrayD fStepsizes;  // step sizes of supercut parameters
 
     Double_t *fLengthUp; //!
@@ -23 +24 @@
     Double_t *fAsymUp;   //!
     Double_t *fAsymLo;   //!
+
+    Double_t *fConcUp;   //!
+    Double_t *fConcLo;   //!
+    Double_t *fLeakage1Up;   //!
+    Double_t *fLeakage1Lo;   //!
+
     Double_t *fAlphaUp;  //!
+
 
 public:
@@ -31 +39 @@
 
     Bool_t SetParameters(const TArrayD &d);
+    Bool_t SetStepsizes(const TArrayD &d);
+
     const TArrayD &GetParameters() const { return fParameters; }
+    const TArrayD &GetStepsizes()  const { return fStepsizes;  }
 
     const Double_t *GetLengthUp() const { return fLengthUp; }
@@ -41 +52 @@
     const Double_t *GetAsymUp() const   { return fAsymUp; }
     const Double_t *GetAsymLo() const   { return fAsymLo; }
+
+    const Double_t *GetConcUp() const   { return fConcUp; }
+    const Double_t *GetConcLo() const   { return fConcLo; }
+
+    const Double_t *GetLeakage1Up() const   { return fLeakage1Up; }
+    const Double_t *GetLeakage1Lo() const   { return fLeakage1Lo; }
+
     const Double_t *GetAlphaUp() const  { return fAlphaUp; }
 
@@ -47 +65 @@
 
 #endif
+
+
+

trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc

@@ -44 +44 @@
 #include "MHillasExt.h"
 #include "MHillasSrc.h"
+#include "MNewImagePar.h"
 #include "MMcEvt.hxx"
 #include "MCerPhotEvt.h"
@@ -65 +66 @@
 
 MCT1SupercutsCalc::MCT1SupercutsCalc(const char *hilname, 
-                                     const char *hilsrcname, const char *name, const char *title)
+                                     const char *hilsrcname, 
+                                     const char *name, const char *title)
   : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname),
-    fSuperName("MCT1Supercuts")
+    fHilExtName("MHillasExt"), fNewParName("MNewImagePar"), 
+    fSuperName("MCT1Supercuts") 
 {
     fName  = name  ? name  : "MCT1SupercutsCalc";
@@ -101 +104 @@
         return kFALSE;
     }
-
 
     if (fMatrix)
@@ -114 +116 @@
     }
 
+    fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt");
+    if (!fHilExt)
+    {
+        *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl;
+        return kFALSE;
+    }
+
     fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
     if (!fHilSrc)
     {
         *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar");
+    if (!fNewPar)
+    {
+        *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl;
         return kFALSE;
     }
@@ -172 +188 @@
 Double_t MCT1SupercutsCalc::GetVal(Int_t i) const
 {
-    return (*fMatrix)[fMap[i]];
+
+    Double_t val = (*fMatrix)[fMap[i]];
+
+
+    //*fLog << "MCT1SupercutsCalc::GetVal; i, fMatrix, fMap, val = "
+    //    << i << ",  " << fMatrix << ",  " << fMap[i] << ",  " << val << endl;
+
+
+    return val;
 }
 
@@ -187 +211 @@
 {
     if (fMatrix)
-        return;
+      return;
 
     fMatrix = mat;
@@ -199 +223 @@
     fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
     fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)");
+    fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
+    fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc");
+    fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1");
 }
 
 // ---------------------------------------------------------------------------
 //
-// Evaluate dynamical supercuts for CT1 Mkn421 2001 (Daniel Kranich)
-// optimized for mkn 421 2001 data
+// Evaluate dynamical supercuts 
 // 
 //          set hadronness to 0.25 if cuts are fullfilled
@@ -222 +248 @@
     const Double_t dist0   = fMatrix ? GetVal(6) : fHilSrc->GetDist();
 
+    Double_t help;
+    if (!fMatrix)
+      help = TMath::Sign(fHilExt->GetM3Long(), 
+                         fHilSrc->GetCosDeltaAlpha());
+    const Double_t asym0   = fMatrix ? GetVal(8) : help;
+    const Double_t conc    = fMatrix ? GetVal(9) : fNewPar->GetConc();
+    const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1();
+
     const Double_t newdist = dist0 * fMm2Deg;
 
@@ -236 +270 @@
     const Double_t length  = length0 * fMm2Deg;
     const Double_t width   = width0  * fMm2Deg;
-
-    if (newdist < 1.05                                                     &&
+    const Double_t asym    = asym0   * fMm2Deg;
+
+    /*
+    *fLog << "newdist, length, width, asym, dist, conc, leakage = " 
+          << newdist << ",  " << length << ",  " << width << ",  "
+          << asym    << ",  " << dist   << ",  " << conc  << ",  " << leakage
+          << endl;
+  
+    *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = " 
+          << CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetLengthUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetLengthLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetWidthUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetWidthLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetConcUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetConcLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetLeakage1Up(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetLeakage1Lo(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << endl;
+    */
+
+
+    if (
+        //dist    < 1.05                                                     &&
+        //newdist < 1.05                                                     &&
+
         newdist < CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) &&
         newdist > CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) &&
-        dist    < 1.05                                                     &&
+
         length  < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) &&
         length  > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) &&
+
         width   < CtsMCut (fSuper->GetWidthUp(),  dmls, dmcza, dmls2, dd2) &&
         width   > CtsMCut (fSuper->GetWidthLo(),  dmls, dmcza, dmls2, dd2) &&
-        //asym  < CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2) &&
-        //asym  > CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2) &&
+
+        asym    < CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2) &&
+        asym    > CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2) &&
+
         dist    < CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) &&
-        dist    > CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2)  )
-        fHadronness->SetHadronness(0.25);
+        dist    > CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) &&
+
+        conc    < CtsMCut (fSuper->GetConcUp(),   dmls, dmcza, dmls2, dd2) &&
+        conc    > CtsMCut (fSuper->GetConcLo(),   dmls, dmcza, dmls2, dd2) &&
+
+        leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) &&
+        leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2)  ) 
+
+      fHadronness->SetHadronness(0.25);
     else
-        fHadronness->SetHadronness(0.75);
+      fHadronness->SetHadronness(0.75);
+
+    //*fLog << "SChadroness = " << fHadronness->GetHadronness() << endl;
 
     fHadronness->SetReadyToSave();
@@ -257 +338 @@
     return kTRUE;
 }
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h

@@ -13 +13 @@
 class MHillas;
 class MHillasSrc;
+class MHillasExt;
+class MNewImagePar;
 class MMcEvt;
 class MCerPhotEvt;
@@ -25 +27 @@
     MHillas       *fHil;
     MHillasSrc    *fHilSrc;
+    MHillasExt    *fHilExt;
+    MNewImagePar  *fNewPar;
     MMcEvt        *fMcEvt;
     MHadronness   *fHadronness; //! output container for hadronness
@@ -32 +36 @@
     TString  fHilName;
     TString  fHilSrcName;
+    TString  fHilExtName;
+    TString  fNewParName;
     TString  fSuperName;        // name of container for supercut parameters
 
     Double_t fMm2Deg;           //!
 
-    Int_t     fMap[8];          //!
+    Int_t     fMap[11];         //!
     MHMatrix *fMatrix;          //!
 
@@ -63 +69 @@
 
 #endif
+
+
+
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc

@@ -103 +103 @@
     //  3   maximum output, showing progress of minimizations
     //
-    minuit.SetPrintLevel(1);
+    minuit.SetPrintLevel(-1);
 
     //..............................................
@@ -175 +175 @@
 
     //..............................................
+    // This doesn't seem to have any effect
     // Set maximum number of iterations (default = 500)
     //Int_t maxiter = 100000;
-    minuit.SetMaxIterations(10);
+    //minuit.SetMaxIterations(maxiter);
 
     //..............................................
@@ -211 +212 @@
         if (nulloutput)
             fLog->SetNullOutput(kTRUE);
-        Double_t tmp = 0;
-        minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize);
+        Int_t    maxcalls  = 3000;
+        Double_t tolerance = 0.1;
+        Double_t tmp[2];
+        tmp[0] = maxcalls;
+        tmp[1] = tolerance;
+        minuit.mnexcm("SIMPLEX", &tmp[0], 2, fErrMinimize);
         if (nulloutput)
             fLog->SetNullOutput(kFALSE);
-        //*fLog << "return from SIMPLEX" << endl;
+        *fLog << "return from SIMPLEX" << endl;
     }
 
@@ -269 +274 @@
     //            2    values, errors, step sizes, internal values
     //            3    values, errors, step sizes, 1st derivatives
-    //            4    values, paraboloc errors, MINOS errors
+    //            4    values, parabolic errors, MINOS errors
 
     //Int_t nkode = 4;
@@ -292 +297 @@
     return kTRUE;
 }
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/mfilter/MFEventSelector2.cc

@@ -188 +188 @@
     }
 
+    *fLog << "-------------------------" << endl;
+    *fLog << "MFEventSelector2::ReadDistribution; read input file to generate the original distribution" << endl;
+
     MEvtLoop run(GetName());
     MParList plist;
@@ -215 +218 @@
     }
 
+    tlist.PrintStatistics(0, kTRUE);
+
+    *fLog << "MFEventSelector2::ReadDistribution; Original distribution has " 
+          << fHistOrig->GetHist().GetEntries() << " entries" << endl;
+    *fLog << "-------------------------" << endl;
+
     return read.Rewind();
 }
@@ -259 +268 @@
     if (fNumMax>0)
     {
-        cout << "SCALE: " << fNumMax/hn.Integral() << endl;
-        cout << "SCALE: " << fNumMax << endl;
-        cout << "SCALE: " << hn.Integral() << endl;
+        *fLog << "MFEventSelector2::PrepareHistograms; SCALE : fNumMax = "
+              << fNumMax << ",  hn.Integral() = "      << hn.Integral()
+              << ",  fNumMax/hn.Integral() = " << fNumMax/hn.Integral() 
+              << endl;
+
         hn.Scale(fNumMax/hn.Integral());
     }
@@ -329 +340 @@
 Int_t MFEventSelector2::PreProcess(MParList *parlist)
 {
+    memset(fErrors, 0, sizeof(fErrors));
+
     MTaskList *tasklist = (MTaskList*)parlist->FindObject("MTaskList");
     if (!tasklist)
@@ -383 +396 @@
     const Double_t valz=fDataZ.GetValue();
 
+
     // get corresponding bin number
     const Int_t bin = fHistNom->FindFixBin(valx, valy, valz)-1;
@@ -389 +403 @@
     if (bin<0)
         return kTRUE;
+
 
     if (gRandom->Rndm()*fIs[bin]<=fNom[bin])
@@ -404 +419 @@
     fIs[bin]-=1;
 
+    //----------------------
+    Int_t rc;
+    if (fResult)
+    {
+        rc = 0;
+        fErrors[rc]++;
+        return kTRUE;
+    }
+    rc = 1;
+    fErrors[rc]++;
+    //----------------------
+
     return kTRUE;
 }
@@ -413 +440 @@
 Int_t MFEventSelector2::PostProcess()
 {
+    //---------------------------------
+    if (GetNumExecutions() != 0)
+    {
+      *fLog << inf << endl;
+      *fLog << GetDescriptor() << " execution statistics:" << endl;
+      *fLog << dec << setfill(' ');
+      *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
+            << (int)(fErrors[1]*100/GetNumExecutions()) 
+            << "%) Events not selected" << endl;
+
+      *fLog << " " << fErrors[0] << " (" 
+            << (int)(fErrors[0]*100/GetNumExecutions()) 
+            << "%) Events selected!" << endl;
+      *fLog << endl;
+    }
+    //---------------------------------
+
     if (!fCanvas || !fDisplay)
         return kTRUE;

trunk/MagicSoft/Mars/mfilter/MFEventSelector2.h

@@ -42 +42 @@
 
     Bool_t fResult;
+    Int_t fErrors[2];
 
     TH1   &InitHistogram(MH3* &hist);

trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc

@@ -61 +61 @@
     //
     fName  = name  ? name  : "MHCT1Supercuts";
-    fTitle = title ? title : "Container for supercuts";
+    fTitle = title ? title : "Container for histograms for the supercuts";
 
 

trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc

@@ -353 +353 @@
   fHistOrig = fhist;
 
-  fHist = (TH1*)fHistOrig->Clone("|alpha| plot");
+  fHist = (TH1*)fHistOrig->Clone();
+  fHist->SetName(fhist->GetName());
   if ( !fHist )
   {
@@ -362 +363 @@
   
   fHist->Sumw2();
-  fHist->SetNameTitle("Alpha", "alpha plot");
+  //fHist->SetNameTitle("Alpha", "alpha plot");
   fHist->SetXTitle("|alpha|  [\\circ]");
   fHist->SetYTitle("Counts");
+  fHist->UseCurrentStyle();
 
   fAlphamin = alphamin;
@@ -502 +504 @@
   fHistOrig = fhist;
 
-  fHist = (TH1*)fHistOrig->Clone("|alpha| plot");
+  fHist = (TH1*)fHistOrig->Clone();
+  fHist->SetName(fhist->GetName());
   fHist->Sumw2();
-  fHist->SetNameTitle("alpha", "alpha plot");
+  //fHist->SetNameTitle("alpha", "alpha plot");
   fHist->SetXTitle("|alpha|  [\\circ]");
   fHist->SetYTitle("Counts");
-
+  fHist->UseCurrentStyle();
 
   fAlphamin = alphamin;
@@ -754 +757 @@
 
   nrebin += 1;
-  //if (fHist)
-    delete fHist;
-    fHist = NULL;
+  TString histname = fHist->GetName();
+  delete fHist;
+  fHist = NULL;
 
   *fLog << "MHFindSignificance::FitPolynomial; rebin the |alpha| plot, grouping "
@@ -766 +769 @@
   fHist = new TH1F;
   fHist->Sumw2();
-  fHist->SetNameTitle("Rebinned", "Rebinned alpha plot");
-
+  fHist->SetNameTitle(histname, histname);
+  fHist->UseCurrentStyle();
 
   // do rebinning such that x0 remains a lower bin edge
@@ -911 +914 @@
       }
 
-      *fLog << "FitPolynomial : before CallMinuit()" << endl;
+      //*fLog << "FitPolynomial : before CallMinuit()" << endl;
 
       MMinuitInterface inter;
@@ -918 +921 @@
                                          kFALSE);
 
-      *fLog << "FitPolynomial : after CallMinuit()" << endl;
+      //*fLog << "FitPolynomial : after CallMinuit()" << endl;
 
       if (rc != 0)
@@ -1909 +1912 @@
   // get errors
 
+  /*
   Double_t eplus; 
   Double_t eminus; 
@@ -1943 +1947 @@
     }
   }  
+  */
 
   //----------------------------------------
+  /*
   *fLog << "Covariance matrix :" << endl;
   for (Int_t j=0; j<=fDegree; j++)
@@ -1966 +1972 @@
     *fLog << endl;
   }
+  */
 
   *fLog << "---------------------------" << endl;

trunk/MagicSoft/Mars/mhist/MHMatrix.h

@@ -30 +30 @@
     static const TString gsDefTitle; //! Default Title
 
-    Int_t   fNumRow;    //! Number of dimensions of histogram
+    Int_t   fNumRow;    // Number of dimensions of histogram
     Int_t   fRow;       //! Present row
     TMatrix fM;         // Matrix to be filled

trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc

@@ -219 +219 @@
   //
   *fLog << inf << "--------------------------------------" << endl;
-  *fLog << inf << "Start minimization" << endl;
+  *fLog << inf << "Start minimization for the energy estimator" << endl;
 
   gMinuit = new TMinuit(12);
@@ -316 +316 @@
   //    eest.Print();
   eest.StopMapping();
-  *fLog << inf << "Minimization finished" << endl;
+  *fLog << inf << "Minimization for the energy estimator finished" << endl;
 
 }