
//#include "MagicEgyEst.C"


void InitBinnings(MParList *plist)
{
        gLog << "InitBinnings" << endl;

        //--------------------------------------------
        MBinning *binse = new MBinning("BinningE");
        //binse->SetEdgesLog(30, 1.0e2, 1.0e5);

	//This is Daniel's binning in energy:
        binse->SetEdgesLog(14, 296.296, 86497.6);
        plist->AddToList(binse);

        //--------------------------------------------

        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 *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,  50.0);
        plist->AddToList(binsb);

        MBinning *binssig = new MBinning("BinningSigma");
        binssig->SetEdges( 100,  0.0, 120.0);
        plist->AddToList(binssig);

        MBinning *binth = new MBinning("BinningTheta");
        // this is Daniel's binning in theta
        //Double_t yedge[8] = 
        //  {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99};
        // this is our binning
        //Double_t yedge[9] = 
        //               {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
        //TArrayD yed;
        //yed.Set(9,yedge);
        //binth->SetEdges(yed);
        
        binth->SetEdges(     1, 0.0, 90.0);
        //binth->SetEdgesCos( 10, 0.0, 90.0);
        plist->AddToList(binth);

        //MBinning *bincosth = new MBinning("BinningCosTheta");
        //Double_t yedge[9] = 
        //               {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
        //Double_t zedge[9]; 
        //for (Int_t i=0; i<9; i++)
	//{
	//zedge[8-i] = cos(yedge[i]/kRad2Deg);
	//}
        //TArrayD zed;
        //zed.Set(9,zedge);
        //bincosth->SetEdges(zed);
        //plist->AddToList(bincosth);

        MBinning *binsdiff = new MBinning("BinningDiffsigma2");
        binsdiff->SetEdges(100, -500.0, 1500.0);
        plist->AddToList(binsdiff);

        // robert ----------------------------------------------
        MBinning *binsalphaf = new MBinning("BinningAlphaFlux");
        binsalphaf->SetEdges(100, -100, 100);
        plist->AddToList(binsalphaf);

	MBinning *binsdifftime = new MBinning("BinningTimeDiff");
	binsdifftime->SetEdges(50, 0, 10);
	plist->AddToList(binsdifftime);

	MBinning *binstime = new MBinning("BinningTime");
	binstime->SetEdges(50, 44500, 61000);
	plist->AddToList(binstime);
}


void DeleteBinnings(MParList *plist)
{
        gLog << "DeleteBinnings" << endl;

        TObject *bin;

        //--------------------------------------------
        bin = plist->FindObject("BinningE");
        if (bin) delete 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("BinningM3Long");
        if (bin) delete bin;

        bin = plist->FindObject("BinningM3Trans");
        if (bin) delete bin;

        //.....
        bin = plist->FindObject("BinningSigmabar");
        if (bin) delete bin;

        bin = plist->FindObject("BinningSigma");
        if (bin) delete bin;

        bin = plist->FindObject("BinningTheta");
        if (bin) delete bin;

        bin = plist->FindObject("BinningCosTheta");
        if (bin) delete bin;

        bin = plist->FindObject("BinningDiffsigma2");
        if (bin) delete bin;


        // robert ----------------------------------------------
        bin = plist->FindObject("BinningAlphaFlux");
        if (bin) delete bin;

        bin = plist->FindObject("BinningTimeDiff");
        if (bin) delete bin;

        bin = plist->FindObject("BinningTime");
        if (bin) delete bin;
}



//************************************************************************
void ONOFFAnalysis()
{
      gLog.SetNoColors();

      if (gRandom)
        delete gRandom;
      gRandom = new TRandom3(0);

      //-----------------------------------------------
      TString tag = "080000";

      const char *offfile  = "*";

      const char *onfile  = "*";

      // Pratik
      //const char *mcfile = "MCGamma_calibrate";
      // Keiichi
      const char *mcfile = "calibrated_MCdata2";

      //-----------------------------------------------

      // path for input for Mars

     if (tag == "080000")  
     { 
       TString inPathON  = "/home/pcmagic14/wittek/CalibData/CrabSept2004/2004_09_21/";
       TString inPathOFF = "/home/pcmagic14/wittek/CalibData/OffSept2004/2004_09_18/";
       // Pratik
       //TString inPathMC  = "/home/pcmagic21/pratik/mcdata/gamma/MCcalibrate/";
       // Keiichi
       TString inPathMC  = "/.magic/data21a/mase/Mars/Mars041103/DataCalibUV/";

     }

      // path for output from Mars
     TString outPath = "/.magic/data21a/wittek/";
     outPath += tag;
     outPath += "/";

      //-----------------------------------------------

      //TEnv env("macros/CT1env.rc");
      //Bool_t printEnv = kFALSE;

    //************************************************************************

    // Job A : 
    //  - produce MHSigmaTheta plots for ON, OFF and MC data
    //  - write out (or read in) these MHSigmaTheta plots
    //  - read ON (or OFF or MC) data
    //  - pad the events; 
    //  - write root file for ON (or OFF or MC) data (ON1.root, ...);

    Bool_t JobA    = kTRUE;  
    Bool_t GPadON  = kFALSE;    // \  generate Pad histograms 
    Bool_t GPadOFF = kFALSE;    //  | and write them onto disk
    Bool_t GPadMC  = kFALSE;    // /
    Bool_t Merge   = kFALSE;   // read the Pad histograms, merge them
                               // and write them onto disk
    Bool_t Wout    = kTRUE;   // \  read in merged padding histograms and
                               //  | write out root file of padded data
                               // /  (ON1.root, OFF1.root or MC1.root) 
    //TString typeInput("ON");
    TString typeInput("OFF");
    //TString typeInput("MC");


    // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file 
    //  - if CTrainRF = TRUE : create matrices of training events
    //                         and root files of training and test events
    //  - if RTrainRF = TRUE : read in training matrices for hadrons and gammas
    //  - if RTree    = TRUE : read in trees, otherwise create trees
    //  - calculate hadroness for method of RANDOM FOREST
    //  - update the input files with the hadronesses (ON2.root, OFF2.root
    //     or MC2.root)

    Bool_t JobB_RF_UP  = kFALSE;  
    Bool_t CTrainRF    = kFALSE;  // create  matrices of training events
                                  // and root files of training and test events
    Bool_t RTrainRF    = kFALSE;  // read in matrices of training events
    Bool_t RTree       = kFALSE;  // read in trees (otherwise grow them)
    Bool_t WRF         = kFALSE;  // update input root file ?


    // Job B_SC_UP : read ON2.root (or MC2.root) file 
    //  - depending on WParSC : create (or read in) supercuts parameter values
    //  - calculate hadroness for the SUPERCUTS
    //  - update the input files with the hadroness (==>ON3.root or MC3.root)

    Bool_t JobB_SC_UP  = kFALSE;
    Bool_t CMatrix     = kFALSE;  // create training and test matrices 
    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file
    Bool_t WOptimize   = kFALSE;  // do optimization using the training sample
                                  // and write supercuts parameter values 
                                  // onto the file parSCfile
    Bool_t RTest       = kFALSE;  // test the supercuts using the test matrix
    Bool_t WSC         = kFALSE;  // update input root file ?


    // Job C: 
    //  - read ON3.root and MC3.root files
    //    which should have been updated to contain the hadronnesses  
    //    for the method of 
    //              RF
    //              SUPERCUTS and
    //  - produce Neyman-Pearson plots

    Bool_t JobC  = kFALSE;  


    // Job D :  
    //  - select g/h separation method XX
    //  - read ON3 (or MC3) root file
    //  - apply cuts in hadronness
    //  - make plots

    Bool_t JobD  = kFALSE;  



    // Job E_XX : extended version of E_XX (including flux plots)  
    //  - select g/h separation method XX
    //  - read MC root file
    //  - calculate eff. collection area
    //  - optimize energy estimator
    //  - read ON root file 
    //  - apply final cuts
    //  - calculate flux
    //  - write root file for ON data after final cuts 


    Bool_t JobE_XX  = kFALSE;  
    Bool_t CCollArea= kFALSE;  // calculate eff. collection areas
    Bool_t OEEst    = kFALSE;  // optimize energy estimator
    Bool_t WEX      = kFALSE;  // update root file  ?
    Bool_t WRobert  = kFALSE;  // write out Robert's file  ?



    //************************************************************************

    
  //---------------------------------------------------------------------
  // Job A
  //=========

    //  - produce the histograms "sigmabar versus Theta", etc. 
    //    for ON, OFF and MC data (to be used for the padding)
    //
    //  - write root file of padded ON (OFF, MC) events (ON1.root, ...) 
    //    (after the standard cuts, before the g/h separation)


 if (JobA)
 {
    gLog << "=====================================================" << endl;
    gLog << "Macro ONOFFAnalysis : Start of Job A" << endl;
    gLog << "" << endl;
    gLog << "Macro ONOFFAnalysis : JobA, GPadON, GPadOFF, GPadMC, Merge, Wout = " 
         << (JobA    ? "kTRUE" : "kFALSE")  << ",  " 
         << (GPadON  ? "kTRUE" : "kFALSE")  << ",  " 
         << (GPadOFF ? "kTRUE" : "kFALSE")  << ",  " 
         << (GPadMC  ? "kTRUE" : "kFALSE")  << ",  " 
         << (Merge   ? "kTRUE" : "kFALSE")  << ",  " 
         << (Wout    ? "kTRUE" : "kFALSE")  << endl;
    
    //--------------------------------------------------

    TString fNamePedPhotCam("MPedPhotCamFromData");
    // for Keiichi's file
    //TString fNamePedPhotCam("MPedPhotCam");


    //************************************************************
    // generate histograms to be used in the padding
    // 
    // read ON, OFF and MC data files
    // generate (or read in) the padding histograms for ON, OFF and MC data
    //

    MPad pad;
    pad.SetName("MPad");

    //--------------------------------------------------
    // names of ON and OFF files to be read
    // for generating the histograms to be used in the padding 
    TString fileON  = inPathON;
    fileON += onfile;
    fileON += ".root";

    TString fileOFF = inPathOFF;
    fileOFF += offfile;
    fileOFF += ".root";

    TString fileMC = inPathMC;
    fileMC += mcfile;
    fileMC += ".root";

    //--------------------------------------------------
    // name of files to contain the paddding histograms of ON, OFF and MC data
      TString NamePadON(outPath);
      NamePadON += "PadON";
      NamePadON += ".root";

      TString NamePadOFF(outPath);
      NamePadOFF += "PadOFF";
      NamePadOFF += ".root";

      TString NamePadMC(outPath);
      NamePadMC += "PadMC";
      NamePadMC += ".root";

    // name of file to conatin the merged histograms for the padding
    TString outNameSigTh = outPath;
    outNameSigTh += "SigmaTheta";
    outNameSigTh += ".root";

    //--------------------------------------------------

    if (GPadON || GPadOFF || GPadMC)
    {
      // generate the padding histograms
      gLog << "=====================================================" << endl;
      gLog << "Start generating the padding histograms" << endl;


    gLog << "fileON, fileOFF, fileMC = " << fileON << ",  " 
         << fileOFF << ",  " << fileMC   << endl;



    //--------------------------------------------------
      MMakePadHistograms makepad;
      makepad.SetMaxEvents(10000);
      makepad.SetNamePedPhotCam(fNamePedPhotCam);
      makepad.SetPedestalLevel(2.0);
      makepad.SetUseInterpolation(kTRUE);
      makepad.SetProcessPedestal(kTRUE);
      makepad.SetProcessTime(kFALSE);

      //-----------------------------------------
      // ON events

      if (GPadON)
      {
        makepad.SetDataType("ON");
        makepad.SetNameInputFile(fileON);
        makepad.SetNameOutputFile(NamePadON);
        makepad.MakeHistograms();
      }

      //-----------------------------------------
      // OFF events

      if (GPadOFF)
      {
        makepad.SetDataType("OFF");
        makepad.SetNameInputFile(fileOFF);
        makepad.SetNameOutputFile(NamePadOFF);
        makepad.MakeHistograms();
      }

      //-----------------------------------------
      // MC events

      if (GPadMC)
      {
        makepad.SetDataType("MC");
        makepad.SetNameInputFile(fileMC);
        makepad.SetNameOutputFile(NamePadMC);
        makepad.MakeHistograms();
      }

      //-----------------------------------------


      gLog << "" << endl;
      gLog << "End of generating the padding histograms" << endl;
      gLog << "=====================================================" << endl;
    }

    //************************************************************

    if (Merge)
    {
      gLog << "=====================================================" << endl;
      gLog << "Start of merging the padding histograms" << endl;
      gLog << "" << endl;

      pad.MergeONOFFMC(NamePadON, NamePadOFF, NamePadMC, outNameSigTh);
      //pad.MergeONOFFMC(NamePadON, "", NamePadMC, outNameSigTh);
      //pad.MergeONOFFMC(NamePadON, NamePadOFF, "", outNameSigTh);
      //pad.MergeONOFFMC("", NamePadOFF, NamePadMC, outNameSigTh);

      gLog << "" << endl;
      gLog << "End of merging the padding histograms" << endl;
      gLog << "=====================================================" << endl;
    }
    // end of Merge



    //************************************************************

  if (Wout)
  {
    // read the target padding histograms ---------------------------
    pad.ReadPaddingDist(outNameSigTh);


    gLog << "=====================================================" << endl;
    gLog << "Start the padding" << endl;

    //--------------------------------------------------
    // type of data to be padded 
    gLog << "typeInput = " << typeInput << endl;

    //-------------------------------------------
    // name of input root file
    if (typeInput == "ON")
      TString filenamein(fileON);
    else if (typeInput == "OFF")
      TString filenamein(fileOFF);
    else if (typeInput == "MC")
      TString filenamein(fileMC);
    gLog << "data to be padded : " << filenamein << endl;

    // name of output root file
    TString outNameImage = outPath;
    outNameImage += "Hillas";
    outNameImage += typeInput;
    outNameImage += "1.root";
    gLog << "padded data to be written onto : " << outNameImage << endl;

    //-----------------------------------------------------------
    pad.SetDataType(typeInput);
    pad.SetNamePedPhotCam(fNamePedPhotCam);

    MTaskList tliston;
    MParList pliston;

    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");

    //-------------------------------------------
    // create the tasks which should be executed 
    //

    MReadMarsFile read("Events", filenamein);
    read.DisableAutoScheme();

    MGeomApply        apply;

    // a way to find out whether one is dealing with MC :
    //MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5);  // MC
    //f1.SetName("Select MC");
    //MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5);  // data
    //f2.SetName("Select Data");



      MBadPixelsCalc badcalc;
      badcalc.SetNamePedPhotContainer(fNamePedPhotCam);
      badcalc.SetPedestalLevel(2.0);
      badcalc.SetName("BadCalc");

      MBadPixelsTreat badtreat;
      badtreat.SetNamePedPhotCam(fNamePedPhotCam);
      badtreat.SetUseInterpolation(kTRUE);
      badtreat.SetProcessPedestal(kTRUE);
      badtreat.SetProcessTimes(kFALSE);
      badtreat.SetName("BadTreat");

      MFSelBasic selbasic;
      MContinue contbasic(&selbasic);
      contbasic.SetName("SelBasic");

      MHBadPixels bad("BadPixels");
      bad.SetNamePedPhotCam(fNamePedPhotCam);
      MFillH fillbad("BadPixels[MHBadPixels]", "MBadPixelsCam");
      fillbad.SetName("FillBad");

      MHSigmaTheta sigth("SigmaTheta");
      sigth.SetNamePedPhotCam(fNamePedPhotCam);
      MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "");
      fillsigtheta.SetName("FillSigTheta");    


    MImgCleanStd    clean(3.0, 2.5);
    //clean.SetMethod(MImgCleanStd::kDemocratic);
    clean.SetCleanRings(3); 
    clean.SetName("Clean");
    clean.SetNamePedPhotCam(fNamePedPhotCam);

    // calculation of  image parameters ---------------------
    TString fHilName    = "MHillas";
    TString fHilNameExt = "MHillasExt";
    TString fHilNameSrc = "MHillasSrc";
    TString fImgParName = "MNewImagePar";

    MHillasCalc    hcalc;
    hcalc.SetNameHillas(fHilName);
    hcalc.SetNameHillasExt(fHilNameExt);
    hcalc.SetNameNewImagePar(fImgParName);

    MHillasSrcCalc hsrccalc("MSrcPosCam", fHilNameSrc);
    hsrccalc.SetInput(fHilName);

    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");
    // --------------------------------------------------

    MFSelStandard selstandard(fHilNameSrc);
    selstandard.SetHillasName(fHilName);
    selstandard.SetImgParName(fImgParName);
    selstandard.SetCuts(200, 5, 200, 0.0, 5.0, 0.1, 0.07);
    MContinue contstandard(&selstandard);
    contstandard.SetName("SelStandard");

    /*  
      MWriteRootFile write(outNameImage);

      write.AddContainer("MRawRunHeader", "RunHeaders");
      //write.AddContainer("MMcRunHeader",  "RunHeaders", kFALSE);
      //write.AddContainer("MTime",         "Events");
      write.AddContainer("MPointingPos", "Events");
      write.AddContainer("MSrcPosCam",    "Events");
      write.AddContainer("MHillas",       "Events");
      write.AddContainer("MHillasExt",    "Events");
      write.AddContainer("MHillasSrc",    "Events");
      write.AddContainer("MNewImagePar",  "Events");
    */      

    //*****************************
    // entries in MParList
    
    pliston.AddToList(&tliston);
    InitBinnings(&pliston);
    pliston.AddToList(&bad);
    pliston.AddToList(&sigth);

    //*****************************
    // entries in MTaskList
    
    tliston.AddToList(&read);
    //tliston.AddToList(&f1);
    //tliston.AddToList(&f2);
    tliston.AddToList(&apply);

    tliston.AddToList(&badcalc);
    tliston.AddToList(&badtreat);

    tliston.AddToList(&contbasic);
    tliston.AddToList(&pad); 

    
    tliston.AddToList(&fillbad);
    tliston.AddToList(&fillsigtheta);
    tliston.AddToList(&clean);

    tliston.AddToList(&hcalc);
    tliston.AddToList(&hsrccalc);

    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    tliston.AddToList(&contstandard);
    //tliston.AddToList(&write);
    

    //*****************************

    //-------------------------------------------
    // Execute event loop
    //
    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&pliston);
    //evtloop.ReadEnv(env, "", printEnv);
    evtloop.SetProgressBar(&bar);
    //  evtloop.Write();

    //Int_t maxevents = -1;
    Int_t maxevents = 2000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // Display the histograms

    pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone();
    pliston.FindObject("BadPixels",  "MHBadPixels")->DrawClone();

    pliston.FindObject("MHHillas")->DrawClone();
    pliston.FindObject("MHHillasExt")->DrawClone();
    pliston.FindObject("MHHillasSrc")->DrawClone();
    pliston.FindObject("MHNewImagePar")->DrawClone();
    pliston.FindObject("MHStarMap")->DrawClone();

    //DeleteBinnings(&pliston);

      gLog << "End of padding" << endl;
      gLog << "=====================================================" << endl;
  }  


    gLog << "Macro ONOFFAnalysis : End of Job A" << endl;
    gLog << "===================================================" << endl;
 }





  //---------------------------------------------------------------------
  // Job B_RF_UP
  //============


    //  - create (or read in) the matrices of training 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 ONOFFAnalysis : Start of Job B_RF_UP" << endl;

    gLog << "" << endl;
    gLog << "Macro ONOFFAnalysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = "
         << (JobB_RF_UP ? "kTRUE" : "kFALSE")  << ",  " 
         << (RTrainRF ?   "kTRUE" : "kFALSE")  << ",  " 
         << (CTrainRF ?   "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;


    TString hadRFName = "HadRF";
    Float_t maxhadronness =  0.23;
    Float_t maxalpha      =  20.0;
    Float_t maxdist       =  10.0;

    TString fHilName    = "MHillas"; 
    TString fHilNameExt = "MHillasExt"; 
    TString fHilNameSrc = "MHillasSrc"; 
    TString fImgParName = "MNewImagePar"; 


    TString extin  = "1.root";
    TString extout = "2.root";

    //--------------------------------------------
    // for the analysis using ON data only set typeMatrixHadrons = "ON"
    //                        ON and OFF data                    = "OFF"
    TString typeMatrixHadrons = "OFF";
    gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl;


    // file to be updated (ON, OFF or MC)

    //TString typeInput = "ON";
    TString typeInput = "OFF";
    //TString typeInput = "MC";
    gLog << "typeInput = " << typeInput << endl;

    // name of input root file
    TString NameData = outPath;
    NameData += typeInput;
    TString inNameData(NameData);
    inNameData += extin;
    gLog << "inNameData = " << inNameData << endl; 

    // name of output root file
    TString outNameData(NameData);
    outNameData += extout;
    gLog << "outNameData = " << outNameData << endl; 

    //--------------------------------------------
    // files to be read for generating 
    //    - the matrices of training events
    //    - and the root files of training and test events


    // "hadrons" :
    TString filenameHad = outPath;
    filenameHad += typeMatrixHadrons;
    filenameHad += extin;
    Int_t howManyHadronsTrain = 12000;
    Int_t howManyHadronsTest  = 12000;
    gLog << "filenameHad = "    << filenameHad << ",   howManyHadronsTrain = "
         << howManyHadronsTrain << ",   howManyHadronsTest = "
         << howManyHadronsTest  << endl; 
    

    // "gammas" :
    TString filenameMC = outPath;
    filenameMC += "MC";
    filenameMC += extin;
    Int_t howManyGammasTrain = 12000;
    Int_t howManyGammasTest  = 12000;
    gLog << "filenameMC = "    << filenameMC << ",   howManyGammasTrain = "
         << howManyGammasTrain << ",   howManyGammasTest = "
         << howManyGammasTest  << endl; 
    
    //--------------------------------------------
    // files for the matrices of training events 

    TString NameGammas = outPath;
    NameGammas += "RFmatrix_gammas_Train_";
    NameGammas += "MC";
    NameGammas += extin;

    TString NameHadrons = outPath;
    NameHadrons += "RFmatrix_hadrons_Train_";
    NameHadrons += typeMatrixHadrons;
    NameHadrons += extin;


    //--------------------------------------------
    // root files for the training events 

    TString NameGammasTrain = outPath;
    NameGammasTrain += "RF_gammas_Train_";
    NameGammasTrain += "MC";
    TString inNameGammasTrain(NameGammasTrain);    
    inNameGammasTrain += extin;
    TString outNameGammasTrain(NameGammasTrain);    
    outNameGammasTrain += extout;


    TString NameHadronsTrain = outPath;
    NameHadronsTrain += "RF_hadrons_Train_";
    NameHadronsTrain += typeMatrixHadrons;
    TString inNameHadronsTrain(NameHadronsTrain);    
    inNameHadronsTrain += extin;
    TString outNameHadronsTrain(NameHadronsTrain);    
    outNameHadronsTrain += extout;


    //--------------------------------------------
    // root files for the test events 

    TString NameGammasTest = outPath;
    NameGammasTest += "RF_gammas_Test_";
    NameGammasTest += "MC";
    TString inNameGammasTest(NameGammasTest);    
    inNameGammasTest += extin;
    TString outNameGammasTest(NameGammasTest);    
    outNameGammasTest += extout;

    TString NameHadronsTest = outPath;
    NameHadronsTest += "RF_hadrons_Test_";
    NameHadronsTest += typeMatrixHadrons;
    TString inNameHadronsTest(NameHadronsTest);    
    inNameHadronsTest += extin;
    TString outNameHadronsTest(NameHadronsTest);    
    outNameHadronsTest += extout;

    //--------------------------------------------------------------------


    MHMatrix matrixg("MatrixGammas");
    matrixg.EnableGraphicalOutput();

    matrixg.AddColumn("cos((MPointingPos.fZd)/kRad2Deg)");
    matrixg.AddColumn("MSigmabar.fSigmabar");
    matrixg.AddColumn("log10(MHillas.fSize)");
    matrixg.AddColumn("MHillasSrc.fDist");
    matrixg.AddColumn("MHillas.fWidth");
    matrixg.AddColumn("MHillas.fLength");
    matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
    matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
    matrixg.AddColumn("MNewImagePar.fConc");
    matrixg.AddColumn("MNewImagePar.fLeakage1");

    MHMatrix matrixh("MatrixHadrons");
    matrixh.EnableGraphicalOutput();

    matrixh.AddColumns(matrixg.GetColumns());

    //--------------------------------------------
    // file of trees of the random forest 

    TString outRF = outPath;
    outRF += "RF.root";


   //*************************************************************************
   // read in matrices of training events
if (RTrainRF)
  {
    const char* mtxName = "MatrixGammas";

    gLog << "" << endl;
    gLog << "========================================================" << endl;
    gLog << "Get matrix for (gammas)" << endl;
    gLog << "matrix name        = " << mtxName << endl;
    gLog << "name of root file  = " << NameGammas << endl;
    gLog << "" << endl;


    // read in the object with the name 'mtxName' from file 'NameGammas'
    //
    TFile fileg(NameGammas); 

    matrixg.Read(mtxName);
    matrixg.Print("SizeCols");


    //***************************************************************** 

    const char* mtxName = "MatrixHadrons";

    gLog << "" << endl;
    gLog << "========================================================" << endl;
    gLog << " Get matrix for (hadrons)" << endl;
    gLog << "matrix name        = " << mtxName << endl;
    gLog << "name of root file  = " << NameHadrons << endl;
    gLog << "" << endl;


    // read in the object with the name 'mtxName' from file 'NameHadrons'
    //
    TFile fileh(NameHadrons); 

    matrixh.Read(mtxName);
    matrixh.Print("SizeCols");
  }


   //*************************************************************************
   // create matrices of training events
   // and root files of training and test events
 
if (CTrainRF)
  {
    gLog << "" << endl;
    gLog << "========================================================" << endl;
    gLog << " Create matrices of training events and root files of training and test events" 
         << endl;
    gLog << " Gammas :" << endl;
    gLog << "---------" << endl;

    MParList  plistg;
    MTaskList tlistg;

    MReadMarsFile  readg("Events", filenameMC);
    readg.DisableAutoScheme();

    TString mgname("costhg");
    MBinning bing("Binning"+mgname);
    bing.SetEdges(10, 0., 1.0);

    MH3 gref("cos((MPointingPos.fZd)/kRad2Deg)");
    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(howManyGammasTrain+howManyGammasTest);
    selectorg.SetName("selectGammasTrainTest");
    selectorg.SetInverted();
    //selectorg.SetUseOrigDistribution(kTRUE);

    MContinue contg(&selectorg);
    contg.SetName("ContGammas");

    Double_t probg = ( (Double_t) howManyGammasTrain )
                   / ( (Double_t)(howManyGammasTrain+howManyGammasTest) );
    MFRandomSplit splitg(probg);

    MFillH fillmatg("MatrixGammas");
    fillmatg.SetFilter(&splitg);
    fillmatg.SetName("fillGammas");

    //-----------------------
    // for writing the root files of training and test events
    // for gammas
    
    MWriteRootFile writetraing(inNameGammasTrain, "RECREATE");
    writetraing.SetName("WriteGammasTrain");
    writetraing.SetFilter(&splitg);

      writetraing.AddContainer("MRawRunHeader", "RunHeaders");
      writetraing.AddContainer("MTime",         "Events");
      writetraing.AddContainer("MPointingPos",  "Events");
      writetraing.AddContainer("ThetaOrig",     "Events");
      writetraing.AddContainer("MSrcPosCam",    "Events");
      writetraing.AddContainer("MSigmabar",     "Events");
      writetraing.AddContainer("MHillas",       "Events");
      writetraing.AddContainer("MHillasExt",    "Events");
      writetraing.AddContainer("MHillasSrc",    "Events");
      writetraing.AddContainer("MNewImagePar",  "Events");

    MContinue contgtrain(&splitg);
    contgtrain.SetName("ContGammaTrain");

    MWriteRootFile writetestg(inNameGammasTest, "RECREATE");
    writetestg.SetName("WriteGammasTest");

      writetestg.AddContainer("MRawRunHeader", "RunHeaders");
      writetestg.AddContainer("MTime",         "Events");
      writetestg.AddContainer("MPointingPos",  "Events");
      writetestg.AddContainer("ThetaOrig",     "Events");
      writetestg.AddContainer("MSrcPosCam",    "Events");
      writetestg.AddContainer("MSigmabar",     "Events");
      writetestg.AddContainer("MHillas",       "Events");
      writetestg.AddContainer("MHillasExt",    "Events");
      writetestg.AddContainer("MHillasSrc",    "Events");
      writetestg.AddContainer("MNewImagePar",  "Events");

    //-----------------------
    
    //*****************************   fill gammas   ***  
    // entries in MParList
    
    plistg.AddToList(&tlistg);
    InitBinnings(&plistg);

    plistg.AddToList(&matrixg);

    //*****************************
    // entries in MTaskList
    
    tlistg.AddToList(&readg);
    tlistg.AddToList(&contg);

    tlistg.AddToList(&splitg);
    tlistg.AddToList(&fillmatg);
    tlistg.AddToList(&writetraing);
    tlistg.AddToList(&contgtrain);

    tlistg.AddToList(&writetestg);

    //*****************************

    MProgressBar matrixbar;
    MEvtLoop evtloopg;
    evtloopg.SetName("FillGammaMatrix");
    evtloopg.SetParList(&plistg);
    //evtloopg.ReadEnv(env, "", printEnv);
    evtloopg.SetProgressBar(&matrixbar);

    Int_t maxevents = -1;
    if (!evtloopg.Eventloop(maxevents))
        return;

    tlistg.PrintStatistics(0, kTRUE);

    matrixg.Print("SizeCols");
    Int_t generatedgTrain = matrixg.GetM().GetNrows();
    if ( fabs(generatedgTrain-howManyGammasTrain) > 
                                    3.0*sqrt(howManyGammasTrain) )
    {
      gLog << "ONOFFAnalysis.C : no.of generated gamma training events (" 
           << generatedgTrain << ") is incompatible with the no.of requested events (" 
           << howManyGammasTrain << ")" << endl;   
    }


    Int_t generatedgTest = writetestg.GetNumExecutions();
    if ( fabs(generatedgTest-howManyGammasTest) > 
                                    3.0*sqrt(howManyGammasTest) )
    {
      gLog << "ONOFFAnalysis.C : no.of generated gamma test events (" 
           << generatedgTest << ") is incompatible with the no.of requested events (" 
           << howManyGammasTest << ")" << endl;   
    }

    //*****************************   fill hadrons   ***  
    gLog << "---------------------------------------------------------------"
         << endl;
    gLog << " Hadrons :" << endl;
    gLog << "----------" << endl;

    MParList  plisth;
    MTaskList tlisth;

    MReadMarsFile  readh("Events", filenameHad);
    readh.DisableAutoScheme();

    TString mhname("costhh");
    MBinning binh("Binning"+mhname);
    binh.SetEdges(10, 0., 1.0);

    //MH3 href("cos((MPointingPos.fZd)/kRad2Deg)");
    //href.SetName(mhname);
    //MH::SetBinning(&href.GetHist(), &binh);
    //for (Int_t i=1; i<=href.GetNbins(); i++)
    //  href.GetHist().SetBinContent(i, 1.0);

    //use the original distribution from the gammas
    MH3 &href = *(selectorg.GetHistOrig());

    MFEventSelector2 selectorh(href);
    selectorh.SetNumMax(howManyHadronsTrain+howManyHadronsTest);
    selectorh.SetName("selectHadronsTrainTest");
    selectorh.SetInverted();

    MContinue conth(&selectorh);
    conth.SetName("ContHadrons");

    Double_t probh = ( (Double_t) howManyHadronsTrain )
                   / ( (Double_t)(howManyHadronsTrain+howManyHadronsTest) );
    MFRandomSplit splith(probh);

    MFillH fillmath("MatrixHadrons");
    fillmath.SetFilter(&splith);
    fillmath.SetName("fillHadrons");

    //-----------------------
    // for writing the root files of training and test events
    // for hadrons
    
    MWriteRootFile writetrainh(inNameHadronsTrain, "RECREATE");
    writetrainh.SetName("WriteHadronsTrain");
    writetrainh.SetFilter(&splith);

      writetrainh.AddContainer("MRawRunHeader", "RunHeaders");
      writetrainh.AddContainer("MTime",         "Events");
      writetrainh.AddContainer("MPointingPos",  "Events");
      writetrainh.AddContainer("ThetaOrig",     "Events");
      writetrainh.AddContainer("MSrcPosCam",    "Events");
      writetrainh.AddContainer("MSigmabar",     "Events");
      writetrainh.AddContainer("MHillas",       "Events");
      writetrainh.AddContainer("MHillasExt",    "Events");
      writetrainh.AddContainer("MHillasSrc",    "Events");
      writetrainh.AddContainer("MNewImagePar",  "Events");

    MContinue conthtrain(&splith);

    MWriteRootFile writetesth(inNameHadronsTest, "RECREATE");
    writetesth.SetName("WriteHadronsTest");

      writetesth.AddContainer("MRawRunHeader", "RunHeaders");
      writetesth.AddContainer("MTime",         "Events");
      writetesth.AddContainer("MPointingPos",  "Events");
      writetesth.AddContainer("ThetaOrig",     "Events");
      writetesth.AddContainer("MSrcPosCam",    "Events");
      writetesth.AddContainer("MSigmabar",     "Events");
      writetesth.AddContainer("MHillas",       "Events");
      writetesth.AddContainer("MHillasExt",    "Events");
      writetesth.AddContainer("MHillasSrc",    "Events");
      writetesth.AddContainer("MNewImagePar",  "Events");


    //*****************************  
    // entries in MParList
    
    plisth.AddToList(&tlisth);
    InitBinnings(&plisth);

    plisth.AddToList(&matrixh);

    //*****************************
    // entries in MTaskList
    
    tlisth.AddToList(&readh);
    tlisth.AddToList(&conth);

    tlisth.AddToList(&splith);
    tlisth.AddToList(&fillmath);
    tlisth.AddToList(&writetrainh);
    tlisth.AddToList(&conthtrain);

    tlisth.AddToList(&writetesth);

    //*****************************

    MProgressBar matrixbar;
    MEvtLoop evtlooph;
    evtlooph.SetName("FillHadronMatrix");
    evtlooph.SetParList(&plisth);
    //evtlooph.ReadEnv(env, "", printEnv);
    evtlooph.SetProgressBar(&matrixbar);

    Int_t maxevents = -1;
    if (!evtlooph.Eventloop(maxevents))
        return;

    tlisth.PrintStatistics(0, kTRUE);

    matrixh.Print("SizeCols");
    Int_t generatedhTrain = matrixh.GetM().GetNrows();
    if ( fabs(generatedhTrain-howManyHadronsTrain) > 
                                    3.0*sqrt(howManyHadronsTrain) )
    {
      gLog << "ONOFFAnalysis.C : no.of generated hadron training events (" 
           << generatedhTrain << ") is incompatible with the no.of requested events (" 
           << howManyHadronsTrain << ")" << endl;   
    }


    Int_t generatedhTest = writetesth.GetNumExecutions();
    if ( fabs(generatedhTest-howManyHadronsTest) > 
                                    3.0*sqrt(howManyHadronsTest) )
    {
      gLog << "ONOFFAnalysis.C : no.of generated gamma test events (" 
           << generatedhTest << ") is incompatible with the no.of requested events (" 
           << howManyHadronsTest << ")" << endl;   
    }


    //*****************************************************  


    // write out matrices of training events 

    gLog << "" << endl;
    gLog << "========================================================" << endl;
    gLog << "Write out matrices of training events" << endl;


      //-------------------------------------------
      // "gammas"
      gLog << "Gammas :" << endl;    
      matrixg.Print("SizeCols");

      TFile writeg(NameGammas, "RECREATE", "");
      matrixg.Write();

      gLog << "" << endl;
      gLog << "Macro ONOFFAnalysis : matrix of training events for gammas written onto file "
           << NameGammas << endl;

      //-------------------------------------------
      // "hadrons"
      gLog << "Hadrons :" << endl;    
      matrixh.Print("SizeCols");

      TFile writeh(NameHadrons, "RECREATE", "");
      matrixh.Write();

      gLog << "" << endl;
      gLog << "Macro ONOFFAnalysis : matrix of training events for hadrons written onto file "
           << NameHadrons << endl;

  }
   //**********   end of creating matrices of training 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.SetName("ReadRFTrees");
      evtlooptr.SetParList(&plisttr);
      if (!evtlooptr.Eventloop())
        return;

      tlisttr.PrintStatistics(0, kTRUE);

      gLog << "ONOFFAnalysis : RF trees were read in from file "
           << outRF << endl;

    // get adresses of objects which are used in the next eventloop
    fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
    if (!fRanForest)
    {
        gLog << err << dbginf << "MRanForest not found... aborting." << endl;
        return kFALSE;
    }

    fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
    if (!fRanTree)                                  
    {                                                                          
        gLog << 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 ONOFFAnalysis : 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");

    MFillH fillh2("MHRanForestGini");

    // list of tasks for the loop over the trees
    
    tlist2.AddToList(&rfgrow2);
    tlist2.AddToList(&rfwrite2);
    tlist2.AddToList(&fillh2);

    //-------------------
    // Execute tree loop
    //
    MEvtLoop treeloop;
    treeloop.SetName("GrowRFTrees");
    treeloop.SetParList(&plist2);

    if ( !treeloop.Eventloop() )
        return;

    tlist2.PrintStatistics(0, kTRUE);

    plist2.FindObject("MHRanForestGini")->DrawClone();


    // get adresses of objects which are used in the next eventloop
    fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
    if (!fRanForest)
    {
        gLog << err << dbginf << "MRanForest not found... aborting." << endl;
        return kFALSE;
    }

    fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
    if (!fRanTree)                                  
    {                                                                          
        gLog << err << dbginf << "MRanTree not found... aborting." << endl;    
        return kFALSE;
    }

    }
    // end of growing the trees of the random forest
    //-----------------------------------------------------------------


    //-----------------------------------------------------------------
    // Update the root files with the RF hadronness
    //

 if (WRF)
  {
    //TString fileName(inNameHadronsTrain);
    //TString outName(outNameHadronsTrain);

    //TString fileName(inNameHadronsTest);
    //TString outName(outNameHadronsTest);

    //TString fileName(inNameGammasTrain);
    //TString outName(outNameGammasTrain);

    //TString fileName(inNameGammasTest);
    //TString outName(outNameGammasTest);

    TString fileName(inNameData);
    TString outName(outNameData);



    gLog << "" << endl;
    gLog << "========================================================" << endl;
    gLog << "Update root file '" <<  fileName 
         << "' with the RF hadronness; ==> " << outName << endl;

   
    MTaskList tliston;
    MParList pliston;


    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");

    //-------------------------------------------
    // create the tasks which should be executed 
    //

    MReadMarsFile read("Events", fileName);
    read.DisableAutoScheme();


    //.......................................................................
    // calculate hadronnes for method of RANDOM FOREST


    MRanForestCalc rfcalc;
    rfcalc.SetHadronnessName(hadRFName);


    //.......................................................................

      //MWriteRootFile write(outName, "UPDATE");
      MWriteRootFile write(outName, "RECREATE");

      write.AddContainer("MRawRunHeader", "RunHeaders");
      write.AddContainer("MTime",         "Events");
      write.AddContainer("MPointingPos",  "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(hadRFName,       "Events");

    //-----------------------------------------------------------------


    MFSelFinal 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");

    MFSelFinal selfinal(fHilNameSrc);
    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    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);
    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);

    pliston.AddToList(&binsalphaabs);
    pliston.AddToList(&alphaabs);


    //*****************************
    // entries in MTaskList
    
    tliston.AddToList(&read);

    tliston.AddToList(&rfcalc);
    tliston.AddToList(&fillranfor);
    tliston.AddToList(&fillhadrf);

    tliston.AddToList(&write);
    tliston.AddToList(&contfinalgh);

    tliston.AddToList(&alpha);
    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    tliston.AddToList(&contfinal);

    //*****************************

    //-------------------------------------------
    // Execute event loop
    //
    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetName("UpdateRootFile");
    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("hadRF", "MHHadronness")->Print();

    pliston.FindObject("MHHillas")->DrawClone();
    pliston.FindObject("MHHillasExt")->DrawClone();
    pliston.FindObject("MHHillasSrc")->DrawClone();
    pliston.FindObject("MHNewImagePar")->DrawClone();
    pliston.FindObject("MHStarMap")->DrawClone();


     //-------------------------------------------
    // fit alpha distribution to get the number of excess events and
    // calculate significance of gamma signal in the alpha plot
  
    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
    TH1  &alphaHist = absalpha->GetHist();
    alphaHist.SetXTitle("|alpha|  [\\circ]");
    alphaHist.SetName("alpha-macro");

    Double_t alphasig = 13.1;
    Double_t alphamin = 30.0;
    Double_t alphamax = 90.0;
    Int_t    degree   =    2;
    Double_t significance = -99.0;
    Bool_t   drawpoly  = kTRUE;
    Bool_t   fitgauss  = kTRUE;
    Bool_t   print     = kTRUE;

    MHFindSignificance findsig;
    findsig.SetRebin(kTRUE);
    findsig.SetReduceDegree(kFALSE);

    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
                        alphasig, drawpoly, fitgauss, print);
    significance = findsig.GetSignificance();
    Float_t alphasi = findsig.GetAlphasi();

    gLog << "For file '" << fileName << "' : " << 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 ONOFFAnalysis : End of Job B_RF_UP" << endl;
    gLog << "=======================================================" << endl;
 }
  //---------------------------------------------------------------------


  //---------------------------------------------------------------------
  // Job B_SC_UP
  //============

    //  - create (or read in) optimum supercuts parameter values
    //
    //  - calculate the hadroness for the supercuts
    //
    //  - update input root file, including the hadroness


 if (JobB_SC_UP)
 {
    gLog << "=====================================================" << endl;
    gLog << "Macro ONOFFAnalysis : Start of Job B_SC_UP" << endl;

    gLog << "" << endl;
    gLog << "Macro ONOFFAnalysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = "
         << (JobB_SC_UP ? "kTRUE" : "kFALSE")  << ",  " 
         << (CMatrix    ? "kTRUE" : "kFALSE")  << ",  "
         << (RMatrix    ? "kTRUE" : "kFALSE")  << ",  "
         << (WOptimize  ? "kTRUE" : "kFALSE")  << ",  "
         << (RTest      ? "kTRUE" : "kFALSE")  << ",  "
         << (WSC        ? "kTRUE" : "kFALSE")  << endl;


    //--------------------------------------------
    // file which contains the initial parameter values for the supercuts 
    // if parSCinit ="" the initial values are taken from the constructor of
    //                  MSupercuts

    TString parSCinit = outPath;
    //parSCinit += "parSC_060204";
    //parSCinit = "parSC_240204a";
    parSCinit = "";

    if (parSCinit != "")
    {
      gLog << "Initial values of parameters are taken from file '" 
           << parSCinit << "'" << endl;
    }
    else
    {
      gLog << "Initial values of parameters are taken from the constructor of MSupercuts" 
           << endl;
    }


    //---------------
    // file onto which the optimal parameter values for the supercuts 
    // are written

    TString parSCfile = outPath;
    parSCfile += "parSC_050304a";

    gLog << "parSCfile = " << parSCfile << endl;

    //--------------------------------------------
    // file to be updated (either ON or MC)

    TString typeInput = "ON";
    //TString typeInput = "OFF";
    //TString typeInput = "MC";
    gLog << "typeInput = " << typeInput << endl;

    if (typeInput == "ON")
      TString file(onfile);
    else if (typeInput == "OFF")
      TString file(offfile);
    else if (typeInput == "MC")
      TString file(mcfile);

    // name of input root file
    TString filenameData = outPath;
    filenameData += file;
    filenameData += "Hillas";
    filenameData += typeInput;
    filenameData += "1c.root";
    gLog << "filenameData = " << filenameData << endl; 

    // name of output root file
    TString outNameImage = outPath;
    outNameImage += file;
    outNameImage += "Hillas";
    outNameImage += typeInput;
    outNameImage += "2c.root";
    

    //TString outNameImage = filenameData;

    gLog << "outNameImage = " << outNameImage << endl; 

    //--------------------------------------------
    // files to be read for optimizing the supercuts
    // 
    // for the training
    TString filenameTrain = outPath;
    filenameTrain += onfile;
    filenameTrain += "Hillas";
    filenameTrain += typeInput;
    filenameTrain += "1.root";
    Int_t howManyTrain = 20000;
    gLog << "filenameTrain = " << filenameTrain << ",   howManyTrain = "
         << howManyTrain  << endl; 

    // for testing
    TString filenameTest = outPath;
    filenameTest += onfile;
    filenameTest += "Hillas";
    filenameTest += typeInput;
    filenameTest += "1.root";
    Int_t howManyTest = 20000;

    gLog << "filenameTest = " << filenameTest << ",   howManyTest = "
         << howManyTest  << endl; 
    

    //--------------------------------------------
    // files to contain the matrices (generated from filenameTrain and
    //                                               filenameTest)
    // 
    // for the training
    TString fileMatrixTrain = outPath;
    fileMatrixTrain += "MatrixTrainSC";
    fileMatrixTrain += ".root";
    gLog << "fileMatrixTrain = " << fileMatrixTrain << endl; 

    // for testing
    TString fileMatrixTest = outPath;
    fileMatrixTest += "MatrixTestSC";
    fileMatrixTest += ".root";
    gLog << "fileMatrixTest = " << fileMatrixTest << endl; 

    

    //---------------------------------------------------------------------
    // Training and test matrices :
    // - either create them and write them onto a file
    // - or read them from a file


    MFindSupercuts findsuper;
    findsuper.SetFilenameParam(parSCfile);
    findsuper.SetHadronnessName("HadSC");
    //findsuper.SetUseOrigDistribution(kTRUE);

    //--------------------------
    // create matrices and write them onto files 
    if (CMatrix)
    {
      TString mname("costheta");
      MBinning bin("Binning"+mname);
      bin.SetEdges(10, 0., 1.0);

      MH3 mh3("cos((MPointingPos.fZd)/kRad2Deg)");
      mh3.SetName(mname);
      MH::SetBinning(&mh3.GetHist(), &bin);
      //for (Int_t i=1; i<=mh3.GetNbins(); i++)
      //  mh3.GetHist().SetBinContent(i, 1.0);


      if (filenameTrain == filenameTest)
      {
        if ( !findsuper.DefineTrainTestMatrix(
                              filenameTrain,   mh3, 
                              howManyTrain,    howManyTest,  
                              fileMatrixTrain, fileMatrixTest)     )
        {
          *fLog << "MagicAnalysis.C : DefineTrainTestMatrix failed" << endl;
          return;
        }

      }
      else
      {
        if ( !findsuper.DefineTrainMatrix(filenameTrain, mh3,
                                          howManyTrain,  fileMatrixTrain) )
        {
          *fLog << "MagicAnalysis.C : DefineTrainMatrix failed" << endl;
          return;
        }

	if ( !findsuper.DefineTestMatrix( filenameTest, mh3, 
                                          howManyTest,  fileMatrixTest)  )
        {
          *fLog << "MagicAnalysis.C : DefineTestMatrix failed" << endl;
          return;
        }
      }
     }

    //--------------------------
    // read matrices from files
    //                              

    if (RMatrix)
      findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
    //--------------------------



    //---------------------------------------------------------------------
    // optimize supercuts using the training sample
    // 
    // the initial values are taken 
    //     - from the file parSCinit (if != "")
    //     - or from the arrays params and steps (if their sizes are != 0)
    //     - or from the MSupercuts constructor


if (WOptimize)
  {
    gLog << "MagicAnalysis.C : optimize the supercuts using the training matrix" 
         << endl;

    TArrayD params(0);
    TArrayD steps(0);
  

    if (parSCinit == "")
    {
      /*
      Double_t vparams[104] = {
      // LengthUp
	0.2,       0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // LengthLo
        0.,        0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // WidthUp
        0.1,       0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // WidthLo
        0.,        0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // DistUp
        0.,        0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // DistLo
        0.,        0.,       0.,       0.,        0.,        0.0,
	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.02,      0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // LengthLo
        0.,        0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // WidthUp
        0.01,      0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // WidthLo
        0.,        0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // DistUp
        0.,        0.,       0.,       0.,        0.,        0.0,
	0.,        0.,
      // DistLo
        0.,        0.,       0.,       0.,        0.,        0.0,
	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                                                 };
      */


      Double_t vparams[104] = {
      // LengthUp
	0.315585,  0.001455, 0.203198, 0.005532, -0.001670, -0.020362,
	0.007388, -0.013463,
      // LengthLo
        0.151530,  0.028323, 0.510707, 0.053089,  0.013708,  2.357993,
	0.000080, -0.007157,
      // WidthUp
        0.145412, -0.001771, 0.054462, 0.022280, -0.009893,  0.056353,
        0.020711, -0.016703,
      // WidthLo
        0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101,
        0.006126, -0.002849,
      // DistUp
        1.787943,  0.0,      2.942310, 0.199815,  0.0,       0.249909,
        0.189697,  0.0,
      // DistLo
        0.589406,  0.0,     -0.083964,-0.007975,  0.0,       0.045374,
       -0.001750,  0.0,
      // AsymUp
        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // AsymLo
       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // ConcUp
        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // ConcLo
       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // Leakage1Up
        1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // Leakage1Lo
       -1.e10,     0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // AlphaUp
	13.12344,  0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0                                                 };

      Double_t vsteps[104] = {
      // LengthUp
        0.03,      0.0002,   0.02,     0.0006,    0.0002,    0.002,
        0.0008,    0.002,
      // LengthLo
        0.02,      0.003,    0.05,     0.006,     0.002,     0.3,
        0.0001,    0.0008,
      // WidthUp
        0.02,      0.0002,   0.006,    0.003,     0.002,     0.006,
        0.002,     0.002,
      // WidthLo
        0.009,     0.0007,   0.008,    0.0004,    0.0005,    0.002,
        0.0007,    0.003,
      // DistUp
        0.2,       0.0,      0.3,      0.02,      0.0,       0.03,
        0.02,      0.0
      // DistLo
        0.06,      0.0,      0.009,    0.0008,    0.0,       0.005,
        0.0002,    0.0
      // AsymUp  
        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // AsymLo  
        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // ConcUp  
        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // ConcLo  
        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // Leakage1Up  
        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // Leakage1Lo  
        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0,
      // AlphaUp  
        0.0,       0.0,      0.0,      0.0,       0.0,       0.0,
        0.0,       0.0                                                 };

 
      params.Set(104, vparams);
      steps.Set (104, vsteps );
    }

    Bool_t rf;
    rf = findsuper.FindParams(parSCinit, params, steps);

    if (!rf) 
    {
       gLog << "MagicAnalysis.C : optimization of supercuts failed" << endl;
       return;
    }
  }

    //--------------------------------------
    // test the supercuts on the test sample
    //    

 if (RTest)
 {
    gLog << "MagicAnalysis.C : test the supercuts on the test matrix" << endl;
    Bool_t rt = findsuper.TestParams();
    if (!rt) 
    {
       gLog << "MagicAnalysis.C : test of supercuts on the test matrix failed" 
            << endl;
    }

 }


    //-----------------------------------------------------------------
    // Update the input files with the SC hadronness
    //

 if (WSC)
 {
    gLog << "" << endl;
    gLog << "========================================================" << endl;
    gLog << "Update input file '" <<  filenameData 
         << "' with the SC hadronness" << endl;


    //----------------------------------------------------
    // read in optimum parameter values for the supercuts

    TFile inparam(parSCfile);
    MSupercuts scin;
    scin.Read("MSupercuts");
    inparam.Close();

    gLog << "Parameter values for supercuts were read in from file '"
         << parSCfile << "'" << endl;

    TArrayD supercutsPar;
    supercutsPar =  scin.GetParameters();

    TArrayD supercutsStep;
    supercutsStep =  scin.GetStepsizes();

    gLog << "Parameter values for supercuts : " << endl;
    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
    {
      gLog << supercutsPar[i] << ",  ";
    }
    gLog << endl;

    gLog << "Step values for supercuts : " << endl;
    for (Int_t i=0; i<supercutsStep.GetSize(); i++)
    {
      gLog << supercutsStep[i] << ",  ";
    }
    gLog << endl;


    //----------------------------------------------------
    MTaskList tliston;
    MParList pliston;

    // set the parameters of the supercuts
    MSupercuts supercuts;
    supercuts.SetParameters(supercutsPar);
    gLog << "parameter values for the supercuts used for updating the input file ' " 
         << filenameData << "'" << endl;
    supercutsPar = supercuts.GetParameters();
    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
    {
      gLog << supercutsPar[i] << ",  ";
    }
    gLog << endl;


    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "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"; 


    //.......................................................................
    // calculation of hadroness for the supercuts
    // (=0.25 if fullfilled, =0.75 otherwise)

    TString hadSCName = "HadSC";
    MSupercutsCalc sccalc(fHilName, fHilNameSrc);
    sccalc.SetHadronnessName(hadSCName);


    //.......................................................................


      //MWriteRootFile write(outNameImage, "UPDATE");
      //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");

    
      MWriteRootFile write(outNameImage, "RECREATE");

      write.AddContainer("MRawRunHeader", "RunHeaders");
      write.AddContainer("MTime",         "Events");
      write.AddContainer("MPointingPos",  "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("HadRF",         "Events");
      write.AddContainer(hadSCName,       "Events");
    

    //-----------------------------------------------------------------
    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");

    Float_t maxhadronness =  0.40;
    Float_t maxalpha      =  20.0;
    Float_t maxdist       =  10.0;

    MFSelFinal selfinalgh(fHilNameSrc);
    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    selfinalgh.SetHadronnessName(hadSCName);
    selfinalgh.SetName("SelFinalgh");
    MContinue contfinalgh(&selfinalgh);
    contfinalgh.SetName("ContSelFinalgh");

    MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
    fillhadsc.SetName("HhadSC");

    MFSelFinal selfinal(fHilNameSrc);
    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    selfinal.SetHadronnessName(hadSCName);
    selfinal.SetName("SelFinal");
    MContinue contfinal(&selfinal);
    contfinal.SetName("ContSelFinal");

    TString mh3name = "abs(Alpha)";
    MBinning binsalphaabs("Binning"+mh3name);
    binsalphaabs.SetEdges(50, -2.0, 98.0);

    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
    alphaabs.SetName(mh3name);

    TH1  &alphahist = alphaabs->GetHist();

    MFillH alpha(&alphaabs);
    alpha.SetName("FillAlphaAbs");

    MFillH hfill1("MHHillas",    fHilName);
    hfill1.SetName("HHillas");

    MFillH hfill2("MHStarMap",   fHilName);
    hfill2.SetName("HStarMap");

    MFillH hfill3("MHHillasExt",    fHilNameSrc);
    hfill3.SetName("HHillasExt");
    
    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    hfill4.SetName("HHillasSrc");    

    MFillH hfill5("MHNewImagePar", fImgParName);
    hfill5.SetName("HNewImagePar");

    //*****************************
    // entries in MParList

    pliston.AddToList(&tliston);
    InitBinnings(&pliston);

    pliston.AddToList(&supercuts);

    pliston.AddToList(&binsalphaabs);
    pliston.AddToList(&alphaabs);

    //*****************************
    // entries in MTaskList
    
    tliston.AddToList(&read);

    tliston.AddToList(&sccalc);
    tliston.AddToList(&fillhadsc);

    tliston.AddToList(&write);
    tliston.AddToList(&contfinalgh);

    tliston.AddToList(&alpha);
    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    tliston.AddToList(&contfinal);

    //*****************************

    //-------------------------------------------
    // Execute event loop
    //
    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&pliston);
    evtloop.SetProgressBar(&bar);

    Int_t maxevents = -1;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // Display the histograms
    //
    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();

    pliston.FindObject("MHHillas")->DrawClone();
    pliston.FindObject("MHHillasExt")->DrawClone();
    pliston.FindObject("MHHillasSrc")->DrawClone();
    pliston.FindObject("MHNewImagePar")->DrawClone();
    pliston.FindObject("MHStarMap")->DrawClone();

     //-------------------------------------------
    // fit alpha distribution to get the number of excess events and
    // calculate significance of gamma signal in the alpha plot
  
    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
    TH1  &alphaHist = absalpha->GetHist();
    alphaHist.SetXTitle("|alpha|  [\\circ]");
    alphaHist.SetName("alpha-macro");

    Double_t alphasig = 13.1;
    Double_t alphamin = 30.0;
    Double_t alphamax = 90.0;
    Int_t    degree   =    2;
    Double_t significance = -99.0;
    Bool_t   drawpoly  = kTRUE;
    Bool_t   fitgauss  = kTRUE;
    Bool_t   print     = kTRUE;

    MHFindSignificance findsig;
    findsig.SetRebin(kTRUE);
    findsig.SetReduceDegree(kFALSE);

    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
                        alphasig, drawpoly, fitgauss, print);
    significance = findsig.GetSignificance();
    Float_t alphasi = findsig.GetAlphasi();

    gLog << "For file '" << filenameData << "' : " << endl;
    gLog << "Significance of gamma signal after supercuts : "
         << significance << " (for |alpha| < " << alphasi << " degrees)" 
         << endl;

    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);

    //-------------------------------------------

    DeleteBinnings(&pliston);
 }


    gLog << "Macro ONOFFAnalysis : End of Job B_SC_UP" << endl;
    gLog << "=======================================================" << endl;
 }
  //---------------------------------------------------------------------



  //---------------------------------------------------------------------
  // Job C  
  //======

    //  - read ON1 and MC1 data files  
    //    which should have been updated to contain the hadronnesses
    //    for the method of Random Forest and for the SUPERCUTS
    //  - produce Neyman-Pearson plots
 
 if (JobC)
 {
    gLog << "=====================================================" << endl;
    gLog << "Macro ONOFFAnalysis : Start of Job C" << endl;

    gLog << "" << endl;
    gLog << "Macro ONOFFAnalysis : JobC = " 
         << (JobC ? "kTRUE" : "kFALSE")  << endl;



    TString ext("2.root");
    TString extout("2.root");

    TString typeHadrons("OFF");
    TString typeGammas("MC");

    //--------------------------------------------
    // name of input data file
    TString NameData = outPath;
    NameData += typeHadrons;
    TString inNameData(NameData);
    inNameData += ext;
    gLog << "inNameData   = " << inNameData   << endl;

    // name of input MC file
    TString NameMC = outPath;
    NameMC += typeGammas;
    TString inNameMC(NameMC);
    inNameMC += ext;
    gLog << "inNameMC   = " << inNameMC   << endl;


    //--------------------------------------------
    // root files for the training events 



    TString NameGammasTrain = outPath;
    NameGammasTrain += "RF_gammas_Train_";
    NameGammasTrain += typeGammas;
    TString outNameGammasTrain(NameGammasTrain);    
    outNameGammasTrain += extout;


    TString NameHadronsTrain = outPath;
    NameHadronsTrain += "RF_hadrons_Train_";
    NameHadronsTrain += typeHadrons;
    TString outNameHadronsTrain(NameHadronsTrain);    
    outNameHadronsTrain += extout;


    //--------------------------------------------
    // root files for the test events 

    TString NameGammasTest = outPath;
    NameGammasTest += "RF_gammas_Test_";
    NameGammasTest += typeGammas;
    TString outNameGammasTest(NameGammasTest);    
    outNameGammasTest += extout;

    TString NameHadronsTest = outPath;
    NameHadronsTest += "RF_hadrons_Test_";
    NameHadronsTest += typeHadrons;
    TString outNameHadronsTest(NameHadronsTest);    
    outNameHadronsTest += extout;

    //--------------------------------------------------------------------

    //TString filenameData(inNameData);
    //TString filenameMC(inNameMC);

    //TString filenameData(outNameHadronsTrain);
    //TString filenameMC(outNameGammasTrain);

    TString filenameData(outNameHadronsTest);
    TString filenameMC(outNameGammasTest);

    gLog << "filenameData = " << filenameData << endl;
    gLog << "filenameMC   = " << filenameMC   << endl;

    //-----------------------------------------------------------------

    MTaskList tliston;
    MParList pliston;


    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");

    //-------------------------------------------
    // create the tasks which should be executed 
    //

    MReadMarsFile read("Events", filenameMC);
    read.AddFile(filenameData);
    read.DisableAutoScheme();


    //.......................................................................
    // names of hadronness containers

    TString hadSCName = "HadSC";
    TString hadRFName = "HadRF";

    //.......................................................................


    TString fHilName    = "MHillas"; 
    TString fHilNameExt = "MHillasExt"; 
    TString fHilNameSrc = "MHillasSrc"; 
    TString fImgParName = "MNewImagePar"; 

    Float_t maxhadronness =  0.40;
    Float_t maxalpha      =  20.0;
    Float_t maxdist       =  10.0;

    MFSelFinal selfinalgh(fHilNameSrc);
    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    selfinalgh.SetHadronnessName(hadRFName);
    selfinalgh.SetName("SelFinalgh");
    MContinue contfinalgh(&selfinalgh);
    contfinalgh.SetName("ContSelFinalgh");


    //MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
    //fillhadsc.SetName("HhadSC");
    MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
    fillhadrf.SetName("HhadRF");

    MFSelFinal selfinal(fHilNameSrc);
    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    selfinal.SetHadronnessName(hadRFName);
    selfinal.SetName("SelFinal");
    MContinue contfinal(&selfinal);
    contfinal.SetName("ContSelFinal");


    MFillH hfill1("MHHillas",    fHilName);
    hfill1.SetName("HHillas");

    MFillH hfill2("MHStarMap",   fHilName);
    hfill2.SetName("HStarMap");

    MFillH hfill3("MHHillasExt",    fHilNameSrc);
    hfill3.SetName("HHillasExt");
    
    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    hfill4.SetName("HHillasSrc");    

    MFillH hfill5("MHNewImagePar", fImgParName);
    hfill5.SetName("HNewImagePar");


    //*****************************
    // entries in MParList

    pliston.AddToList(&tliston);
    InitBinnings(&pliston);


    //*****************************
    // entries in MTaskList
    
    tliston.AddToList(&read);

    //tliston.AddToList(&fillhadsc);
    tliston.AddToList(&fillhadrf);
   
    tliston.AddToList(&contfinalgh);
    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    tliston.AddToList(&contfinal);

    //*****************************

    //-------------------------------------------
    // Execute event loop
    //
    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&pliston);
    evtloop.SetProgressBar(&bar);

    Int_t maxevents = -1;
    //Int_t maxevents = 35000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // Display the histograms
    //

    //pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();

    pliston.FindObject("MHHillas")->DrawClone();
    pliston.FindObject("MHHillasExt")->DrawClone();
    pliston.FindObject("MHHillasSrc")->DrawClone();
    pliston.FindObject("MHNewImagePar")->DrawClone();
    pliston.FindObject("MHStarMap")->DrawClone();

    DeleteBinnings(&pliston);

    gLog << "Macro ONOFFAnalysis : End of Job C" << endl;
    gLog << "===================================================" << endl;
 }



  //---------------------------------------------------------------------
  // Job D
  //======

    //  - select g/h separation method XX
    //  - read ON2 (or MC2) root file 
    //  - apply cuts in hadronness
    //  - make plots


 if (JobD)
 {
    gLog << "=====================================================" << endl;
    gLog << "Macro ONOFFAnalysis : Start of Job D" << endl;

    gLog << "" << endl;
    gLog << "Macro ONOFFAnalysis : JobD = " 
         << (JobD        ? "kTRUE" : "kFALSE")  << endl;


    // type of data to be analysed
    TString typeData = "ON";
    //TString typeData = "OFF";
    //TString typeData = "MC";
    gLog << "typeData = " << typeData << endl;

    TString ext      = "3.root";


    //------------------------------
    // selection of g/h separation method
    // and definition of final selections

    //TString XX("SC");
    TString XX("RF");
    TString fhadronnessName("Had");
    fhadronnessName += XX;
    gLog << "fhadronnessName = " << fhadronnessName << endl;

    // maximum values of the hadronness, |ALPHA| and DIST
    Float_t maxhadronness   = 0.233;
    Float_t maxalpha        = 20.0;
    Float_t maxdist         = 10.0;
    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
         << maxhadronness << ",  " << maxalpha << ",  " 
         << maxdist << endl;


    //------------------------------
    // name of data file to be analysed
    TString filenameData(outPath);
    filenameData += typeData;
    filenameData += ext;
    gLog << "filenameData = " << filenameData << endl;



    //*************************************************************************
    //
    // Analyse the data
    //

    MTaskList tliston;
    MParList pliston;

    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");


    TString fHilName    = "MHillas"; 
    TString fHilNameExt = "MHillasExt"; 
    TString fHilNameSrc = "MHillasSrc"; 
    TString fImgParName = "MNewImagePar"; 

    //-------------------------------------------
    // create the tasks which should be executed 
    //

    MReadMarsFile read("Events", filenameData);
    read.DisableAutoScheme();


    //-----------------------------------------------------------------
    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");

    MFSelFinal selfinalgh(fHilNameSrc);
    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    selfinalgh.SetHadronnessName(fhadronnessName);
    selfinalgh.SetName("SelFinalgh");
    MContinue contfinalgh(&selfinalgh);
    contfinalgh.SetName("ContSelFinalgh");

    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
    fillhadsc.SetName("HhadSC");
    MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
    fillhadrf.SetName("HhadRF");

    TString mh3name = "abs(Alpha)";
    MBinning binsalphaabs("Binning"+mh3name);
    binsalphaabs.SetEdges(50, -2.0, 98.0);

    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
    alphaabs.SetName(mh3name);

    TH1  &alphahist = alphaabs->GetHist();

    MFillH alpha(&alphaabs);
    alpha.SetName("FillAlphaAbs");

    MFillH hfill1("MHHillas",    fHilName);
    hfill1.SetName("HHillas");

    MFillH hfill2("MHStarMap",   fHilName);
    hfill2.SetName("HStarMap");

    MFillH hfill3("MHHillasExt",   fHilNameSrc);
    hfill3.SetName("HHillasExt");    

    MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    hfill4.SetName("HHillasSrc");    

    MFillH hfill5("MHNewImagePar", fImgParName);
    hfill5.SetName("HNewImagePar");    

    MFSelFinal selfinal(fHilNameSrc);
    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    selfinal.SetHadronnessName(fhadronnessName);
    selfinal.SetName("SelFinal");
    MContinue contfinal(&selfinal);
    contfinal.SetName("ContSelFinal");


    //*****************************
    // entries in MParList

    pliston.AddToList(&tliston);
    InitBinnings(&pliston);
    pliston.AddToList(&binsalphaabs);
    pliston.AddToList(&alphaabs);

    //*****************************
    // entries in MTaskList
    
    tliston.AddToList(&read);

    tliston.AddToList(&contfinalgh);

    tliston.AddToList(&fillhadsc);
    tliston.AddToList(&fillhadrf);

    tliston.AddToList(&alpha);
    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    tliston.AddToList(&contfinal);

    //*****************************

    //-------------------------------------------
    // Execute event loop
    //
    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&pliston);
    evtloop.SetProgressBar(&bar);

    Int_t maxevents = -1;
    //Int_t maxevents = 10000;
    if ( !evtloop.Eventloop(maxevents) )
        return;

    tliston.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // Display the histograms
    //

    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();

    pliston.FindObject("MHHillas")->DrawClone();
    pliston.FindObject("MHHillasExt")->DrawClone();
    pliston.FindObject("MHHillasSrc")->DrawClone();
    pliston.FindObject("MHNewImagePar")->DrawClone();
    pliston.FindObject("MHStarMap")->DrawClone();


    //-------------------------------------------

    // fit alpha distribution to get the number of excess events and
    // calculate significance of gamma signal in the alpha plot
  
    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
    TH1  &alphaHist = absalpha->GetHist();
    alphaHist.SetXTitle("|alpha|  [\\circ]");
    alphaHist.SetName("alpha-JobD");

    Double_t alphasig = 13.1;
    Double_t alphamin = 30.0;
    Double_t alphamax = 90.0;
    Int_t    degree   =    2;
    Double_t significance = -99.0;
    Bool_t   drawpoly  = kTRUE;
    Bool_t   fitgauss  = kTRUE;
    Bool_t   print     = kTRUE;

    MHFindSignificance findsig;
    findsig.SetRebin(kTRUE);
    findsig.SetReduceDegree(kFALSE);

    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
                        alphasig, drawpoly, fitgauss, print);
    significance = findsig.GetSignificance();
    Float_t alphasi = findsig.GetAlphasi();

    gLog << "For file '" << filenameData << "' : " << endl;
    gLog << "Significance of gamma signal after supercuts : "
         << significance << " (for |alpha| < " << alphasi << " degrees)" 
         << endl;

    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);

    //-------------------------------------------


    DeleteBinnings(&pliston);

    gLog << "Macro ONOFFAnalysis : End of Job D" << endl;
    gLog << "=======================================================" << endl;
 }
  //---------------------------------------------------------------------





  //---------------------------------------------------------------------
  // Job E_XX
  //=========

    //  - select g/h separation method XX
    //  - read MC_XX2.root file 
    //  - calculate eff. collection area
    //  - read ON_XX2.root file 
    //  - apply final cuts
    //  - calculate flux
    //  - write root file for ON data after final cuts (ON_XX3.root))


 if (JobE_XX)
 {
    gLog << "=====================================================" << endl;
    gLog << "Macro ONOFFAnalysis : Start of Job E_XX" << endl;

    gLog << "" << endl;
    gLog << "Macro ONOFFAnalysis : JobE_XX, CCollArea, OEEst, WEX = " 
         << (JobE_XX ? "kTRUE" : "kFALSE")  << ",  " 
         << (CCollArea?"kTRUE" : "kFALSE")  << ",  " 
         << (OEEst ?   "kTRUE" : "kFALSE")  << ",  " 
         << (WEX     ? "kTRUE" : "kFALSE")  << endl;


    // type of data to be analysed
    //TString typeData = "ON";
    //TString typeData = "OFF";
    TString typeData = "MC";
    gLog << "typeData = " << typeData << endl;

    TString typeMC   = "MC";
    TString ext      = "3.root";
    TString extout   = "4.root";

    //------------------------------
    // selection of g/h separation method
    // and definition of final selections

    //TString XX("SC");
    TString XX("RF");
    TString fhadronnessName("Had");
    fhadronnessName += XX;
    gLog << "fhadronnessName = " << fhadronnessName << endl;

    // maximum values of the hadronness, |ALPHA| and DIST
    Float_t maxhadronness   = 0.23;
    Float_t maxalpha        = 20.0;
    Float_t maxdist         = 10.0;
    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
         << maxhadronness << ",  " << maxalpha << ",  " 
         << maxdist << endl;

    //------------------------------
    // name of MC file to be used for optimizing the energy estimator
    TString filenameOpt(outPath);
    filenameOpt += typeMC;
    filenameOpt += ext; 
    gLog << "filenameOpt = " << filenameOpt << endl;

    //------------------------------
    // name of file containing the parameters of the energy estimator
    TString energyParName(outPath);
    energyParName += "energyest_";
    energyParName += XX;
    energyParName += ".root";
    gLog << "energyParName = " << energyParName << endl;

    //------------------------------
    // name of MC file to be used for calculating the eff. collection areas
    TString filenameArea(outPath);
    filenameArea += typeMC;
    filenameArea += ext; 
    gLog << "filenameArea = " << filenameArea << endl;

    //------------------------------
    // name of file containing the eff. collection areas
    TString collareaName(outPath);
    collareaName += "area_";
    collareaName += XX;
    collareaName += ".root";
    gLog << "collareaName = " << collareaName << endl;

    //------------------------------
    // name of data file to be analysed
    TString filenameData(outPath);
    filenameData += typeData;
    filenameData += ext;
    gLog << "filenameData = " << filenameData << endl;

    //------------------------------
    // name of output data file (after the final cuts)
    TString filenameDataout(outPath);
    filenameDataout += typeData;
    filenameDataout += "_";
    filenameDataout += XX;
    filenameDataout += extout;
    gLog << "filenameDataout = " << filenameDataout << endl;

    //------------------------------
    // name of file containing histograms for flux calculastion
    TString filenameResults(outPath);
    filenameResults += typeData;
    filenameResults += "Results_";
    filenameResults += XX;
    filenameResults += extout;
    gLog << "filenameResults = " << filenameResults << endl;


    //====================================================================

    MHMcCollectionArea collarea;
    collarea.SetEaxis(MHMcCollectionArea::kLinear);

    MParList  parlist;
    InitBinnings(&parlist);

  if (CCollArea)
  {
    gLog << "-----------------------------------------------" << endl;
    gLog << "Start calculation of effective collection areas" << endl;


    MTaskList tasklist;

    //---------------------------------------
    // Setup the tasks to be executed
    //
    MReadMarsFile reader("Events", filenameArea);
    reader.DisableAutoScheme();

    MFSelFinal cuthadrons;
    cuthadrons.SetHadronnessName(fhadronnessName);
    cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);

    MContinue conthadrons(&cuthadrons);


    MFillH filler("MHMcCollectionArea", "MPointingPos");
    filler.SetName("CollectionArea");

    //********************************
    // entries in MParList

    parlist.AddToList(&tasklist);

    parlist.AddToList(&collarea);

    //********************************
    // entries in MTaskList

    tasklist.AddToList(&reader);   
    tasklist.AddToList(&conthadrons);
    tasklist.AddToList(&filler);

    //********************************

    //-----------------------------------------
    // Execute event loop
    //
    MEvtLoop magic;
    magic.SetParList(&parlist);

    MProgressBar bar;
    magic.SetProgressBar(&bar);
    if (!magic.Eventloop())
        return;

    tasklist.PrintStatistics(0, kTRUE);

    // Calculate effective collection areas 
    // and display the histograms
    //
    //MHMcCollectionArea *collarea = 
    //     (MHMcCollectionArea*)parlist.FindObject("MHMcCollectionArea");
    collarea.CalcEfficiency();
    collarea.DrawClone();



    //---------------------------------------------
    // Write histograms to a file 
    //

    TFile f(collareaName, "RECREATE");
    //collarea.GetHist()->Write();
    //collarea.GetHAll()->Write();
    //collarea.GetHSel()->Write();
    collarea.Write();

    f.Close();

    gLog << "Collection area plots written onto file " << collareaName << endl;

    gLog << "Calculation of effective collection areas done" << endl;
    gLog << "-----------------------------------------------" << endl;    
    //------------------------------------------------------------------
  }

  if (!CCollArea)
  {
    gLog << "-----------------------------------------------" << endl;
    gLog << "Read in effective collection areas from file " 
         << collareaName << endl;

    TFile collfile(collareaName);
    collfile.ls();
    collarea.Read("MHMcCollectionArea");
    collarea.DrawClone();

    gLog << "Effective collection areas were read in from file " 
         << collareaName << endl;
    gLog << "-----------------------------------------------" << endl;    
  }


    // save binnings for call to MagicEEst
    MBinning *binsE     = (MBinning*)parlist.FindObject("BinningE");
    if (!binsE)
	{
          gLog << "Object 'BinningE' not found in MParList" << endl;
          return;
	}
    MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
    if (!binsTheta)
	{
          gLog << "Object 'BinningTheta' not found in MParList" << endl;
          return;
	}

    //-------------------------------------
    TString fHilName    = "MHillas"; 
    TString fHilNameExt = "MHillasExt"; 
    TString fHilNameSrc = "MHillasSrc"; 
    TString fImgParName = "MNewImagePar"; 


 if (OEEst)
 { 
   //===========================================================
    //
    // Optimization of energy estimator
    //
    gLog << "Macro ONOFFAnalysis.C : calling MagicEEst" << endl;

    TString inpath("");
    TString outpath("");
    Int_t howMany = 2000;
    MagicEEst(inpath,   filenameOpt,   outpath, energyParName, 
            fHilName, fHilNameSrc,   fhadronnessName,
            howMany,  maxhadronness, maxalpha, maxdist,
            binsE, binsTheta);
    gLog << "Macro ONOFFAnalysis.C : returning from MagicEEst" << endl;
 }

 if (WEX)
 {
    //-----------------------------------------------------------
    //
    // Read in parameters of energy estimator ("MMcEnergyEst")
    //                   and migration matrix ("MHMcEnergyMigration")
    //
    gLog << "================================================================"
         << endl;
    gLog << "Macro ONOFFAnalysis.C : read parameters of energy estimator and migration matrix from file '"
         << energyParName << "'" << endl;
    TFile enparam(energyParName);
    enparam.ls();
    MMcEnergyEst mcest("MMcEnergyEst"); 
    mcest.Read("MMcEnergyEst");

    //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
    gLog << "Parameters of energy estimator were read in" << endl;


    gLog << "Read in Migration matrix" << endl;   

    MHMcEnergyMigration mighiston("MHMcEnergyMigration");
    mighiston.Read("MHMcEnergyMigration");
    //MHMcEnergyMigration &mighiston = 
    //      *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));

    gLog << "Migration matrix was read in" << endl;


    TArrayD parA(mcest.GetNumCoeffA());
    TArrayD parB(mcest.GetNumCoeffB());
    for (Int_t i=0; i<parA.GetSize(); i++)
      parA[i] = mcest.GetCoeff(i);
    for (Int_t i=0; i<parB.GetSize(); i++)
      parB[i] = mcest.GetCoeff( i+parA.GetSize() );

    //*************************************************************************
    //
    // Analyse the data
    //
    gLog << "============================================================"
         << endl;
    gLog << "Analyse the data" << endl;

    MTaskList tliston;
    MParList pliston;

    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");


    //-------------------------------------------
    // create the tasks which should be executed 
    //

    MReadMarsFile read("Events", filenameData);
    read.DisableAutoScheme();

    //.......................................................................

      gLog << "MagicAnalysis.C : write root file '" << filenameDataout 
           << "'" << endl;
   
      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));


      MWriteRootFile write(filenameDataout, "RECREATE");

      write.AddContainer("MRawRunHeader", "RunHeaders");
      write.AddContainer("MTime",         "Events");
      write.AddContainer("MPointingPos",  "Events");
      write.AddContainer("ThetaOrig",     "Events");
      write.AddContainer("MSrcPosCam",    "Events");
      write.AddContainer("MSigmabar",     "Events");
      write.AddContainer("MHillas",       "Events");
      write.AddContainer("MHillasExt",    "Events");
      write.AddContainer("MHillasSrc",    "Events");
      write.AddContainer("MNewImagePar",  "Events");

      //write.AddContainer("HadNN",         "Events");
      write.AddContainer("HadSC",         "Events");
      write.AddContainer("HadRF",         "Events");

      write.AddContainer("MEnergyEst",    "Events");


    //-----------------------------------------------------------------
    // geometry is needed in  MHHillas... classes 
    MGeomCam *fGeom = 
             (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");

    MFSelFinal selfinalgh(fHilNameSrc);
    selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    selfinalgh.SetHadronnessName(fhadronnessName);
    selfinalgh.SetName("SelFinalgh");
    MContinue contfinalgh(&selfinalgh);
    contfinalgh.SetName("ContSelFinalgh");

    //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
    //fillhadnn.SetName("HhadNN");
    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
    fillhadsc.SetName("HhadSC");
    MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
    fillhadrf.SetName("HhadRF");

    //---------------------------
    // calculate estimated energy

    MEnergyEstParam eeston(fHilName);
    eeston.Add(fHilNameSrc);

    eeston.SetCoeffA(parA);
    eeston.SetCoeffB(parB);

    //---------------------------
    // calculate estimated energy using Daniel's parameters

    //MEnergyEstParamDanielMkn421 eeston(fHilName);
    //eeston.Add(fHilNameSrc);
    //eeston.SetCoeffA(parA);
    //eeston.SetCoeffB(parB);


    //---------------------------


    MFillH hfill1("MHHillas",    fHilName);
    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");    

    //---------------------------
    // new from Robert

    MFillH hfill6("MHTimeDiffTheta", "MPointingPos");
    hfill6.SetName("HTimeDiffTheta");

    MFillH hfill6a("MHTimeDiffTime", "MPointingPos");
    hfill6a.SetName("HTimeDiffTime");

    MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
    hfill7.SetName("HAlphaEnergyTheta");

    MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
    hfill7a.SetName("HAlphaEnergyTime");

    MFillH hfill7b("MHThetabarTime", fHilNameSrc);
    hfill7b.SetName("HThetabarTime");

    MFillH hfill7c("MHEnergyTime", "MMcEvt");
    hfill7c.SetName("HEnergyTime");


    //---------------------------

    MFSelFinal selfinal(fHilNameSrc);
    selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    selfinal.SetHadronnessName(fhadronnessName);
    selfinal.SetName("SelFinal");
    MContinue contfinal(&selfinal);
    contfinal.SetName("ContSelFinal");


    //*****************************
    // entries in MParList

    pliston.AddToList(&tliston);
    InitBinnings(&pliston);


    //*****************************
    // entries in MTaskList
    
    tliston.AddToList(&read);

    // robert      
    tliston.AddToList(&hfill6);   //timediff
    tliston.AddToList(&hfill6a);   //timediff

    tliston.AddToList(&contfinalgh);
    tliston.AddToList(&eeston);

    tliston.AddToList(&write);

    //tliston.AddToList(&fillhadnn);
    tliston.AddToList(&fillhadsc);
    tliston.AddToList(&fillhadrf);

    tliston.AddToList(&hfill1);
    tliston.AddToList(&hfill2);
    tliston.AddToList(&hfill3);
    tliston.AddToList(&hfill4);
    tliston.AddToList(&hfill5);

    //robert
    tliston.AddToList(&hfill7);
    tliston.AddToList(&hfill7a);
    tliston.AddToList(&hfill7b);
    tliston.AddToList(&hfill7c);

    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("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;

    TFile outfile(filenameResults,"recreate");

    MHHillasSrc* hillasSrc = 
      (MHHillasSrc*)(pliston->FindObject("MHHillasSrc"));
        TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
    alphaHist->Write();
    gLog << "Alpha plot has been written out" << endl;    


    MHAlphaEnergyTheta* aetH = 
      (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
    TH3D* aetHist = (TH3D*)(aetH->GetHist());
    aetHist->SetName("aetHist");
    aetHist->Write();
    gLog << "AlphaEnergyTheta plot has been written out" << endl;    

    MHAlphaEnergyTime* aetH2 = 
      (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime"));
    TH3D* aetHist2 = (TH3D*)(aetH2->GetHist());
    aetHist2->SetName("aetimeHist");
//     aetHist2->DrawClone();
    aetHist2->Write();
    gLog << "AlphaEnergyTime plot has been written out" << endl;    

    MHThetabarTime* tbt = 
      (MHThetabarTime*)(pliston->FindObject("MHThetabarTime"));
    TProfile* tbtHist = (TProfile*)(tbt->GetHist());
    tbtHist->SetName("tbtHist");
    tbtHist->Write();
    gLog << "ThetabarTime plot has been written out" << endl;    

    MHEnergyTime* ent = 
      (MHEnergyTime*)(pliston->FindObject("MHEnergyTime"));
    TH2D* entHist = (TH2D*)(ent->GetHist());
    entHist->SetName("entHist");
    entHist->Write();
    gLog << "EnergyTime plot has been written out" << endl;    
    
    MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
    TH2D* timeHist = (TH2D*)(time->GetHist());
    timeHist->SetName("MHTimeDiffTheta");
    timeHist->SetTitle("Time diffs");
    timeHist->Write();
    gLog << "TimeDiffTheta plot has been written out" << endl;    


    MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
    TH2D* timeHist2 = (TH2D*)(time2->GetHist());
    timeHist2->SetName("MHTimeDiffTime");
    timeHist2->SetTitle("Time diffs");
    timeHist2->Write();
    gLog << "TimeDiffTime plot has been written out" << endl;    

//rwagner write also collareas to same file
    collarea->GetHist()->Write();
    collarea->GetHAll()->Write();
    collarea->GetHSel()->Write();
    gLog << "Effective collection areas have been written out" << endl;        

//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data
//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file

    gLog << "before closing outfile" << endl;

    //outfile.Close();
    gLog << "Results were written onto file '" << filenameResults 
         << "'" << endl;
    gLog << "=======================================================" << endl;
  }

  }

    gLog << "Macro ONOFFAnalysis : End of Job E_XX" << endl;
    gLog << "=======================================================" << endl;
 }
  //---------------------------------------------------------------------

}







   





















