Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C
===================================================================
--- trunk/MagicSoft/Mars/macros/CT1Analysis.C	(revision 2065)
+++ 	(revision )
@@ -1,2347 +1,0 @@
-
-#include "CT1EEst.C"
-
-void InitBinnings(MParList *plist)
-    {
-      cout << "InitBinnings" << endl;
-
-      MBinning *binssize = new MBinning("BinningSize");
-      binssize->SetEdgesLog(50, 10, 1.0e5);
-      plist->AddToList(binssize);
-
-      MBinning *binsdistc = new MBinning("BinningDist");
-      binsdistc->SetEdges(50, 0, 1.4);
-      plist->AddToList(binsdistc);
-
-      MBinning *binswidth = new MBinning("BinningWidth");
-      binswidth->SetEdges(50, 0, 1.0);
-      plist->AddToList(binswidth);
-
-      MBinning *binslength = new MBinning("BinningLength");
-      binslength->SetEdges(50, 0, 1.0);
-      plist->AddToList(binslength);
-
-      MBinning *binsalpha = new MBinning("BinningAlpha");
-      binsalpha->SetEdges(100, -100, 100);
-      plist->AddToList(binsalpha);
-
-      MBinning *binsasym = new MBinning("BinningAsym");
-      binsasym->SetEdges(50, -1.5, 1.5);
-      plist->AddToList(binsasym);
-
-      MBinning *binsht = new MBinning("BinningHeadTail");
-      binsht->SetEdges(50, -1.5, 1.5);
-      plist->AddToList(binsht);
-
-      MBinning *binsm3l = new MBinning("BinningM3Long");
-      binsm3l->SetEdges(50, -1.5, 1.5);
-      plist->AddToList(binsm3l);
-
-      MBinning *binsm3t = new MBinning("BinningM3Trans");
-      binsm3t->SetEdges(50, -1.5, 1.5);
-      plist->AddToList(binsm3t);
-
-   
-      //.....
-      MBinning *binsb = new MBinning("BinningSigmabar");
-      binsb->SetEdges( 100,  0.0,  5.0);
-      plist->AddToList(binsb);
-
-      Int_t nbins = 8;
-      TArrayD xbins(nbins+1);
-      xbins[0] =  0.0;
-      xbins[1] =  7.5;
-      xbins[2] = 17.5; 
-      xbins[3] = 23.5;
-      xbins[4] = 29.5;
-      xbins[5] = 35.5;
-      xbins[6] = 42.0;
-      xbins[7] = 55.0;
-      xbins[8] = 68.0;
-
-      MBinning *binth = new MBinning("BinningTheta");
-      binth->SetEdges(xbins);
-      plist->AddToList(binth);
-
-      MBinning *binsdiff = new MBinning("BinningDiffsigma2");
-      binsdiff->SetEdges(100, -5.0, 20.0);
-      plist->AddToList(binsdiff);
-    }
-
-void DeleteBinnings(MParList *plist)
-    {
-      cout << "DeleteBinnings" << endl;
-
-      TObject *bin;
-
-      bin = plist->FindObject("BinningSize");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningDist");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningWidth");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningLength");
-      if (!bin) delete bin;
-
-      bin = plist->FindObject("BinningAlpha");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningAsym");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningHeadTail");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningM3Long");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningM3Trans");
-      if (!bin) delete bin;
-
-      //.....
-      bin = plist->FindObject("BinningSigmabar");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningTheta");
-      if (!bin) delete bin;
-      bin = plist->FindObject("BinningDiffsigma2");
-      if (!bin) delete bin;
-    }
-
-void CT1Analysis()
-{
-    //-----------------------------------------------
-    const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; 
-
-    //const char *onfile  = "~magican/ct1test/wittek/mkn421_on.preproc"; 
-    const char *onfile  = "~magican/ct1test/wittek/mkn421_00-01"; 
-
-    const char *mcfile  = "~magican/ct1test/wittek/mkn421_mc_pedrms_0.001.preproc";
-    //-----------------------------------------------
-
-    // path for input for Mars
-    TString inPath = "~magican/ct1test/wittek/marsoutput/optionB/";
-
-    // path for output from Mars
-    TString outPath = "~magican/ct1test/wittek/marsoutput/optionB/";
-
-
-
-    //************************************************************************
-
-    // Job A_ON : read ON data
-    //  - generate sigmabar vs. Theta plot; 
-    //  - write root file for ON data (ON1.root);
-    //  - generate hadron matrix for g/h separation
-
-    Bool_t JobA_ON = kFALSE;  
-    Bool_t WHistON = kFALSE;   // write out histogram sigmabar vs. Theta ?
-    Bool_t WLookON = kFALSE;   // write out look-alike events for hadrons ?
-    Bool_t WON1    = kFALSE;   // write out root file ON1.root ?
-
-
-    // Job A_MC : read MC gamma data, 
-    //  - read sigmabar vs. Theta plot from ON data  
-    //  - do padding; 
-    //  - write root file for MC gammas (MC1.root);
-    //  - generate gamma matrix for g/h separation)
-
-    Bool_t JobA_MC  = kFALSE;  
-    Bool_t WLookMC  = kFALSE;  // write out look-alike events for gammas ?
-    Bool_t WMC1     = kFALSE;  // write out root file MC1.root ?
-
-
-    // Job B_NN_UP : read ON1.root (or MC1.root) file 
-    //  - read hadron and gamma matrix
-    //  - calculate hadroness for method of Nearest Neighbors
-    //    and for the supercuts
-    //  - update the input files with the hadronesses (ON2.root or MC2.root)
-
-    Bool_t JobB_NN_UP  = kFALSE;  
-    Bool_t WNN         = kFALSE;  // write root file ?
-
-
-    // Job B_EST_UP : 
-    //  - read MC2.root file 
-    //  - select g/h separation method XX
-    //  - optimize energy estimation for events passing the final cuts
-    //  - write parameters of energy estimator onto file
-    //  - update ON2.root and MC2.root files with estimated energy
-    //    (ON_XX2.root and MC_XX2.root)
-
-    Bool_t JobB_EST_UP  = kFALSE;  
-    Bool_t WESTUP       = kFALSE;  // update root file ?
-
-
-    // Job C_NN : read ON1.root and MC1.root files  
-    //  - read look-alike events for hadrons and gammas 
-    //  - calculate hadronness for method of NEAREST NEIGHBOR
-    //  - produce Neyman-Pearson plot
-
-    Bool_t JobC_NN  = kFALSE;  
-
-
-    // Job C_RF : read ON1.root and MC1.root files  
-    //  - read look-alike events for hadrons and gammas 
-    //  - calculate hadronness for method of RANDOM FOREST
-    //  - produce Neyman-Pearson plot
-
-    Bool_t JobC_RF  = kFALSE;  
-
-
-    // Job D_XX :  
-    //  - select g/h separation method XX
-    //  - read MC_XX2.root file 
-    //  - calculate eff. collection area
-    //  - read ON_XX2.root file 
-    //  - apply final cuts
-    //  - calculate flux
-    //  - write root file for ON data after final cuts (ON3.root))
-    Bool_t JobD_XX  = kTRUE;  
-    Bool_t WXXD     = kFALSE;  // write out root file ON3.root ?
-
-
-    //************************************************************************
-
-    
-  //---------------------------------------------------------------------
-  // Job A_ON
-  //=========
-    // read ON data file 
-
-    //  - produce the 2D-histogram "sigmabar versus Theta" 
-    //    (SigmaTheta_ON.root) for ON data
-    //    (to be used for the padding of the MC gamma data)
-
-    //  - write a file of look-alike events (matrix_hadrons_ON.root)
-    //    (to be used later together with the corresponding MC gamma file
-    //     for the g/h separation in the analysis of the ON data)
-
-    //  - write a file of ON events (ON1.root) 
-    //    (after the standard cuts, before the g/h separation)
-    //    (to be used together with the corresponding MC gamma file (MC1.root)
-    //     for the optimization of the g/h separation)
-
-
- if (JobA_ON)
- {
-    cout << "=====================================================" << endl;
-    cout << "Macro CT1Analysis : Start of Job A_ON" << endl;
-    cout << "" << endl;
-    cout << "Macro CT1Analysis : JobA_ON, WHistON, WLookON, WON1 = " 
-         << JobA_ON  << ",  " << WHistON << ",  " << WLookON << ",  " 
-         << WON1 << endl;
-
-    TString typeInput = "ON";
-    cout << "typeInput = " << typeInput << endl;
-
-    // name of input root file
-    TString filenamein(onfile);
-
-    // name of output root file
-    TString outNameImage = outPath;
-    outNameImage += typeInput;
-    outNameImage += "1.root";
-
-
-    //-----------------------------------------------------------
-    MTaskList tliston;
-    MParList pliston;
-    MFilterList flist;
-
-    char *sourceName = "MSrcPosCam";
-    MSrcPosCam source(sourceName);
-
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MCT1ReadPreProc read(filenamein);
-
-    MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", 
-                                               "MCT1PointingCorrCalc");
-
-    MBlindPixelCalc blind;
-
-    MFCT1SelBasic selbasic;
-    MContinue contbasic(&selbasic);
-
-    MSigmabarCalc sigbarcalc;
-
-    MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
-    fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
-
-    MImgCleanStd    clean; 
-
-    // calculate also extended image parameters
-    TString fHilName = "MHillas";
-    MHillasExt hext(fHilName);
-
-    MHillasCalc    hcalc(fHilName);
-
-    TString fHilSrcName = "MHillasSrc";
-    MHillasSrc hilsrc(fHilSrcName);
-
-    MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", "");
-    csrc1.SetInput(fHilName);
-
-    TString fnewparName = "MNewImagePar";
-    MNewImagePar newpar(fnewparName);
-
-    MNewImageParCalc newimage(sourceName, fnewparName);
-    newimage.SetInput(fHilName);
-
-    MFCT1SelStandard selstandard(fHilName, fHilSrcName);
-    MContinue contstandard(&selstandard);
-
-    MFillH hfill1("MHHillas",    fHilName);
-    hfill1.SetTitle("Task to plot the source independent Hillas parameters");
-
-    MFillH hfill2("MHStarMap",   fHilName);
-    hfill2.SetTitle("Task to plot the star map (points along main axis of shower)");
-
-    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
-
-    MFillH hfill3("MHHillasExt",   fHilSrcName);
-    hfill3.SetTitle("Task to plot the extended Hillas parameters");
-
-    MFillH hfill4("MHHillasSrc",   fHilSrcName);
-    hfill4.SetTitle("Task to plot the source dependent Hillas parameters");
-
-    MFillH hfill5("MHNewImagePar", fnewparName);
-    hfill5.SetTitle("Task to plot the new image parameters");
-
-    //+++++++++++++++++++++++++++++++++++++++++++++++++++
-    // prepare writing of look-alike events (for hadrons)
-
-    const char* mtxName = "MatrixHadrons";
-    Int_t howMany = 500000;
-    TString outName = outPath;
-    outName += "matrix_hadrons_";
-    outName += typeInput;
-    outName += ".root";
-
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << "Matrix for (hadrons)" << endl;
-    cout << "input file                    = " << filenamein<< endl;
-    cout << "matrix name                   = " << mtxName << endl;
-    cout << "max. no.of look-alike events  = " << howMany << endl;
-    cout << "name of output root file      = " << outName << endl;
-    cout << "" << endl;
-
-
-    // matrix limitation for look-alike events (approximate number)
-    MFEventSelector selector("", "");
-    selector.SetNumSelectEvts(howMany);
-
-    //
-    // --- setup the matrix and add it to the parlist
-    //
-    MHMatrix *pmatrix = new MHMatrix(mtxName);
-    MHMatrix& matrix = *pmatrix;
-
-    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
-    matrix.AddColumn("MSigmabar.fSigmabar");
-    matrix.AddColumn("log10(MHillas.fSize)");
-    matrix.AddColumn("MHillasSrc.fDist");
-    matrix.AddColumn("MHillas.fWidth");
-    matrix.AddColumn("MHillas.fLength");
-    matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
-    matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)");
-    matrix.AddColumn("MHillas.fConc");
-    matrix.AddColumn("MNewImagePar.fLeakage1");
-
-    matrix.Print("cols");
-
-    MFillH fillmatg(mtxName);
-    fillmatg.SetFilter(&flist);
-    //+++++++++++++++++++++++++++++++++++++++++++++++++++
-
-    if (WON1)
-    {    
-      MWriteRootFile write(outNameImage);
-
-      write.AddContainer("MRawRunHeader", "RunHeaders");
-      write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      write.AddContainer("MSrcPosCam",    "Events");
-      write.AddContainer("MSigmabar",     "Events");
-      write.AddContainer("MHillas",       "Events");
-      write.AddContainer("MHillasSrc",    "Events");
-      write.AddContainer("MNewImagePar",  "Events");
-    }
-
-    //*****************************
-    // entries in MFilterList
-
-    flist.AddToList(&selector);
-
-
-    //*****************************
-    // entries in MParList
-    
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-    pliston.AddToList(&source);
-    pliston.AddToList(&sigbarcalc);
-
-    pliston.AddToList(&hext);
-    pliston.AddToList(&hilsrc);
-    pliston.AddToList(&newpar);
-    pliston.AddToList(&hhilext);
-
-    pliston.AddToList(pmatrix);
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-    tliston.AddToList(&pointcorr);
-    tliston.AddToList(&blind);
-
-    tliston.AddToList(&contbasic);
-    tliston.AddToList(&sigbarcalc);
-    tliston.AddToList(&fillsigtheta);
-    tliston.AddToList(&clean);
-
-    tliston.AddToList(&hcalc);
-    tliston.AddToList(&csrc1);
-    tliston.AddToList(&newimage);
-
-    tliston.AddToList(&flist);
-    tliston.AddToList(&fillmatg);    
-
-    tliston.AddToList(&hfill1);
-    tliston.AddToList(&hfill2);
-    tliston.AddToList(&hfill3);
-    tliston.AddToList(&hfill4);
-    tliston.AddToList(&hfill5);
-
-    tliston.AddToList(&contstandard);
-    if (WON1)
-      tliston.AddToList(&write);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-    //evtloop.Write();
-
-    //Int_t maxevents = 1000000000;
-    Int_t maxevents = 20000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&pliston);
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-    TObject *c;
-    c = ((MHSigmaTheta*)pliston.FindObject("SigmaTheta", "MHSigmaTheta"))
-                                                            ->DrawClone();
-
-    c = pliston.FindObject("MHHillas")->DrawClone();
-    c = pliston.FindObject("MHHillasExt")->DrawClone();
-    c = pliston.FindObject("MHHillasSrc")->DrawClone();
-    c = pliston.FindObject("MHNewImagePar")->DrawClone();
-    c = pliston.FindObject("MHStarMap")->DrawClone();
-
-
-    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-    //
-    // ----------------------------------------------------------
-    //  Definition of the reference sample of look-alike events
-    //  (this is a subsample of the original sample)
-    // ----------------------------------------------------------
-    //
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    // Select a maximum of nmaxevts events from the sample of look-alike 
-    // events. They will form the reference sample.
-    Int_t nmaxevts  = 2000;
-
-    // target distribution for the variable in column refcolumn (start at 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->SetXTitle("cos( \\Theta)");
-    thsh->SetYTitle("Counts");
-    Float_t dbin = (tobin-frombin)/nbins;
-    Float_t lbin = frombin +dbin*0.5;
-    for (Int_t j=1; j<=nbins; j++) 
-    {
-      thsh->Fill(lbin,1.0);
-      lbin += dbin;
-    }
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << "Macro CT1Analysis : define reference sample for hadrons" << endl; 
-    cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
-         << refcolumn << ",  " << nmaxevts << endl;    
-
-    matrix.EnableGraphicalOutput();
-    Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts);
-
-    if ( !matrixok ) 
-    {
-      cout << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
-      return;
-    }
-    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-    //-------------------------------------------
-    // write out look-alike events (for hadrons)
-    //
-    if (WLookON)
-    {
-      TFile *writejens = new TFile(outName, "RECREATE", "");
-      matrix.Write();
-
-      cout << "" << endl;
-      cout << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
-           << outName << endl;
-
-      writejens->Close();
-      delete writejens;
-    }
-
-
-    //-------------------------------------------
-    // Write histograms onto a file
-  if (WHistON)
-  {
-      MHSigmaTheta *sigtheta = 
-            (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
-      if (!sigtheta)
-	{
-          cout << "Object 'SigmaTheta' not found" << endl;
-          return;
-	}
-      TH2D *fHSigTh    = sigtheta->GetSigmaTheta();
-      TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta();
-      TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
-
-      TString outNameSigTh = outPath;
-      outNameSigTh += "SigmaTheta_";
-      outNameSigTh += typeInput;
-      outNameSigTh += ".root";
-
-      TFile outfile(outNameSigTh, "RECREATE");
-      fHSigTh->Write();
-      fHSigPixTh->Write();
-      fHDifPixTh->Write();
-      outfile.Close();
-     
-      cout << "" << endl;
-      cout << "File " << outNameSigTh << " was written out" << endl;
-  }
-
-
-    cout << "Macro CT1Analysis : End of Job A_ON" << endl;
-    cout << "===================================================" << endl;
- }
-
-
-  //---------------------------------------------------------------------
-   // Job A_MC
-   //=========
-
-    // read MC gamma data  
-
-    //    - to pad them
-    //      (using the 2D-histogram "sigmabar versus Theta" 
-    //       (SigmaTheta_ON.root)  of the ON data)
-
-    //    - to write a file of look-alike events (matrix_gammas_MC.root)
-    //      (to be used later together with the corresponding hadron file
-    //       for the g/h separation in the analysis of the ON data)
-
-    //    - to write a file of padded MC gamma events (MC1.root)
-    //      (after the standard cuts, before the g/h separation)
-    //      (to be used together with the corresponding hadron file
-    //       for the optimization of the g/h separation)
-
-
- if (JobA_MC)
- {
-    cout << "=====================================================" << endl;
-    cout << "Macro CT1Analysis : Start of Job A_MC" << endl;
-
-    cout << "" << endl;
-    cout << "Macro CT1Analysis : JobA_MC, WLookMC, WMC1 = " 
-         << JobA_MC  << ",  " << WLookMC << ",  " << WMC1 << endl;
-
-    TString typeInput = "MC";
-    cout << "typeInput = " << typeInput << endl;
-
-
-    // name of input root file
-    TString filenamein(mcfile);
-
-    // name of output root file
-    TString outNameImage = outPath;
-    outNameImage += typeInput;
-    outNameImage += "1.root";
-
-
-    // use for padding sigmabar vs. Theta from ON data
-    TString typeHist = "ON";
-    cout << "typeHist = " << typeHist << endl;
-
-    //------------------------------------
-    // Get the histograms "2D-ThetaSigmabar"
-    // and                "3D-ThetaPixSigma"
-    // and                "3D-ThetaPixDiff"
-
-
-      TString outNameSigTh = outPath;
-      outNameSigTh += "SigmaTheta_";
-      outNameSigTh += typeHist;
-      outNameSigTh += ".root";
-
-
-      cout << "Reading in file " << outNameSigTh << endl;
-
-      TFile *infile = new TFile(outNameSigTh);
-      infile->ls();
-
-      TH2D *fHSigmaTheta = 
-      (TH2D*) gROOT.FindObject("2D-ThetaSigmabar");
-      if (!fHSigmaTheta)
-	{
-          cout << "Object '2D-ThetaSigmabar' not found on root file" << endl;
-          return;
-	}
-      cout << "Object '2D-ThetaSigmabar' was read in" << endl;
-
-      TH3D *fHSigmaPixTheta = 
-      (TH3D*) gROOT.FindObject("3D-ThetaPixSigma");
-      if (!fHSigmaPixTheta)
-	{
-          cout << "Object '3D-ThetaPixSigma' not found on root file" << endl;
-          return;
-	}
-      cout << "Object '3D-ThetaPixSigma' was read in" << endl;
-
-      TH3D *fHDiffPixTheta = 
-      (TH3D*) gROOT.FindObject("3D-ThetaPixDiff");
-      if (!fHSigmaPixTheta)
-	{
-          cout << "Object '3D-ThetaPixDiff' not found on root file" << endl;
-          return;
-	}
-      cout << "Object '3D-ThetaPixDiff' was read in" << endl;
-
-    //------------------------------------
-
-    MTaskList tlist;
-    MParList plist;
-    MFilterList flist;
-
-    char *sourceName = "MSrcPosCam";
-    MSrcPosCam source(sourceName);
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MCT1ReadPreProc read(filenamein);
-
-    MBlindPixelCalc blind;
-
-    MFCT1SelBasic selbasic;
-    MContinue contbasic(&selbasic);
-
-
-    // There are 2 options for Thomas Schweizer's padding
-    //     fPadFlag = 1   get Sigmabar from fHSigmaTheta
-    //                    and Sigma    from fHDiffPixTheta
-    //     fPadFlag = 2   get Sigma    from fHSigmaPixTheta
-    
-    MPadSchweizer padthomas("MPadSchweizer","Task for the padding (Schweizer)",
-                             fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta);
-    padthomas.SetPadFlag(1);
-
-    //...........................................
-
-    MSigmabarCalc sigbarcalc;
-
-    MFillH fillsigtheta ("MCSigmaTheta[MHSigmaTheta]", "MMcEvt");
-    fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
-
-    MImgCleanStd    clean; 
-
-    // calculate also extended image parameters
-    TString fHilName = "MHillas";
-    MHillasExt hext(fHilName);
-
-    MHillasCalc    hcalc(fHilName);
-
-    TString fHilSrcName = "MHillasSrc";
-    MHillasSrc hilsrc(fHilSrcName);
-
-    MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", "");
-    csrc1.SetInput(fHilName);
-
-    TString fnewparName = "MNewImagePar";
-    MNewImagePar newpar(fnewparName);
-
-    MNewImageParCalc newimage(sourceName, fnewparName);
-    newimage.SetInput(fHilName);
-
-    MFCT1SelStandard selstandard(fHilName, fHilSrcName);
-    MContinue contstandard(&selstandard);
-
-    MFillH hfill1("MHHillas",    fHilName);
-    hfill1.SetTitle("Task to plot the source independent Hillas parameters");
-
-    MFillH hfill2("MHStarMap",   fHilName);
-    hfill2.SetTitle("Task to plot the star map (points along main axis of shower)");
-
-    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
-
-    MFillH hfill3("MHHillasExt",   fHilSrcName);
-    hfill3.SetTitle("Task to plot the extended Hillas parameters");
-
-    MFillH hfill4("MHHillasSrc",   fHilSrcName);
-    hfill4.SetTitle("Task to plot the source dependent Hillas parameters");
-
-    MFillH hfill5("MHNewImagePar", fnewparName);
-    hfill5.SetTitle("Task to plot the new image parameters");
-
-
-    //+++++++++++++++++++++++++++++++++++++++++++++++++++
-    // prepare writing of look-alike events (for gammas)
-
-    const char* mtxName = "MatrixGammas";
-    Int_t howMany = 500000;
-
-    TString outName = outPath;
-    outName += "matrix_gammas_";
-    outName += typeInput;
-    outName += ".root";
-
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << "Matrix for (gammas)" << endl;
-    cout << "input file                    = " << filenamein<< endl;
-    cout << "matrix name                   = " << mtxName << endl;
-    cout << "max. no.of look-alike events  = " << howMany << endl;
-    cout << "name of output root file      = " << outName << endl;
-    cout << "" << endl;
-
-
-    // matrix limitation for look-alike events (approximate number)
-    MFEventSelector selector("", "");
-    selector.SetNumSelectEvts(howMany);
-
-    //
-    // --- setup the matrix and add it to the parlist
-    //
-    MHMatrix *pmatrix = new MHMatrix(mtxName);
-    MHMatrix& matrix = *pmatrix;
-
-    matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
-    matrix.AddColumn("MSigmabar.fSigmabar");
-    matrix.AddColumn("log10(MHillas.fSize)");
-    matrix.AddColumn("MHillasSrc.fDist");
-    matrix.AddColumn("MHillas.fWidth");
-    matrix.AddColumn("MHillas.fLength");
-    matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
-    matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)");
-    matrix.AddColumn("MHillas.fConc");
-    matrix.AddColumn("MNewImagePar.fLeakage1");
-
-    matrix.Print("cols");
-
-    MFillH fillmatg(mtxName);
-    fillmatg.SetFilter(&flist);
-    //+++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-    if (WMC1)
-    {    
-      MWriteRootFile write(outNameImage);
-
-      write.AddContainer("MRawRunHeader", "RunHeaders");
-      write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      write.AddContainer("MSrcPosCam",    "Events");
-      write.AddContainer("MSigmabar",     "Events");
-      write.AddContainer("MHillas",       "Events");
-      write.AddContainer("MHillasSrc",    "Events");
-      write.AddContainer("MNewImagePar",  "Events");
-    }
-
-    //*****************************
-    // entries in MFilterList
-
-    flist.AddToList(&selector);
-
-
-    //*****************************
-    // entries in MParList
-
-    plist.AddToList(&tlist);
-    InitBinnings(&plist);
-
-    plist.AddToList(&source);
-    plist.AddToList(&sigbarcalc);
-
-    plist.AddToList(&hext);
-    plist.AddToList(&hilsrc);
-    plist.AddToList(&newpar);
-    plist.AddToList(&hhilext);
-
-    plist.AddToList(pmatrix); 
-
-
-    //*****************************
-    // entries in MTaskList
-
-    tlist.AddToList(&read);
-    tlist.AddToList(&blind);
-    tlist.AddToList(&padthomas);
-
-    tlist.AddToList(&contbasic);
-    tlist.AddToList(&sigbarcalc);
-    tlist.AddToList(&fillsigtheta);
-    tlist.AddToList(&clean);
-
-    tlist.AddToList(&hcalc);
-    tlist.AddToList(&csrc1);
-    tlist.AddToList(&newimage);
-
-    tlist.AddToList(&flist);
-    tlist.AddToList(&fillmatg);
-
-    tlist.AddToList(&hfill1);
-    tlist.AddToList(&hfill2);
-    tlist.AddToList(&hfill3);
-    tlist.AddToList(&hfill4);
-    tlist.AddToList(&hfill5);
-
-    tlist.AddToList(&contstandard);
-    if (WMC1)
-      tlist.AddToList(&write);
-
-    //*****************************
-
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&plist);
-    evtloop.SetProgressBar(&bar);
-    //evtloop.Write();
-
-    //Int_t maxevents = 1000000000;
-    Int_t maxevents = 1000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tlist.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&plist);
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-    TObject *c;
-    c = ((MHSigmaTheta*)plist.FindObject("MCSigmaTheta", "MHSigmaTheta"))
-                                                            ->DrawClone();
-    c = plist.FindObject("MHHillas")->DrawClone();
-    c = plist.FindObject("MHHillasExt")->DrawClone();
-    c = plist.FindObject("MHHillasSrc")->DrawClone();
-    c = plist.FindObject("MHNewImagePar")->DrawClone();
-    c = plist.FindObject("MHStarMap")->DrawClone();
-
-
-    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-    //
-    // ----------------------------------------------------------
-    //  Definition of the reference sample of look-alike events
-    //  (this is a subsample of the original sample)
-    // ----------------------------------------------------------
-    //
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    // Select a maximum of nmaxevts events from the sample of look-alike 
-    // events. They will form the reference sample.
-    Int_t nmaxevts  = 2000;
-
-    // target distribution for the variable in column refcolumn (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);
-    Float_t dbin = (tobin-frombin)/nbins;
-    Float_t lbin = frombin +dbin*0.5;
-    for (Int_t j=1; j<=nbins; j++) 
-    {
-      thsh->Fill(lbin,1.0);
-      lbin += dbin;
-    }
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << "Macro CT1Analysis : define reference sample for gammas" << endl; 
-    cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
-         << refcolumn << ",  " << nmaxevts << endl;    
-
-    matrix.EnableGraphicalOutput();
-    Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts);
-    if ( !matrixok ) 
-    {
-      cout << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
-      return;
-    }
-    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-    //-------------------------------------------
-    // write out look-alike events (for gammas)
-    //
-    if (WLookMC)
-    {
-      TFile *writejens = new TFile(outName, "RECREATE", "");
-      matrix.Write();
-      cout << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
-           << outName << endl;
-
-      writejens->Close();
-      delete writejens;
-    }
-
-    cout << "Macro CT1Analysis : End of Job A_MC" 
-         << endl;
-    cout << "=========================================================" 
-         << endl;
- }
-
-
-
-  //---------------------------------------------------------------------
-  // Job B_NN_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)
- {
-    cout << "=====================================================" << endl;
-    cout << "Macro CT1Analysis : Start of Job B_NN_UP" << endl;
-
-    cout << "" << endl;
-    cout << "Macro CT1Analysis : JobB_NN_UP, WNN = " 
-         << JobB_NN_UP  << ",  " << WNN << endl;
-
-
-    //TString typeInput = "ON";
-    TString typeInput = "MC";
-    cout << "typeInput = " << typeInput << endl;
-
-    // name of input root file
-    TString filenameData = outPath;
-    filenameData += typeInput;
-    filenameData += "1.root";
-    cout << "filenameData = " << filenameData << endl; 
-
-    // name of output root file
-    TString outNameImage = outPath;
-    outNameImage += typeInput;
-    outNameImage += "2.root";
-
-    //$$$$$$$$$$$$$$$$$$$$$$$$$$
-    outNameImage = filenameData;
-    //$$$$$$$$$$$$$$$$$$$$$$$$$$
-
-    cout << "outNameImage = " << outNameImage << endl; 
-
-    TString typeMatrixHadrons = "ON";
-    cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
-
-   //*************************************************************************
-   // read in matrices of look-alike events
-
-    const char* mtxName = "MatrixGammas";
-
-    TString outName = outPath;
-    outName += "matrix_gammas_";
-    outName += "MC";
-    outName += ".root";
-
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << "Get matrix for (gammas)" << endl;
-    cout << "matrix name        = " << mtxName << endl;
-    cout << "name of root file  = " << outName << endl;
-    cout << "" << endl;
-
-
-    // read in the object with the name 'mtxName' from file 'outName'
-    //
-    TFile *fileg = new TFile(outName); 
-
-    MHMatrix matrixg;
-    matrixg.Read(mtxName);
-    matrixg.Print("cols");
-
-    delete fileg;
-
-    //***************************************************************** 
-
-    const char* mtxName = "MatrixHadrons";
-
-    TString outName = outPath;
-    outName += "matrix_hadrons_";
-    outName += typeMatrixHadrons;
-    outName += ".root";
-
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << " Get matrix for (hadrons)" << endl;
-    cout << "matrix name        = " << mtxName << endl;
-    cout << "name of root file  = " << outName << endl;
-    cout << "" << endl;
-
-
-    // read in the object with the name mtxName
-    //
-    TFile *fileh = new TFile(outName); 
-
-    MHMatrix matrixh;
-    matrixh.Read(mtxName);
-    matrixh.Print("cols");
-
-    delete fileh;
-
-   //*************************************************************************
-
-
-    MTaskList tliston;
-    MParList pliston;
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameData);
-    read.DisableAutoScheme();
-
-    TString fHilName    = "MHillas"; 
-    TString fHilSrcName = "MHillasSrc"; 
-    TString fnewparName = "MNewImagePar"; 
-
-
-    //.......................................................................
-    // 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, fHilSrcName);
-    sccalc.SetHadronnessName(hadSCName);
-
-    //.......................................................................
-
-    if (WNN)
-    {    
-      //MWriteRootFile write(outNameImage);
-      MWriteRootFile write(outNameImage, "UPDATE");
-
-      //write.AddContainer("MRawRunHeader", "RunHeaders");
-      //write.AddContainer("MTime",         "Events");
-      //write.AddContainer("MMcEvt",        "Events");
-      //write.AddContainer("MSrcPosCam",    "Events");
-      //write.AddContainer("MSigmabar",     "Events");
-      //write.AddContainer("MHillas",       "Events");
-      //write.AddContainer("MHillasSrc",    "Events");
-      //write.AddContainer("MNewImagePar",  "Events");
-
-      write.AddContainer(hadNNName,       "Events");
-      write.AddContainer(hadSCName,       "Events");
-    }
-
-    //-----------------------------------------------------------------
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-
-    Float_t hadcut   =  0.40;
-    Float_t alphacut = 100.0;
-    MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
-    selfinalgh.SetCuts(hadcut, alphacut);
-    selfinalgh.SetHadronnessName(hadNNName);
-    MContinue contfinalgh(&selfinalgh);
-
-    MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
-    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
-
-    Float_t hadcut   = 0.40;
-    Float_t alphacut = 20.0;
-    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
-    selfinal.SetCuts(hadcut, alphacut);
-    selfinal.SetHadronnessName(hadNNName);
-    MContinue contfinal(&selfinal);
-
-
-    MFillH hfill1("MHHillas",    fHilName);
-    MFillH hfill2("MHStarMap",   fHilName);
-
-    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
-
-    MFillH hfill3("MHHillasExt",   fHilSrcName);
-    MFillH hfill4("MHHillasSrc",   fHilSrcName);
-    MFillH hfill5("MHNewImagePar", fnewparName);
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-    pliston.AddToList(&matrixg);
-    pliston.AddToList(&matrixh);
-
-    pliston.AddToList(&hhilext); 
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-
-    tliston.AddToList(&nncalc);
-    tliston.AddToList(&sccalc);
-    tliston.AddToList(&fillhadnn);
-    tliston.AddToList(&fillhadsc);
-
-    if (WNN)
-      tliston.AddToList(&write);
-
-    tliston.AddToList(&contfinalgh);
-    tliston.AddToList(&hfill1);
-    tliston.AddToList(&hfill2);
-    tliston.AddToList(&hfill3);
-    tliston.AddToList(&hfill4);
-    tliston.AddToList(&hfill5);
-
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-    //evtloop.Write();
-
-    if ( !evtloop.Eventloop() )
-      //if ( !evtloop.Eventloop(1000) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&pliston);
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-    TObject *c;
-
-    c = pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
-    c = pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
-
-    c = pliston.FindObject("MHHillas")->DrawClone();
-    c = pliston.FindObject("MHHillasExt")->DrawClone();
-    c = pliston.FindObject("MHHillasSrc")->DrawClone();
-    c = pliston.FindObject("MHNewImagePar")->DrawClone();
-    c = pliston.FindObject("MHStarMap")->DrawClone();
-
-
-    cout << "Macro CT1Analysis : End of Job B_NN_UP" << endl;
-    cout << "=======================================================" << endl;
- }
-  //---------------------------------------------------------------------
-
-
-
-  //---------------------------------------------------------------------
-  // Job B_EST_UP
-  //============
-
-    //  - read MC2.root file 
-    //  - select g/h separation method XX
-    //  - optimize energy estimator for events passing final cuts
-    //  - write parameters of energy estimator onto file "energyest_XX.root"
-    //
-    //  - read ON2.root and MC2.root files 
-    //  - update input root file with the estimated energy
-    //    (ON_XX2.root, MC_XX2.root)
-
-
- if (JobB_EST_UP)
- {
-    cout << "=====================================================" << endl;
-    cout << "Macro CT1Analysis : Start of Job B_EST_UP" << endl;
-
-    cout << "" << endl;
-    cout << "Macro CT1Analysis : JobB_EST_UP, WESTUP = " 
-         << JobB_EST_UP  << ",  " << WESTUP << endl;
-
-
-    TString typeON = "ON";
-    TString typeMC = "MC";
-    TString ext    = "3.root";
-
-    //------------------------------
-    // name of MC file to be used for optimizing the energy estimator
-    TString filenameOpt(outPath);
-    filenameOpt += typeMC;
-    filenameOpt += ext; 
-    cout << "filenameOpt = " << filenameOpt << endl;
-
-    //------------------------------
-    // selection of g/h separation method
-    // and definition of final selections
-
-    TString XX("NN");
-    //TString XX("SC");
-    //TString XX("RF");
-    TString fhadronnessName("Had");
-    fhadronnessName += XX;
-    cout << "fhadronnessName = " << fhadronnessName << endl;
-
-    // maximum values of the hadronness and of |alpha|
-    Float_t maxhadronness   = 0.40;
-    Float_t maxalpha        = 20.0;
-    cout << "Maximum values of hadronness and |alpha| = "
-         << maxhadronness << ",  " << maxalpha << endl;
-
-    // name of file containing the parameters of the energy estimator
-    TString energyParName(outPath);
-    energyParName += "energyest_";
-    energyParName += XX;
-    energyParName += ".root";
-    cout << "energyParName = " << energyParName << endl;
-
-
-    //------------------------------
-    // name of ON file to be updated
-    TString filenameON(outPath);
-    filenameON += typeON;
-    filenameON += ext;
-    cout << "filenameON = " << filenameON << endl;
-
-    // name of MC file to be updated
-    TString filenameMC(outPath);
-    filenameMC += typeMC;
-    filenameMC += ext;
-    cout << "filenameMC = " << filenameMC << endl;
-
-    //------------------------------
-    // name of updated ON file 
-    TString filenameONup(outPath);
-    filenameONup += typeON;
-    filenameONup += "_";
-    filenameONup += XX;
-    filenameONup += ext;
-    cout << "filenameONup = " << filenameONup << endl;
-
-    // name of updated MC file 
-    TString filenameMCup(outPath);
-    filenameMCup += typeMC;
-    filenameMCup += "_";
-    filenameMCup += XX;
-    filenameMCup += ext;
-    cout << "filenameMCup = " << filenameMCup << endl;
-
-    //-----------------------------------------------------------
-
-    TString fHilName    = "MHillas"; 
-    TString fHilSrcName = "MHillasSrc"; 
-    TString fnewparName = "MNewImagePar"; 
-
-    //===========================================================
-    //
-    // Optimization of energy estimator
-    //
-
-    TString inpath("");
-    TString outpath("");
-    Int_t howMany = 2000;
-    CT1EEst(inpath,   filenameOpt,   outpath, energyParName, 
-            fHilName, fHilSrcName,   fhadronnessName,
-            howMany,  maxhadronness, maxalpha);
-
-    //-----------------------------------------------------------
-    //
-    // Read in parameters of energy estimator
-    //
-
-    TFile enparam(energyParName);
- 
-    TArrayD parA(5);
-    TArrayD parB(7);
-    for (Int_t i=0; i<parA.GetSize(); i++)
-      parA[i] = EnergyParams(i);
-    for (Int_t i=0; i<parB.GetSize(); i++)
-      parB[i] = EnergyParams( i+parA.GetSize() );
-
-
-   if (WESTUP)
-   {
-    //==========   start update   ============================================
-    //
-    // Update ON and MC root files with the estimated energy
-
-    //---------------------------------------------------
-    // Update ON data
-    //
-    MTaskList tliston;
-    MParList pliston;
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameON);
-    read.DisableAutoScheme();
-
-    //---------------------------
-    // calculate estimated energy
-
-    MEnergyEstParam eest2(fHilName);
-    eest2.Add(fHilSrcName);
-
-    eest2.SetCoeffA(parA);
-    eest2.SetCoeffB(parB);
-
-    //.......................................................................
-
-      MWriteRootFile write(filenameONup);
-
-      write.AddContainer("MRawRunHeader", "RunHeaders");
-      write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      write.AddContainer("MSrcPosCam",    "Events");
-      write.AddContainer("MSigmabar",     "Events");
-      write.AddContainer("MHillas",       "Events");
-      write.AddContainer("MHillasSrc",    "Events");
-      write.AddContainer("MNewImagePar",  "Events");
-
-      write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
-
-      write.AddContainer("MEnergyEst",    "Events");
-
-    //-----------------------------------------------------------------
-
-    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
-    selfinal.SetCuts(maxhadronness, maxalpha);
-    selfinal.SetHadronnessName(fhadronnessName);
-    MContinue contfinal(&selfinal);
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-    tliston.AddToList(&eest2);
-    tliston.AddToList(&write);
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-    //evtloop.Write();
-
-    Int_t maxevents = 1000000000;
-    //Int_t maxevents = 1000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&pliston);
-
-    //---------------------------------------------------
-    //---------------------------------------------------
-    // Update MC data
-    //
-    MTaskList tliston;
-    MParList pliston;
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameMC);
-    read.DisableAutoScheme();
-
-    //---------------------------
-    // calculate estimated energy
-
-    MEnergyEstParam eest2(fHilName);
-    eest2.Add(fHilSrcName);
-
-    eest2.SetCoeffA(parA);
-    eest2.SetCoeffB(parB);
-
-    //.......................................................................
-
-      MWriteRootFile write(filenameMCup);
-
-      write.AddContainer("MRawRunHeader", "RunHeaders");
-      write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      write.AddContainer("MSrcPosCam",    "Events");
-      write.AddContainer("MSigmabar",     "Events");
-      write.AddContainer("MHillas",       "Events");
-      write.AddContainer("MHillasSrc",    "Events");
-      write.AddContainer("MNewImagePar",  "Events");
-
-      write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
-
-      write.AddContainer("MEnergyEst",    "Events");
-
-    //-----------------------------------------------------------------
-
-    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
-    selfinal.SetCuts(maxhadronness, maxalpha);
-    selfinal.SetHadronnessName(fhadronnessName);
-    MContinue contfinal(&selfinal);
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-    tliston.AddToList(&eest2);
-    tliston.AddToList(&write);
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-    //evtloop.Write();
-
-    Int_t maxevents = 1000000000;
-    //Int_t maxevents = 1000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&pliston);
-
-
-    //==========   end update   ============================================
-   }
-    
-    enparam.Close();
-
-    cout << "Macro CT1Analysis : End of Job B_EST_UP" << endl;
-    cout << "=======================================================" << endl;
- }
-  //---------------------------------------------------------------------
-
-
-  //---------------------------------------------------------------------
-  // Job C_NN
-  //=========
-
-    // read ON1 and MC1 data files  
-    //
-    //  - calculate hadronness for method of Nearest Neighbors
-    //  - produce Neyman-Pearson plot
- 
- if (JobC_NN)
- {
-    cout << "=====================================================" << endl;
-    cout << "Macro CT1Analysis : Start of Job C_NN" << endl;
-
-    cout << "" << endl;
-    cout << "Macro CT1Analysis : JobC_NN = " << JobC_NN  << endl;
-
-
-    TString typeMC = "MC";
-    cout << "typeMC = " << typeMC << endl;
-
-    TString typeData = "ON";
-    cout << "typeData = " << typeData << endl;
-
-    TString typeMatrixHadrons = "ON";
-    cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
-
-
-
-   //*************************************************************************
-   // read in matrices of look-alike events
-
-    const char* mtxName = "MatrixGammas";
-
-    TString outName = outPath;
-    outName += "matrix_gammas_";
-    outName += "MC";
-    outName += ".root";
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << "Get matrix for (gammas)" << endl;
-    cout << "matrix name        = " << mtxName << endl;
-    cout << "name of root file  = " << outName << endl;
-    cout << "" << endl;
-
-
-    // read in the object with the name 'mtxName' from file 'outName'
-    //
-    TFile *fileg = new TFile(outName); 
-
-    MHMatrix matrixg;
-    matrixg.Read(mtxName);
-    matrixg.Print("cols");
-
-    delete fileg;
-
-
-    //***************************************************************** 
-
-    const char* mtxName = "MatrixHadrons";
-
-    TString outName = outPath;
-    outName += "matrix_hadrons_";
-    outName += typeMatrixHadrons;
-    outName += ".root";
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << " Get matrix for (hadrons)" << endl;
-    cout << "matrix name        = " << mtxName << endl;
-    cout << "name of root file  = " << outName << endl;
-    cout << "" << endl;
-
-
-    // read in the object with the name mtxName
-    //
-    TFile *fileh = new TFile(outName); 
-
-    MHMatrix matrixh;
-    matrixh.Read(mtxName);
-    matrixh.Print("cols");
-
-    delete fileh;
-
-    //***************************************************************** 
-
-
-    //-----------------------------------------------------------------
-
-    MTaskList tliston;
-    MParList pliston;
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-    TString filenameData = outPath;
-    filenameData += typeData;
-    filenameData += "1.root";
-
-    TString filenameMC = outPath;
-    filenameMC += typeMC;
-    filenameMC += "1.root";
-   
-    cout << "filenameData = " << filenameData << endl;
-    cout << "filenameMC   = " << filenameMC   << endl;
-
-    MReadMarsFile read("Events", filenameMC);
-    read.AddFile(filenameData);
-    read.DisableAutoScheme();
-
-    MFParticleId fgamma ("MMcEvt", '=', kGAMMA);
-    MFParticleId fhadron("MMcEvt", '!', kGAMMA);
-
-    //.......................................................................
-    // calculate hadronnes for method of Nearest Neighbors
-
-    TString hadName = "HadNN";
-    MMultiDimDistCalc nncalc;
-    nncalc.SetUseNumRows(25);
-    nncalc.SetUseKernelMethod(kFALSE);
-    nncalc.SetHadronnessName(hadName);
-
-    //.......................................................................
-
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-    TString fHilName    = "MHillas"; 
-    TString fHilSrcName = "MHillasSrc"; 
-    TString fnewparName = "MNewImagePar"; 
-
-    Float_t hadcut   =  0.36;
-    Float_t alphacut = 100.0;
-    MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
-    selfinalgh.SetCuts(hadcut, alphacut);
-    selfinalgh.SetHadronnessName(hadName);
-    MContinue contfinalgh(&selfinalgh);
-
-    MFillH fillhadnn("MHHadronness", hadName);
-
-    Float_t hadcut   = 0.36;
-    Float_t alphacut = 20.0;
-    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
-    selfinal.SetCuts(hadcut, alphacut);
-    selfinal.SetHadronnessName(hadName);
-    MContinue contfinal(&selfinal);
-
-    MFillH hfill1("MHHillas",    fHilName);
-    MFillH hfill2("MHStarMap",   fHilName);
-
-    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
-
-    MFillH hfill3("MHHillasExt",   fHilSrcName);
-    MFillH hfill4("MHHillasSrc",   fHilSrcName);
-    MFillH hfill5("MHNewImagePar", fnewparName);
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-    pliston.AddToList(&matrixg);
-    pliston.AddToList(&matrixh);
-
-    pliston.AddToList(&hhilext); 
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-
-    tliston.AddToList(&fgamma);
-    tliston.AddToList(&fhadron);
-    tliston.AddToList(&nncalc);
-    tliston.AddToList(&fillhadnn);
-   
-    tliston.AddToList(&contfinalgh);
-    tliston.AddToList(&hfill1);
-    tliston.AddToList(&hfill2);
-    tliston.AddToList(&hfill3);
-    tliston.AddToList(&hfill4);
-    tliston.AddToList(&hfill5);
-
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-    //evtloop.Write();
-
-    Int_t maxevents = 1000000000;
-    //Int_t maxevents = 35000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&pliston);
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-    TObject *c;
-    c = pliston.FindObject("MHHadronness");
-    if (c)
-    {
-      c->DrawClone();
-      c->Print();
-    }
-    c = pliston.FindObject("MHHillas")->DrawClone();
-    c = pliston.FindObject("MHHillasExt")->DrawClone();
-    c = pliston.FindObject("MHHillasSrc")->DrawClone();
-    c = pliston.FindObject("MHNewImagePar")->DrawClone();
-    c = pliston.FindObject("MHStarMap")->DrawClone();
-
-
-
-    cout << "Macro CT1Analysis : End of Job C_NN" << endl;
-    cout << "===================================================" << endl;
- }
-
-
-
-  //---------------------------------------------------------------------
-  // Job C_RF
-  //=========
-
-    // read ON1 and MC1 data files  
-    //
-    //  - calculate hadronness for method of Random Forest
-    //  - produce Neyman-Pearson plot
- 
- if (JobC_RF)
- {
-    cout << "=====================================================" << endl;
-    cout << "Macro CT1Analysis : Start of Job C_RF" << endl;
-
-    cout << "" << endl;
-    cout << "Macro CT1Analysis : JobC_RF = " << JobC_RF  << endl;
-
-
-    TString typeMC = "MC";
-    cout << "typeMC = " << typeMC << endl;
-
-    TString typeData = "ON";
-    cout << "typeData = " << typeData << endl;
-
-    TString typeMatrixHadrons = "ON";
-    cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
-
-   //*************************************************************************
-   // read in matrices of look-alike events
-
-    const char* mtxName = "MatrixGammas";
-
-    TString outName = outPath;
-    outName += "matrix_gammas_";
-    outName += "MC";
-    outName += ".root";
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << "Get matrix for (gammas)" << endl;
-    cout << "matrix name        = " << mtxName << endl;
-    cout << "name of root file  = " << outName << endl;
-    cout << "" << endl;
-
-
-    // read in the object with the name 'mtxName' from file 'outName'
-    //
-    TFile *fileg = new TFile(outName); 
-
-    MHMatrix matrixg;
-    matrixg.Read(mtxName);
-    matrixg.Print("cols");
-
-    delete fileg;
-
-
-    //***************************************************************** 
-
-    const char* mtxName = "MatrixHadrons";
-
-    TString outName = outPath;
-    outName += "matrix_hadrons_";
-    outName += typeMatrixHadrons;
-    outName += ".root";
-
-    cout << "" << endl;
-    cout << "========================================================" << endl;
-    cout << " Get matrix for (hadrons)" << endl;
-    cout << "matrix name        = " << mtxName << endl;
-    cout << "name of root file  = " << outName << endl;
-    cout << "" << endl;
-
-
-    // read in the object with the name mtxName
-    //
-    TFile *fileh = new TFile(outName); 
-
-    MHMatrix matrixh;
-    matrixh.Read(mtxName);
-    matrixh.Print("cols");
-
-    delete fileh;
-
-    //***************************************************************** 
-
-    //-----------------------------------------------------------------
-    // grow the trees of the random forest (event loop = tree loop)
-
-    cout << "Macro CT1Analysis : start growing trees" << endl;
-
-    MTaskList tlist2;
-    MParList plist2;
-    plist2.AddToList(&tlist2);
-
-    plist2.AddToList(&matrixg);
-    plist2.AddToList(&matrixh);
-
-    MRanForestGrow rfgrow2;
-    rfgrow2.SetNumTrees(100);
-    rfgrow2.SetNumTry(4);
-    rfgrow2.SetNdSize(10);
-
-    TString outRF = outPath;
-    outRF += "RF.root";
-    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
-    MRanForest *fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
-    if (!fRanForest)
-    {
-        *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    MRanTree *fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
-    if (!fRanTree)                                  
-    {                                                                          
-        *fLog << err << dbginf << "MRanTree not found... aborting." << endl;    
-        return kFALSE;
-    }
-
-    //-----------------------------------------------------------------
-
-    MTaskList tliston;
-    MParList pliston;
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-    TString filenameData = outPath;
-    filenameData += typeData;
-    filenameData += "1.root";
-
-    TString filenameMC = outPath;
-    filenameMC += typeMC;
-    filenameMC += "1.root";
-   
-    cout << "filenameData = " << filenameData << endl;
-    cout << "filenameMC   = " << filenameMC   << endl;
-
-    MReadMarsFile read("Events", filenameMC);
-    read.AddFile(filenameData);
-    read.DisableAutoScheme();
-
-    MFParticleId fgamma ("MMcEvt", '=', kGAMMA);
-    MFParticleId fhadron("MMcEvt", '!', kGAMMA);
-
-    //.......................................................................
-    // calculate hadronnes for method of Random Forest
-
-    TString hadName = "HadRF";
-    MRanForestCalc rfcalc;
-    rfcalc.SetHadronnessName(hadName);
-
-    //.......................................................................
-
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-    TString fHilName    = "MHillas"; 
-    TString fHilSrcName = "MHillasSrc"; 
-    TString fnewparName = "MNewImagePar"; 
-
-    Float_t hadcut   =  0.36;
-    Float_t alphacut = 100.0;
-    MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
-    selfinalgh.SetCuts(hadcut, alphacut);
-    selfinalgh.SetHadronnessName(hadName);
-    MContinue contfinalgh(&selfinalgh);
-
-    MFillH fillhadrf("MHHadronness", hadName);
-    MFillH fillhrf("MHRanForest");
-
-    Float_t hadcut   = 0.36;
-    Float_t alphacut = 20.0;
-    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
-    selfinal.SetCuts(hadcut, alphacut);
-    selfinal.SetHadronnessName(hadName);
-    MContinue contfinal(&selfinal);
-
-    MFillH hfill1("MHHillas",    fHilName);
-    MFillH hfill2("MHStarMap",   fHilName);
-
-    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
-
-    MFillH hfill3("MHHillasExt",   fHilSrcName);
-    MFillH hfill4("MHHillasSrc",   fHilSrcName);
-    MFillH hfill5("MHNewImagePar", fnewparName);
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-    pliston.AddToList(fRanForest);
-    pliston.AddToList(fRanTree);
-
-    pliston.AddToList(&hhilext); 
-
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-
-    tliston.AddToList(&fgamma);
-    tliston.AddToList(&fhadron);
-    tliston.AddToList(&rfcalc);
-    tliston.AddToList(&fillhadrf);
-    tliston.AddToList(&fillhrf);
-   
-    tliston.AddToList(&contfinalgh);
-    tliston.AddToList(&hfill1);
-    tliston.AddToList(&hfill2);
-    tliston.AddToList(&hfill3);
-    tliston.AddToList(&hfill4);
-    tliston.AddToList(&hfill5);
-
-    tliston.AddToList(&contfinal);
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-    //evtloop.Write();
-
-    Int_t maxevents = 1000000000;
-    //Int_t maxevents = 35000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&pliston);
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-    TObject *c;
-    c = pliston.FindObject("MHRanForest")->DrawClone();    
-    c = pliston.FindObject("MHHadronness");
-    if (c)
-    {
-      c->DrawClone();
-      c->Print();
-    }
-
-    c = pliston.FindObject("MHHillas")->DrawClone();
-    c = pliston.FindObject("MHHillasExt")->DrawClone();
-    c = pliston.FindObject("MHHillasSrc")->DrawClone();
-    c = pliston.FindObject("MHNewImagePar")->DrawClone();
-    c = pliston.FindObject("MHStarMap")->DrawClone();
-
-
-
-    cout << "Macro CT1Analysis : End of Job C_RF" << endl;
-    cout << "===================================================" << endl;
- }
-
-
-  //---------------------------------------------------------------------
-  // Job D_XX
-  //=========
-
-    //  - select g/h separation method XX
-    //  - read MC_XX2.root file 
-    //  - calculate eff. collection area
-    //  - read ON_XX2.root file 
-    //  - apply final cuts
-    //  - calculate flux
-    //  - write root file for ON data after final cuts (ON_XX3.root))
-
-
- if (JobD_XX)
- {
-    cout << "=====================================================" << endl;
-    cout << "Macro CT1Analysis : Start of Job D_XX" << endl;
-
-    cout << "" << endl;
-    cout << "Macro CT1Analysis : JobD_XX, WXXD = " 
-         << JobD_XX  << ",  " << WXXD << endl;
-
-    // type of data to be analysed
-    TString typeData = "ON";
-    //TString typeData = "OFF";
-    //TString typeData = "MC";
-    cout << "typeData = " << typeData << endl;
-
-    TString typeMC   = "MC";
-    TString ext      = "3.root";
-    TString extout   = "4.root";
-
-    //------------------------------
-    // selection of g/h separation method
-    // and definition of final selections
-
-    TString XX("NN");
-    //TString XX("SC");
-    //TString XX("RF");
-    TString fhadronnessName("Had");
-    fhadronnessName += XX;
-    cout << "fhadronnessName = " << fhadronnessName << endl;
-
-    // maximum values of the hadronness and of |alpha|
-    Float_t maxhadronness   = 0.40;
-    Float_t maxalpha        = 20.0;
-    cout << "Maximum values of hadronness and |alpha| = "
-         << maxhadronness << ",  " << maxalpha << endl;
-
-
-    //------------------------------
-    // name of MC file to be used for calculating the eff. collection areas
-    TString filenameArea(outPath);
-    filenameArea += typeMC;
-    filenameArea += "_";
-    filenameArea += XX;
-    filenameArea += ext; 
-    cout << "filenameArea = " << filenameArea << endl;
-
-    //------------------------------
-    // name of file containing the eff. collection areas
-    TString collareaName(outPath);
-    collareaName += "area_";
-    collareaName += XX;
-    collareaName += ".root";
-    cout << "collareaName = " << collareaName << endl;
-
-    //------------------------------
-    // name of data file to be analysed
-    TString filenameData(outPath);
-    filenameData += typeData;
-    filenameData += "_";
-    filenameData += XX;
-    filenameData += ext;
-    cout << "filenameData = " << filenameData << endl;
-
-    //------------------------------
-    // name of output data file (after the final cuts)
-    TString filenameDataout(outPath);
-    filenameDataout += typeData;
-    filenameDataout += "_";
-    filenameDataout += XX;
-    filenameDataout += extout;
-    cout << "filenameDataout = " << filenameDataout << endl;
-
-
-    //====================================================================
-    cout << "-----------------------------------------------" << endl;
-    cout << "Start calculation of effective collection areas" << endl;
-    MParList  parlist;
-    MTaskList tasklist;
-
-    //---------------------------------------
-    // Setup the tasks to be executed
-    //
-    MReadMarsFile reader("Events", filenameArea);
-
-    MFCT1SelFinal cuthadrons;
-    cuthadrons.SetHadronnessName(fhadronnessName);
-    cuthadrons.SetCuts(maxhadronness, maxalpha);
-
-    MContinue conthadrons(&cuthadrons);
-
-    MHMcCT1CollectionArea* collarea = 
-              new MHMcCT1CollectionArea("MHMcCT1CollectionArea","",30,2.,5.);
-    MMcCT1CollectionAreaCalc effi;
-
-    //********************************
-    // entries in MParList
-
-    parlist.AddToList(&tasklist);
-    parlist.AddToList(collarea);
-
-    //********************************
-    // entries in MTaskList
-
-    tasklist.AddToList(&reader);   
-    tasklist.AddToList(&conthadrons);
-    tasklist.AddToList(&effi);
-
-    //********************************
-
-    //-----------------------------------------
-    // Execute event loop
-    //
-    MEvtLoop magic;
-    magic.SetParList(&parlist);
-
-    MProgressBar bar;
-    magic.SetProgressBar(&bar);
-    if (!magic.Eventloop())
-        return;
-
-    tasklist.PrintStatistics(0, kTRUE);
-
-    // Display the histograms
-    //
-    parlist.FindObject("MHMcCT1CollectionArea")->DrawClone("lego");
-
-    //---------------------------------------------
-    // Write histograms to a file 
-    //
-
-    TFile f(collareaName, "RECREATE");
-    collarea->GetHist()->Write();
-    collarea->GetHAll()->Write();
-    collarea->GetHSel()->Write();
-    f.Close();
-
-
-    cout << "Calculation of effective collection areas done" << endl;
-    cout << "-----------------------------------------------" << endl;    
-    //------------------------------------------------------------------
-
-
-    //*************************************************************************
-    //
-    // Analyse the data
-    //
-
-    MTaskList tliston;
-    MParList pliston;
-
-    TString hadName("Had");
-    hadName += XX;
-
-    TString fHilName    = "MHillas"; 
-    TString fHilSrcName = "MHillasSrc"; 
-    TString fnewparName = "MNewImagePar"; 
-
-    //-------------------------------------------
-    // create the tasks which should be executed 
-    //
-
-    MReadMarsFile read("Events", filenameData);
-    read.DisableAutoScheme();
-
-    //.......................................................................
-
-    if (WXXD)
-    {    
-      MWriteRootFile write(filenameDataout);
-
-      write.AddContainer("MRawRunHeader", "RunHeaders");
-      write.AddContainer("MTime",         "Events");
-      write.AddContainer("MMcEvt",        "Events");
-      write.AddContainer("MSrcPosCam",    "Events");
-      write.AddContainer("MSigmabar",     "Events");
-      write.AddContainer("MHillas",       "Events");
-      write.AddContainer("MHillasSrc",    "Events");
-      write.AddContainer("MNewImagePar",  "Events");
-
-      write.AddContainer("HadNN",         "Events");
-      write.AddContainer("HadSC",         "Events");
-
-      write.AddContainer("MEnergyEst",    "Events");
-    }
-
-    //-----------------------------------------------------------------
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-    MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
-    selfinalgh.SetCuts(maxhadronness, 100.0);
-    selfinalgh.SetHadronnessName(fhadronnessName);
-    MContinue contfinalgh(&selfinalgh);
-
-    MFillH fillhadnn("MHadNN[MHHadronness]", "HadNN");
-    MFillH fillhadsc("MHadSC[MHHadronness]", "HadSC");
-    //MFillH fillhadrf("MHadRF[MHHadronness]", "HadRF");
-
-    MFillH hfill1("MHHillas",    fHilName);
-    MFillH hfill2("MHStarMap",   fHilName);
-
-    MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
-
-    MFillH hfill3("MHHillasExt",   fHilSrcName);
-    MFillH hfill4("MHHillasSrc",   fHilSrcName);
-    MFillH hfill5("MHNewImagePar", fnewparName);
-
-    MFCT1SelFinal selfinal(fHilName, fHilSrcName);
-    selfinal.SetCuts(maxhadronness, maxalpha);
-    selfinal.SetHadronnessName(fhadronnessName);
-    MContinue contfinal(&selfinal);
-
-
-    //*****************************
-    // entries in MParList
-
-    pliston.AddToList(&tliston);
-    InitBinnings(&pliston);
-
-    pliston.AddToList(&hhilext); 
-
-    //*****************************
-    // entries in MTaskList
-    
-    tliston.AddToList(&read);
-
-    tliston.AddToList(&contfinalgh);
-
-    if (WXXD)
-      tliston.AddToList(&write);
-
-    tliston.AddToList(&contfinal);
-    tliston.AddToList(&fillhadnn);
-    tliston.AddToList(&fillhadsc);
-    //tliston.AddToList(&fillhadrf);
-    tliston.AddToList(&hfill1);
-    tliston.AddToList(&hfill2);
-    tliston.AddToList(&hfill3);
-    tliston.AddToList(&hfill4);
-    tliston.AddToList(&hfill5);
-
-
-
-    //*****************************
-
-    //-------------------------------------------
-    // Execute event loop
-    //
-    MProgressBar bar;
-    MEvtLoop evtloop;
-    evtloop.SetParList(&pliston);
-    evtloop.SetProgressBar(&bar);
-    //evtloop.Write();
-
-    Int_t maxevents = 1000000000;
-    //Int_t maxevents = 10000;
-    if ( !evtloop.Eventloop(maxevents) )
-        return;
-
-    tliston.PrintStatistics(0, kTRUE);
-    DeleteBinnings(&pliston);
-
-    //-------------------------------------------
-    // Display the histograms
-    //
-    TObject *c;
-
-    c = pliston.FindObject("MHadNN", "MHHadronness");
-    if (c)
-    {
-      c->DrawClone();
-      c->Print();
-    }
-    c = pliston.FindObject("MHadRF", "MHHadronness");
-    if (c)
-    {
-      c->DrawClone();
-      c->Print();
-    }
-    c = pliston.FindObject("MHadSC", "MHHadronness");
-    if (c)
-    {
-      c->DrawClone();
-      c->Print();
-    }
-
-    c = pliston.FindObject("MHHillas")->DrawClone();
-    c = pliston.FindObject("MHHillasExt")->DrawClone();
-    c = pliston.FindObject("MHHillasSrc")->DrawClone();
-    c = pliston.FindObject("MHNewImagePar")->DrawClone();
-    c = pliston.FindObject("MHStarMap")->DrawClone();
-
-
-    cout << "Macro CT1Analysis : End of Job D_XX" << endl;
-    cout << "=======================================================" << endl;
- }
-  //---------------------------------------------------------------------
-
-}
-
-
-
-
-
-
-
