#include <TROOT.h>
#include <TClass.h>
#include <TSystem.h>
#include <TGClient.h>
#include <TApplication.h>
#include <TObjectTable.h>

#include "MHMatrix.h"
#include "MParList.h"
#include "MTaskList.h"
#include "MWriteRootFile.h"
#include "MReadMarsFile.h"
#include "MReadReports.h"
#include "MApplySupercuts.h"
#include "MHFindSignificance.h"
#include "MEvtLoop.h"
#include "MProgressBar.h"
#include "MContinue.h"
#include "MGeomApply.h"
#include "MHadronness.h"
#include "MFilterList.h"
#include "MF.h"
#include "MFindSupercutsONOFFThetaLoop.h"
#include "MRanForestFill.h"

#include "MArgs.h"
#include "MArray.h"
#include "MDirIter.h"

#include "MStatusDisplay.h"

#include "MSequence.h"
#include "MJStar.h"
#include "MH3.h"
#include "MRFEnergyEst.h"

#include "MRanForestCalc.h"
#include "MTaskInteractive.h"
#include "MHillas.h"
//#include "MFHadAlpha.h"

#include "MLog.h"
#include "MLogManip.h"
#include "MPrint.h"

using namespace std;


Int_t PreProcess(MParList *plist);
Int_t Process();

static void StartUpMessage()
{
  gLog << all << endl;

  //                1         2         3         4         5         6
  //       123456789012345678901234567890123456789012345678901234567890
  gLog << "========================================================" << endl;
  gLog << "                   Melibea - MARS V" << MARSVER              << endl;
  gLog << "         MARS -- Supercuts Hillas Image Treatment"        << endl;
  gLog << "               Compiled on <" << __DATE__ << ">"          << endl;
  gLog << "                  Using ROOT v" << ROOTVER                << endl;
  gLog << "========================================================" << endl;
  
  gLog << endl;
}

static void Usage()
{
    //                1         2         3         4         5         6         7         8
    //       12345678901234567890123456789012345678901234567890123456789012345678901234567890
    gLog << all << endl;
    gLog << "Sorry the usage is:" << endl;


}

int main(int argc, char **argv)
{
  StartUpMessage();
  
  //
  // Evaluate arguments
  //
  MArgs arg(argc, argv, kTRUE);
  
  if (arg.HasOnly("-V") || arg.HasOnly("--version"))
    return 0;
  
  if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
    {
      Usage();
      return -1;
    }
  
  gLog.Setup(arg);

  const Bool_t  kOptimize   = arg.HasOnlyAndRemove("--optimize");
  const Bool_t  kSupercuts  = arg.HasOnlyAndRemove("--supercuts");
  const Bool_t  kOverwrite  = arg.HasOnlyAndRemove("-f");
  const Bool_t  kPrintSeq   = arg.HasOnlyAndRemove("--print-seq");
  const Bool_t  kBatch      = arg.HasOnlyAndRemove("-b");
  const Bool_t  kPlot       = arg.HasOnlyAndRemove("-p");
  const Bool_t  kMc         = arg.HasOnlyAndRemove("--mc");
  const Bool_t  kUseSeq     = !arg.HasOnlyAndRemove("--nouseseq");

  // Variable that decides whether the optimized cuts are used 
  // in the test sample.
  const Bool_t TestParams   = !arg.HasOnlyAndRemove("--opttest");  
  const Bool_t CombineCosThetaBinsForTrainSample = arg.HasOnlyAndRemove("--combcosthtrain"); 
  const Bool_t CombineCosThetaBinsForTestSample = arg.HasOnlyAndRemove("--combcosthtest"); 
  const Double_t whichfractiontrain = arg.GetFloatAndRemove("--ontrainfrac=",0.5);
  const Double_t whichfractiontest = arg.GetFloatAndRemove("--ontestfrac=",0.5);
  const Double_t whichfractiontrainOFF = arg.GetFloatAndRemove("--offtrainfrac=",0.5);
  const Double_t whichfractiontestOFF = arg.GetFloatAndRemove("--offtestfrac=",0.5);
  const Double_t gammaeff = arg.GetFloatAndRemove("--geff=",0.6);
  const Double_t alphasig = arg.GetFloatAndRemove("--alphasig=",6);
  const Double_t alphabkgmin = arg.GetFloatAndRemove("--alphabgmin=",20);
  const Double_t alphabkgmax = arg.GetFloatAndRemove("--alphabgmax=",80);
  const Double_t SizeLow = arg.GetFloatAndRemove("--sizelow=",60);
  const Double_t SizeUp = arg.GetFloatAndRemove("--sizeup=",1000000);
  const Double_t LeakageMax = arg.GetFloatAndRemove("--leakagemax=",0.25);
  const Double_t DistMax = arg.GetFloatAndRemove("--distmax=",1.0);
  const Double_t DistMin = arg.GetFloatAndRemove("--distmin=",0.2);
  const Double_t fThetaMax = arg.GetFloatAndRemove("--zdmax=",30.);
  const Double_t fThetaMin = arg.GetFloatAndRemove("--zdmin=",0.);
  const Int_t NAlphaBins = arg.GetIntAndRemove("--nalphabin=",30);
  const Double_t AlphaBinLow = arg.GetFloatAndRemove("--alphabinlow=",0);
  const Double_t AlphaBinUp = arg.GetFloatAndRemove("--alphabinup=",90);
  const Bool_t NormFactorFromAlphaBkg = !arg.HasOnlyAndRemove("--normbgmeth");
  const Bool_t UseFittedQuantities = !arg.HasOnlyAndRemove("--usefittedq");
  const Bool_t NotUseTheta = !arg.HasOnlyAndRemove("--notuseth");
  const Bool_t UseStaticCuts = !arg.HasOnlyAndRemove("--staticcuts");
  const Bool_t ReadInitParamsFromAsciiFile = !arg.HasOnlyAndRemove("--ascii");
  const Bool_t kPrintFiles = !arg.HasOnlyAndRemove("--printseq");
  const Int_t NInitSCPar = arg.GetIntAndRemove("--ninitscpar=",104);
  const Int_t degree    = arg.GetIntAndRemove("--degree=",2);
  const TString InitSCParamAsciiFile = arg.GetStringAndRemove("--inscparascii=","mtemp/mmpi/asciifiles/StartingValuesForSmallSizes1.txt");
  const TString ONDataFilename = arg.GetStringAndRemove("--inondataopt=","./*CrabNebula*.root");
  const TString OFFDataFilename = arg.GetStringAndRemove("--inoffdataopt=","./*OffCrab*.root");
  const TString PathForFiles = arg.GetStringAndRemove("--outopt=","./");

  const TString kInpath     = arg.GetStringAndRemove("--ind=", "");
  const TString kOutpath    = arg.GetStringAndRemove("--out=", "");
  const TString kSCfile     = arg.GetStringAndRemove("--scf=", "");
  const TString kSource     = arg.GetStringAndRemove("--src=", "");
  const TString kDate       = arg.GetStringAndRemove("--date=", "");
  const TString kRFTreeFile = arg.GetStringAndRemove("--rftree=", "");
  const Double_t kAlphaMin  = arg.GetFloatAndRemove("--alphamin=",-1);
  const Double_t kAlphaMax  = arg.GetFloatAndRemove("--alphamax=",-1);
  const Double_t kAlphaSig  = arg.GetFloatAndRemove("--alphasignal=",-1);
  const Int_t    kDegree    = arg.GetIntAndRemove("--poldeg=",-1);
  const Int_t    kNumTrees  = arg.GetIntAndRemove("--ntrees=",100);


  if (arg.GetNumOptions()>0)
    {
      gLog << warn << "WARNING - Unknown commandline options..." << endl;
      arg.Print("options");
      gLog << endl;
      return -1;
    }
  
  if (arg.GetNumArguments()!=1)
    {
      Usage();
      return -1;
    }
  //
  // Setup sequence file and check for its existance
  //
  const TString kSequence = arg.GetArgumentStr(0);
  
  if (gSystem->AccessPathName(kSequence, kFileExists))
    {
      gLog << err << "Sorry, sequence file '" << kSequence << "' doesn't exist." << endl;
      return -1;
    }
  
  //
  // Setup sequence and check its validity
  //
  MSequence seq(kSequence);
  if (kPrintSeq)
    {
      gLog << all;
      gLog.Separator(kSequence);
      seq.Print();
      gLog << endl;
    }
  if (!seq.IsValid())
    {
      gLog << err << "Sequence read but not valid!" << endl << endl;
      return -1;
    }
  
  //
  // Process print options
  //

  MDirIter iter;

  if(kUseSeq) {
    gLog << inf;
    gLog << "Calculate image parameters from sequence ";
    gLog << seq.GetName() << endl;
    gLog << endl;


    const Int_t n0 = seq.SetupDatRuns(iter, kInpath, "I");
    const Int_t n1 = seq.GetNumDatRuns();

    if (n0==0)
    {
        gLog << err << "ERROR - No input files of sequence found!" << endl;
        return kFALSE;
    }
    if (n0!=n1)
    {
        gLog << err << "ERROR - Number of files found (" << n0 << ") doesn't match number of files in sequence (" << n1 << ")" << endl;
        return kFALSE;
    }

   }



    //
    // Initialize root
    //
    MArray::Class()->IgnoreTObjectStreamer();
    MParContainer::Class()->IgnoreTObjectStreamer();

    TApplication app("Star", &argc, argv);
    if (!gROOT->IsBatch() && !gClient || gROOT->IsBatch() && !kBatch)
    {
        gLog << err << "Bombing... maybe your DISPLAY variable is not set correctly!" << endl;
        return 1;
    }
    

    // **************************************************** //
    //           Setup flags for optimization               //
    // **************************************************** //
    
    
    // Name of the root file where alpha distributions, TTree objects
    // with info about the events and cuts applied and  info support histograms 
    // will be stored. 
    // Write only the name of the file. The Path 
    // is the one defined previously
    
    TString RootFilename = kOutpath + "RootFileDynCuts.root";
    TString ONFiles = ONDataFilename + "*_I_*.root";
    TString OFFFiles = OFFDataFilename + "*_I_*.root";
    // Vector containing the theta bins in which data (ON/OFF train/test)
    // will be divided. Actually this vector contains the cosinus of 
    // these theta bins. The dimension of the vector is N+1, where 
    // N is the number of theta bins intended. The first component of the 
    // vector is the low bin edge of the first bin, and the last 
    // vector component the upper bin edge of the last bin.
    
    TArrayD CosThetaRangeVector(2);
    CosThetaRangeVector[1] = 0.866;  //30
    CosThetaRangeVector[0] = 0.5;  // 60
    
    // Object of MCT1FindSupercutsONOFFThetaLoop created, data that was specified 
    // above is introduced and ... and the party starts.
    
    MFindSupercutsONOFFThetaLoop FindSupercuts("MFindSupercutsONOFFThetaLoop", 
					       "Optimizer for the supercuts");
    
    
    // ****************************************************** //
    //        Calculate and/or apply dynamical cuts           //
    // ****************************************************** //
    
    // Hillas data file
    //TString datafile = kInpath + "*_I_CrabNebula_E.root";
    
    TString datafile;
    if (kUseSeq)
       datafile = kInpath + kDate + "*_I_" + kSource + "*.root";
    else
       datafile = kInpath + "*" + kSource + "*.root";

     MReadMarsFile read("Events");
     read.DisableAutoScheme(); 
     if (kUseSeq)
        read.AddFiles(iter);
     else
        read.AddFile(datafile);

     cout << " datafile is : " << datafile << endl;

     MReadReports readreal;
     readreal.AddTree("Events", "MTime.", kTRUE);
     readreal.AddTree("Drive");
     readreal.AddTree("EffectiveOnTime");  
     if (kUseSeq) 
	readreal.AddFiles(iter);        
     else
	readreal.AddFile(datafile); 

//     readreal.AddToBranchList("MReportDrive.*");
//     readreal.AddToBranchList("MRawEvtHeader.*");
//     readreal.AddToBranchList("MEffectiveOnTime.*");  

 
    MGeomApply geomapl;
    
    MApplySupercuts appcuts;
    appcuts.SetSCFilename(kSCfile);
    
    MFilterList flist;
    MF hadrf("MHadronness.fHadronness<0.5");
    MF sizef("MHillas.fSize>300");
    MF leakmax("MNewImagePar.fLeakage1<0.25");
    MF distmax("MHillasSrc.fDist<300");
    MF distmin("MHillasSrc.fDist>60");
    MF widthmin("MHillas.fWidth>15");
    MF widthmax("MHillas.fWidth<120");
    
    TString ThetaCutMinString ("MPointingPos.fZd>"); // new! // for magic
    ThetaCutMinString += fThetaMin;
    MContinue ThetaCutMin(ThetaCutMinString);
    ThetaCutMin.SetInverted();

    TString ThetaCutMaxString ("MPointingPos.fZd<"); // new! // for magic
    ThetaCutMaxString += fThetaMax;
    MContinue ThetaCutMax(ThetaCutMaxString);
    ThetaCutMax.SetInverted();

    MContinue cont("(MHillasSrc.fDist<300) && (MHillasSrc.fDist>60) && (MNewImagePar.fLeakage1<0.25) && (MHillas.fSize>60.)");
    cont.SetInverted();


    MRFEnergyEst rfs;
    //rfs.SetFile("/datm3/data/RFEnergyEst/EForestsBeforeCutsLZA.root");
    //rfs.SetFile("/datm3/data/RFEnergyEst/EForests19990101_10003_I_MCGamTrainHZA_E_10_5.root");
    rfs.SetFile("/home/pcmagic16/mazin/mars/files/EForestsHZA.root");

    
    //TString outname = Form("%s{s/_I_/_Q_}", kOutpath.Data());
    MWriteRootFile &write=*new MWriteRootFile(2,  Form("%s{s/_I_/_Q_}", kOutpath.Data()), "RECREATE");
    // Write the Events tree
    write.AddContainer("MHillas",       "Events");
    write.AddContainer("MHillasExt",    "Events");
    write.AddContainer("MHillasSrc",    "Events");
    write.AddContainer("MImagePar",     "Events");
    write.AddContainer("MNewImagePar",  "Events");
    write.AddContainer("MRawEvtHeader", "Events", kFALSE);
    write.AddContainer("MPointingPos",  "Events");
    write.AddContainer("MHadronness","Events", kFALSE);
    write.AddContainer("MEnergyEst","Events", kFALSE);

    if(kMc)
    {
        // Monte Carlo
        write.AddContainer("MMcEvt",              "Events");
        write.AddContainer("MMcTrig",             "Events");
        // Monte Carlo Headers
        write.AddContainer("MMcTrigHeader",       "RunHeaders");
        write.AddContainer("MMcConfigRunHeader",  "RunHeaders");
        write.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
        write.AddContainer("MMcRunHeader",  "RunHeaders", kFALSE);
    }
    else
    {
        write.AddContainer("MTime",               "Events");
        // Run Header
        write.AddContainer("MRawRunHeader",       "RunHeaders");
        write.AddContainer("MBadPixelsCam",       "RunHeaders", kFALSE);
        write.AddContainer("MGeomCam",            "RunHeaders", kFALSE);
        //write.AddContainer("MObservatory", "RunHeaders");
        // Drive
        //write.AddContainer("MSrcPosCam",   "Drive");
        write.AddContainer("MReportDrive",        "Drive", kFALSE);
        write.AddContainer("MTimeDrive",          "Drive", kFALSE);
        // Effective On Time
        write.AddContainer("MEffectiveOnTime",     "EffectiveOnTime", kFALSE);
        write.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime", kFALSE);
    }


/*
    write.AddContainer("MHillas","Events");
    write.AddContainer("MHillasExt","Events");
    write.AddContainer("MHillasSrc","Events");
    write.AddContainer("MImagePar","Events");
    write.AddContainer("MNewImagePar","Events");
    write.AddContainer("MPointingPos","Events");
    write.AddContainer("MRawEvtHeader","Events");
    write.AddContainer("MTime","Events", kFALSE);
    write.AddContainer("MHadronness","Events");
    write.AddContainer("MEnergyEst","Events");
    // Write the RunHeader tree
    write.AddContainer("MRawRunHeader","RunHeaders", kFALSE);
    write.AddContainer("MBadPixelsCam","RunHeaders",kFALSE);
    write.AddContainer("MGeomCam","RunHeaders",kFALSE);
    // Write the Drive tree
    write.AddContainer("MReportDrive","Drive",kFALSE);
    write.AddContainer("MTimeDrive","Drive",kFALSE);
    // Write the EffectiveOnTime tree
    write.AddContainer("MEffectiveOnTime","EffectiveOnTime",kFALSE);
    write.AddContainer("MTimeEffectiveOnTime","EffectiveOnTime",kFALSE);
    

   write.AddContainer("MMcConfigRunHeader",  "RunHeaders", kFALSE);
   write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
   write.AddContainer("MMcFadcHeader",  "RunHeaders", kFALSE);
   write.AddContainer("MMcTrigHeader",  "RunHeaders", kFALSE);
   write.AddContainer("MMcEvt",        "Events", kFALSE); 
   write.AddContainer("MMcTrig",       "Events", kFALSE);
*/
    cout << datafile << endl;
    MParList  plist;
    MTaskList tlist;
    plist.AddToList(&tlist);
    
    MProgressBar bar;

    switch(kSupercuts)
                     {
    case kTRUE:
    switch(kOptimize){
    case kTRUE:
       //cout << ONFiles << endl;
      FindSupercuts.SetPathForFiles(PathForFiles);
      FindSupercuts.SetDataONOFFRootFilenames(ONFiles, OFFFiles);
      FindSupercuts.SetFractionTrainTestOnOffEvents(whichfractiontrain, 
						    whichfractiontest, 
						    whichfractiontrainOFF, 
						    whichfractiontestOFF);
      FindSupercuts.SetGammaEfficiency(gammaeff);
      FindSupercuts.SetAlphaSig(alphasig);
      // Bkg alpha region is set 
      FindSupercuts.SetAlphaBkgMin(alphabkgmin);
      FindSupercuts.SetAlphaBkgMax(alphabkgmax);
      // alpha bkg and signal region set in object FindSupercuts
      // are re-checked in order to be sure that make sense
      FindSupercuts.CheckAlphaSigBkg();
      // Degree of the polynomials used to fit the ON OFF data is set
      FindSupercuts.SetDegreeON(degree);
      FindSupercuts.SetDegreeOFF(degree);
      // binning for alpha plots is defined
      FindSupercuts.SetAlphaPlotBinining(NAlphaBins, AlphaBinLow, AlphaBinUp);
      // Size range is defined
      FindSupercuts.SetSizeRange(SizeLow, SizeUp);
      FindSupercuts.SetFilters(LeakageMax, DistMax, DistMin);
      FindSupercuts.SetNormFactorFromAlphaBkg(NormFactorFromAlphaBkg);
      FindSupercuts.SetUseFittedQuantities(UseFittedQuantities);
      FindSupercuts.SetVariableUseStaticCuts(UseStaticCuts);
      FindSupercuts.SetVariableNotUseTheta(NotUseTheta);
      FindSupercuts.SetReadMatricesFromFile(kFALSE);// normal case kFALSE
      FindSupercuts.SetTrainParameters(kTRUE);
      FindSupercuts.SetSkipOptimization(kFALSE);
      
      FindSupercuts.SetUseInitialSCParams(kTRUE);
      FindSupercuts.SetTestParameters(TestParams);
      FindSupercuts.SetHadronnessName("MHadSC");
      FindSupercuts.SetHadronnessNameOFF("MHadOFFSC");
      FindSupercuts.SetAlphaDistributionsRootFilename (RootFilename);
      FindSupercuts.SetCosThetaRangeVector (CosThetaRangeVector);

// energy estimation from Thomas Hengstebeck

      
      // Names for all root files (matrices, alpha distributions...) are created 
      FindSupercuts.SetALLNames();
      
      if(ReadInitParamsFromAsciiFile)
	{
	  // Initial SC Parameters and steps are retrieved from Ascii file
	  if(!FindSupercuts.ReadSCParamsFromAsciiFile(InitSCParamAsciiFile,NInitSCPar))
	    {
	      cout << "Initial SC Parameters could not be read from Ascii file "
		   << InitSCParamAsciiFile << endl
		   << "Aborting execution of macro... " << endl;
	      return -1;
	    }
	}
      
      // Finally loop over all theta bins defined is executed
      if (!FindSupercuts.LoopOverThetaRanges())
	{
	  cout << "Function MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges()" << endl
	       << "could not be performed" << endl;
	}
      
      // Nex and Significance are computed vs alphasig
      if (!FindSupercuts.ComputeNexSignificanceVSAlphaSig())
	{
	  cout << "Function MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl
	       << "could not be performed" << endl;
	}
    
      // Several theta bins are combined to produced a single alpha plot (for train and test)
      // with single Nex and significances
      if (CombineCosThetaBinsForTrainSample || CombineCosThetaBinsForTestSample)
	{
	  if(!FindSupercuts.ComputeOverallSignificance(CombineCosThetaBinsForTrainSample, 
						       CombineCosThetaBinsForTestSample))
	    {
	      cout << "Function MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl
		   << "could not be performed" << endl;
	    }
	}
      break;
    case kFALSE:      
      if(kAlphaMin!=-1)
	appcuts.SetBGAlphaMin(kAlphaMin);
      if(kAlphaMax!=-1)
	appcuts.SetBGAlphaMax(kAlphaMax);
      if(kAlphaSig!=-1)
	appcuts.SetAlphaSig(kAlphaSig);
      if(kDegree!=-1)
	appcuts.SetPolDegree(kDegree);
      if(kPlot)
	appcuts.SetPlotON();
      appcuts.SetPlotSizeLow(SizeLow);
      appcuts.SetPlotSizeUp(SizeUp);
      cout << "alphasignal=" << kAlphaSig << endl;
      //flist.AddToList(&hadrf);
      flist.AddToList(&leakmax);
      flist.AddToList(&distmax);
      flist.AddToList(&distmin);
      flist.AddToList(&widthmin);
      flist.AddToList(&widthmax);
      flist.AddToList(&sizef);
      appcuts.SetFilter(&flist);
      //write.SetFilter(&flist);
      
      tlist.AddToList(kMc ? (MTask*)&read : (MTask*)&readreal);
      cout << " is montecarlo  ? " << kMc << endl;

      tlist.AddToList(&geomapl);
      tlist.AddToList(&ThetaCutMin, "Events");
      tlist.AddToList(&ThetaCutMax, "Events");
      tlist.AddToList(&flist, "Events");
      tlist.AddToList(&appcuts, "Events");
      tlist.AddToList(&rfs, "Events");
      tlist.AddToList(&cont, "Events");
      
      //tlist.AddToList(&hsigma);
MPrint print1("MTimeDrive");
MPrint print2("MEffectiveOnTime");
//      tlist.AddToList(&print1, "Drive");
//      tlist.AddToList(&print2, "EffectiveOnTime");

      tlist.AddToList(&write);
      MEvtLoop shitloop;
      shitloop.SetParList(&plist);
      //
      // Start to loop over all events
      //
      MProgressBar bar;
      shitloop.SetProgressBar(&bar);
      
      if (!shitloop.Eventloop())
	return -1;
      
      /*
	
	if (!shitloop.PreProcess())
	return -1;
	
	while (tlist.Process())
	{
	//MHadronness *fH = (MHadronness *)plist.FindObject("MHadronness");
	//cout<< "Hadronness " << fH->GetHadronness() << endl;
	}
	
	if (!shitloop.PostProcess())
	return -1;
      */ 
      tlist.PrintStatistics();   


      break;
    }
    break;
    case kFALSE:
      // PUT THINGS FOR THE RANDOM FOREST


      gLog << endl;
      gLog.Separator("Applying RF to data");
      

      //
      // Create a empty Parameter List and an empty Task List
      // The tasklist is identified in the eventloop by its name
      //
      //MParList  plist;
      
      //MTaskList tlist;
      //plist.AddToList(&tlist);
      //
      // ---------------------------------------------------------------
      // ---------------------------------------------------------------
      // first event loop: the trees of the random forest are read in
      // ---------------------------------------------------------------
      // ---------------------------------------------------------------
      //
      MReadTree readrf("Tree",kRFTreeFile);
      readrf.DisableAutoScheme();
   
      MRanForestFill rffill;
      rffill.SetNumTrees(kNumTrees);
      
      tlist.AddToList(&readrf);
      tlist.AddToList(&rffill);
      
      //
      // Eventloop
      //
      MEvtLoop evtloop;
      evtloop.SetParList(&plist);
      evtloop.SetProgressBar(&bar);
      
      if (!evtloop.Eventloop())
        return 1;


      // ---------------------------------------------------------------
      //
      // Apply RF to data
      // ---------------------------------------------------------------
     // ---------------------------------------------------------------
      MTaskList tlist3;
      plist.Replace(&tlist3);

      // Calc RF
      MRanForestCalc calc;
      
      // Fill Hist (Only for MC file with Gammas & Hadrons mixed)
      //MFillH fillh3a("MHHadronness");
      //MFillH fillh3b("MHRanForest");


    // My intercative task to apply cuts
    //MTaskInteractive  mytask;
    //mytask.SetPreProcess(PreProcess);
    //mytask.SetProcess(Process);

      //
      // Hadronness cut
      //
//      MFHadAlpha cut;
//      cut.SetCutType(MFHadAlpha::kHadCut );
//      cut.SetInverted();
//      MContinue contHad(&cut);




    // Add to TaskList
    tlist3.AddToList(kMc ? (MTask*)&read : (MTask*)&readreal);
    //tlist3.AddToList(&SizeCut);
    tlist3.AddToList(&calc);
    tlist3.AddToList(&rfs);
    // tlist3.AddToList(&fillh3a);
//      tlist3.AddToList(&fillh3b);
  
    //tlist.AddToList(&mytask);
    //tlist.AddToList(&contHad);


//  tlist3.AddToList(&flist);
  tlist3.AddToList(&ThetaCutMax);
  tlist3.AddToList(&ThetaCutMin);
  tlist3.AddToList(&write);




  //
    // Execute your analysis
    //
    //    MProgressBar bar;
    evtloop.SetProgressBar(&bar);
    if (!evtloop.Eventloop())
        return 1;

    tlist3.PrintStatistics();

    //  plist.FindObject("MHRanForest")->DrawClone();
//      plist.FindObject("MHHadronness")->DrawClone();




// END OF PUTTING THINGS FOR THE RANDOM FOREST
      break;
    }

//abelardo
  ////////////////////////////////////////////////////////////////////
  //
  // Second event loop: simply copy the tree MMcEvtBasic in the tree
  // "OriginalMC" from the input file to the output file:
      delete &write;
      if(kMc)
      //if(0==1)
	{
	  MParList  plist2;
	  MTaskList tlist2;

	  plist2.AddToList(&tlist2);

	  MReadMarsFile read2("OriginalMC",datafile);
//  read2.AddFile(AnalysisFilename);
	  read2.DisableAutoScheme();

	  tlist2.AddToList(&read2);

	  MWriteRootFile writeOrig(2, Form("%s{s/_I_/_Q_}", kOutpath.Data()),"UPDATE");
	  writeOrig.AddContainer("MMcEvtBasic", "OriginalMC", kFALSE);

	  tlist2.AddToList(&writeOrig);

	  MEvtLoop evtloop2;
	  bar.SetWindowName("Copying OriginalMC tree...");
	  evtloop2.SetProgressBar(&bar);
	  evtloop2.SetParList(&plist2);

	  if (!evtloop2.Eventloop())
	    return -1;

	  tlist2.PrintStatistics();
        }
	// end ab
    return 0;
}

Int_t OptimizeCuts()
{
  return 0;
}






// ---------------------------------------------------------------------
//
//  Interactive Task
//
//
// Data members of the InteractiveTask
MParList      *fParList;
MTaskList     *fTaskList;
MHadronness   *fHadronness;
MHillas       *fHillas;


Int_t PreProcess(MParList *plist)
{
    // Get parameter and tasklist, see Process
    fParList = plist;
    fTaskList = (MTaskList*)plist->FindObject("MTaskList");

    // Search in the PreProcess for the objects
    fHadronness = (MHadronness*)fParList->FindObject("MHadronness");
    fHillas  = (MHillas*)fParList->FindObject("MHillas");

    return kTRUE;
}

// --------------------------------------------------------------
//
Int_t Process()
{
    Float_t size = fHillas->GetSize()/0.18;  //photons!
    Float_t had  =  fHadronness->GetHadronness();

    if(size<500)
        return kCONTINUE;
    
    if(size>500 && size <1000 && had>0.32)
      return kCONTINUE;
    
    if(size>1000 && size <2000 && had>0.28)
      return kCONTINUE;
    
    if(size>2000 && size <4000 && had>0.12)
      return kCONTINUE;
    
    if(size>4000 && size <8000 && had>0.21)
        return kCONTINUE;

    if(size>8000 && had>0.49)
        return kCONTINUE;

    return kTRUE;
}



