/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Wolfgang Wittek 9/2003 <mailto:wittek@mppmu.mpg.de>
!              David Paneque   11/2003 <mailto:dpaneque@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MFindSupercutsONOFF                                                       //
//                                                                         //
// Class for optimizing the parameters of the supercuts                     //
// Using ON and OFF data                                                   //
//                                                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MFindSupercutsONOFF.h"

#include <math.h>            // fabs 

#include <TFile.h>
#include <TArrayD.h>
#include <TMinuit.h>
#include <TCanvas.h>
#include <TStopwatch.h>
#include <TVirtualFitter.h>

#include "MBinning.h"
#include "MContinue.h"
#include "MSupercuts.h"
#include "MSupercutsCalcONOFF.h"
#include "MDataElement.h"
#include "MDataMember.h"


#include <TPostScript.h>

#include "MEvtLoop.h"
#include "MFCT1SelFinal.h"
#include "MF.h"
#include "MFEventSelector.h"
#include "MFEventSelector2.h"
#include "MFillH.h"
//#include "MGeomCamCT1Daniel.h"
//#include "MGeomCamCT1.h"
#include "MGeomCamMagic.h"
#include "MFRandomSplit.h"
#include "MH3.h"
#include "MHCT1Supercuts.h"
#include "MHFindSignificance.h" // To be removed at some point...
#include "MHFindSignificanceONOFF.h"
#include "MTSupercutsApplied.h"
#include "MHMatrix.h"
#include "MHOnSubtraction.h"
#include "MDataValue.h"
// #include "MDataString.h"

#include "MLog.h"
#include "MLogManip.h"
#include "MMatrixLoop.h"
#include "MMinuitInterface.h"
#include "MParList.h"
#include "MProgressBar.h"
#include "MReadMarsFile.h"
#include "MReadTree.h"
#include "MTaskList.h"


ClassImp(MFindSupercutsONOFF);

using namespace std;



// Function that computes the normalization factor using COUNTED events in alpha histogram 
// for ON data  and alpha histogram for OFF data;
// in alpha region defined by AlphaBkgMin and AlphaBkgMax

// It is defined outside the class MFindSupercutsONOFF so that it can be used 
// in function fcnSupercuts


static Double_t ComputeNormFactorFromAlphaBkg(TH1 *histON, TH1 *histOFF, 
					      Double_t AlphaBkgMin, 
					      Double_t AlphaBkgMax)
{                                                      

    Double_t NormFactor = 0.0;
    Double_t ONEvents = 0.0; 
    Double_t OFFEvents = 0.0;

    const Double_t SmallQuantity = 0.01;
    
    Double_t xlo = 0.0;
    Double_t xup = 0.0;
    Double_t width = 0.0;

    Int_t BinCounterOFF = 0;
    Int_t BinCounterON = 0;
    

    // I make a copy of the histograms so that nothing happens to the 
    // histograms used in argument

    TH1* HistON = (TH1*) histON->Clone();
    TH1* HistOFF = (TH1*) histOFF->Clone();

    if ( !HistON )
    {
	gLog << "ComputeNormFactorFromAlphaBkg; " << endl
	      << "Clone of ON histogram could not be generated" 
	      << endl;
	return 0.0;
    }


    
    if ( !HistOFF )
    {
	gLog << "ComputeNormFactorFromAlphaBkg; " << endl
	      << " Clone of OFF histogram could not be generated" 
	      << endl;
	return 0.0;
    }

  // Calculate the number of OFF events in the Bkg region 
  // defined by AlphaBkgMin and AlphaBkgMax 
  // ___________________________________________________
  

    Int_t nbinsOFF = HistOFF -> GetNbinsX();  
    Double_t binwidthOFF = HistOFF -> GetBinWidth(1);

    for (Int_t i=1; i<=nbinsOFF; i++)
    {
	xlo = HistOFF->GetBinLowEdge(i);
	xup = HistOFF->GetBinLowEdge(i+1);
 
	// bin must be completely contained in the bkg region
	if ( xlo >= (AlphaBkgMin-SmallQuantity)  &&  xup <= (AlphaBkgMax+SmallQuantity)    )
	{
	    width = fabs(xup-xlo);
	    if (fabs(width-binwidthOFF) > SmallQuantity)
	    {
		gLog << "ComputeNormFactorFromAlphaBkg; " 
		      << endl << " HistOFF has variable binning, which is not allowed" 
		      << endl;
		return 0.0;
	    }

	    BinCounterOFF++;

	    OFFEvents  += HistOFF->GetBinContent(i);
	    
	}
    }


    // Calculate the number of ON events in the Bkg region 
  // defined by AlphaBkgMin and AlphaBkgMax 
  // ___________________________________________________
  

    Int_t nbinsON = HistON -> GetNbinsX();  
    Double_t binwidthON = HistON -> GetBinWidth(1);

    for (Int_t i=1; i<=nbinsON; i++)
    {
	xlo = HistON->GetBinLowEdge(i);
	xup = HistON->GetBinLowEdge(i+1);
 
	// bin must be completely contained in the bkg region
	if ( xlo >= (AlphaBkgMin-SmallQuantity)  &&  xup <= (AlphaBkgMax+SmallQuantity)    )
	{
	    width = fabs(xup-xlo);
	    if (fabs(width-binwidthON) > SmallQuantity)
	    {
		gLog << "ComputeNormFactorFromAlphaBkg; " 
		      << endl << " HistON has variable binning, which is not allowed" 
		      << endl;
		return 0.0;
	    }

	    BinCounterON++;
	    ONEvents  += HistON->GetBinContent(i);
	    
	}
    }


  

    // NormFactor is computed

      if (ONEvents < SmallQuantity || OFFEvents < SmallQuantity)
      {
	  gLog << "ComputeNormFactorFromAlphaBkg; " 
		<< endl 
		<< "ONEvents or OFFEvents computed in bkg region are < " 
		<< SmallQuantity << endl;
	  return 0.0;
	  
      }

     
	  
    
    NormFactor = ONEvents/OFFEvents;

    Double_t error = 1/ONEvents + 1/OFFEvents;
    error = TMath::Sqrt(error);
    error = error * NormFactor;

     // tmp info
      gLog << "ComputeNormFactorFromAlphaBkg;" << endl
	   << "ON Events in bkg region = " << ONEvents 
	   << " (" << BinCounterON << " bins)"
	   << " , OFF Events in bkg region = " << OFFEvents 
	   << " (" << BinCounterOFF << " bins)" << endl
	   <<"NormFactor computed from bkg region = " << NormFactor 
	   << " +/- " << error << endl;
      // end temp

    return NormFactor;

}








//------------------------------------------------------------------------
//
// fcnSupercuts 
//
// - calculates the quantity to be minimized (using TMinuit)
//
// - the quantity to be minimized is (-1)*significance of the gamma signal
//   in the alpha distribution (after cuts)
//
// - the parameters to be varied in the minimization are the cut parameters
//   (par)
//
static void fcnSupercuts(Int_t &npar, Double_t *gin, Double_t &f, 
                         Double_t *par, Int_t iflag)
{
    //cout <<  "entry fcnSupercuts" << endl;

    //-------------------------------------------------------------
    // save pointer to the MINUIT object for optimizing the supercuts
    // because it will be overwritten 
    // when fitting the alpha distribution in MHFindSignificance
    TMinuit *savePointer = gMinuit;
    //-------------------------------------------------------------


    MEvtLoop* evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit();

    // Event loops for ON and OFF data are recovered. 

    MEvtLoop* ONDataevtloopfcn = &evtloopfcn[0];
    MEvtLoop* OFFDataevtloopfcn = &evtloopfcn[1];


    // Parameter list from event loops for ON and OFF data are recovered

    MParList *ONDataplistfcn   = (MParList*) ONDataevtloopfcn->GetParList();
    MParList *OFFDataplistfcn   = (MParList*) OFFDataevtloopfcn->GetParList();
    
    // MTaskList *ONDataTasklistfcn   = (MTaskList*) ONDataevtloopfcn->GetTaskList();
    // MTaskList *OFFDataTasklistfcn   = (MTaskList*) OFFDataevtloopfcn->GetTaskList();



    // Container for supercuts is retrieved from ONDataplistfcn. 
    // NOTE: The same supercuts parameter container is used in  OFFDataplistfcn

    MSupercuts *super = (MSupercuts*) ONDataplistfcn->FindObject("MSupercuts");
    if (!super)
    {
        gLog << "fcnSupercuts : MSupercuts object '" << "MSupercuts"
            << "' not found... aborting" << endl;
        return;
    }


    // Normalization factor for train sample (Train ON before cuts/ train OFF before cuts) 
    // is retrieved from the ONDataplistfcn

    Double_t NormalizationFactorTrain = 0.0;
    
    MDataValue* NormFactorContainer = (MDataValue*) ONDataplistfcn->FindObject("NormFactorTrain");
    NormalizationFactorTrain = NormFactorContainer -> GetValue();

    
    gLog << "fcnSupercuts : Normalization factor retrieved from ONDataplistfcn: " << endl
	 << "NormalizationFactorTrain = " << NormalizationFactorTrain << endl;
    

    // Degree of polynomials used to fit the ON and the OFF data are retrieved from 
    // the event loop

    Int_t degreeON;
    Int_t degreeOFF;

    MDataValue* DegreeONContainer = (MDataValue*) ONDataplistfcn->FindObject("DegreeON");
    degreeON = int(DegreeONContainer -> GetValue());

    MDataValue* DegreeOFFContainer = (MDataValue*) OFFDataplistfcn->FindObject("DegreeOFF");
    degreeOFF = int(DegreeOFFContainer -> GetValue());

    gLog << "fcnSupercuts : fDegreeON and fDegreeOFF retrieved from ONDataplistfcn and OFFDataplistfcn: " << endl
         << "fDegreeON, fDegreeOFF: = " << degreeON << " ," << degreeOFF << endl;


    Double_t AlphaSig = 0.0;
    MDataValue* AlphaSigContainer = (MDataValue*) ONDataplistfcn->FindObject("AlphaSigValue");
    AlphaSig = AlphaSigContainer -> GetValue();

    Double_t AlphaBkgMin = 0.0;
    MDataValue* AlphaBkgMinContainer = 
	(MDataValue*) ONDataplistfcn->FindObject("AlphaBkgMinValue");
    AlphaBkgMin = AlphaBkgMinContainer -> GetValue();

    Double_t AlphaBkgMax = 0.0;
    MDataValue* AlphaBkgMaxContainer = 
	(MDataValue*) ONDataplistfcn->FindObject("AlphaBkgMaxValue");
    AlphaBkgMax = AlphaBkgMaxContainer -> GetValue();

    
    gLog << "fcnSupercuts : AlphaSig and AlphaBkgMin-AlphaBkgMax retrieved from ONDataplistfcn: "
	 << endl
	 << "AlphaSig = " << AlphaSig << "; AlphaBkgMin-AlphaBkgMax = " 
	 << AlphaBkgMin << "-" << AlphaBkgMax << endl;



    // Variable fUseFittedQuantities is retrieved from ONDataplistfcn

    Bool_t UseFittedQuantities;
    MDataValue* UseFittedQuantitiesContainer = 
	(MDataValue*) ONDataplistfcn->FindObject("UseFittedQuantitiesValue");
    UseFittedQuantities = bool(UseFittedQuantitiesContainer -> GetValue());
    
    if (UseFittedQuantities)
    {
	gLog << "fcnSupercuts :  UseFittedQuantities variable set to kTRUE" << endl;
    }
    else
    {
	gLog << "fcnSupercuts :  UseFittedQuantities variable set to kFALSE" << endl;
    }
    

    Bool_t UseNormFactorFromAlphaBkg;
    MDataValue* UseNormFactorFromAlphaBkgContainer = 
      (MDataValue*) ONDataplistfcn -> FindObject("UseNormFactorFromAlphaBkgValue");
    UseNormFactorFromAlphaBkg = bool(UseNormFactorFromAlphaBkgContainer -> GetValue());


    if (UseNormFactorFromAlphaBkg)
    {
	gLog << "fcnSupercuts :  UseNormFactorFromAlphaBkg variable set to kTRUE" << endl;
    }
    else
    {
	gLog << "fcnSupercuts :  UseNormFactorFromAlphaBkg variable set to kFALSE" << endl;
    }




    //
    // transfer current parameter values to MSupercuts
    //
    // Attention : npar is the number of variable parameters
    //                  not the total number of parameters
    //
    Double_t fMin, fEdm, fErrdef;
    Int_t     fNpari, fNparx, fIstat;
    gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);

    super->SetParameters(TArrayD(fNparx, par));

    //$$$$$$$$$$$$$$$$$$$$$
    // for testing
    //  TArrayD checkparameters = super->GetParameters();
    //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ",  " 
    //     << fNparx  << endl;
    //gLog << "fcnsupercuts : i, par, checkparameters =" << endl;
    //for (Int_t i=0; i<fNparx; i++)
    //{
    //gLog << i << ",  " << par[i] << ",  " << checkparameters[i] << endl;
    //  }
    //$$$$$$$$$$$$$$$$$$$$$

    //
    // plot alpha for ON data with the current cuts
    //
    if (!ONDataevtloopfcn->Eventloop())
    {
	gLog << "fcnsupercuts  : ONDataevtloopfcn->Eventloop() failed" << endl;
    }
    
    // Somehow (??) I can not use the function MTaskList::PrintStatistics...
    // it produces a segmentation fault...

    //ONDataTasklistfcn->PrintStatistics(0, kFALSE);
    
    MH3* alpha = (MH3*)ONDataplistfcn->FindObject("AlphaFcn", "MH3");
    if (!alpha)
        return;

    TH1 &alphaHist = alpha->GetHist();
    alphaHist.SetName("alpha-fcnSupercuts");


    //
    // plot alpha for OFF data with the current cuts
    //
    if(!OFFDataevtloopfcn->Eventloop())
    {
	gLog << "fcnsupercuts  : OFFDataevtloopfcn->Eventloop() failed" << endl;
    }

    // Somehow (??) I can not use the function MTaskList::PrintStatistics...
    // it produces a segmentation fault...

    //OFFDataTasklistfcn->PrintStatistics(0, kFALSE);

    MH3* alphaOFF = (MH3*) OFFDataplistfcn->FindObject("AlphaOFFFcn", "MH3");
    if (!alphaOFF)
        return;

    TH1 &alphaHistOFF = alphaOFF->GetHist();
    alphaHistOFF.SetName("alphaOFF-fcnSupercuts");





    if (UseNormFactorFromAlphaBkg)
      {
	// Normalization factor computed using alpha bkg region
	Double_t NewNormFactor = 
	    ComputeNormFactorFromAlphaBkg(&alphaHist, &alphaHistOFF, 
					  AlphaBkgMin, AlphaBkgMax);
	
	gLog << "Normalization factor computed from alpha plot (after cuts) " << endl
	      << "using counted number of ON and OFF events in alpha region " << endl
	      << "defined by range " << AlphaBkgMin << "-" << AlphaBkgMax << endl
	      << "Normalization factor = " << NewNormFactor << endl;

	gLog << "Normalization factor used is the one computed in bkg region; " << endl
	      << "i.e. " << NewNormFactor << " instead of " << NormalizationFactorTrain << endl;

	NormalizationFactorTrain = NewNormFactor;


    }







    //-------------------------------------------
    // set Minuit pointer to zero in order not to destroy the TMinuit
    // object for optimizing the supercuts
    gMinuit = NULL;

    //=================================================================
    // fit alpha distribution to get the number of excess events and
    // calculate significance of gamma signal in the alpha plot
  
    const Double_t alphasig = AlphaSig;
    const Double_t alphamin = AlphaBkgMin; 
   // alpha min for bkg region in ON data
    const Double_t alphamax = AlphaBkgMax; // alpha max for bkg region in ON data


    Bool_t drawpoly;
    Bool_t fitgauss;
    Bool_t saveplots;
    

    if (iflag == 3)
    {// Even though minimization finished successfully, I will NOT produce 
	// final plots now... i'll do it later via the function 
	// "MFindSupercutsONOFF::OutputOptimizationOnTrainSample()"

        drawpoly  = kFALSE;
        fitgauss  = kFALSE;
	saveplots = kFALSE;
    }
    else
    {
        drawpoly  = kFALSE;
        fitgauss  = kFALSE;
	saveplots = kFALSE;
	
	
    }


    const Bool_t print = kTRUE;

    MHFindSignificanceONOFF findsig;
    findsig.SetRebin(kTRUE);
    findsig.SetReduceDegree(kFALSE);
    findsig.SetUseFittedQuantities(UseFittedQuantities);
    
    // TPostScript* DummyPs  = new TPostScript("dummy.ps");
    
    TString DummyPs = ("Dummy");
  

    const Bool_t rc = findsig.FindSigmaONOFF(&alphaHist,&alphaHistOFF, 
					     NormalizationFactorTrain, 
					     alphamin, alphamax, 
					     degreeON, degreeOFF,
					     alphasig, drawpoly, fitgauss, 
					     print, saveplots, 
					     DummyPs);
    //DummyPs -> Close();
    //delete DummyPs;

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

    // reset gMinuit to the MINUIT object for optimizing the supercuts 
    gMinuit = savePointer;
    //-------------------------------------------

    if (!rc)
    {
        gLog << "fcnSupercuts : FindSigmaONOFF() failed" << endl;
        f = 1.e10;
        return;
    }

    /*

    // plot some quantities during the optimization
    MHCT1Supercuts *plotsuper = (MHCT1Supercuts*)ONDataplistfcn->FindObject("MHCT1Supercuts");
    if (plotsuper)
        plotsuper->Fill(&findsig);



    */

    //------------------------
    // get significance
    const Double_t significance = findsig.GetSignificance();
    f = significance>0 ? -significance : 0;


    gLog << " ///****************************************** ///" << endl;
    gLog << "Significance (Li&Ma)is now: " << f << endl;
    gLog << " ///****************************************** ///" << endl;

    //------------------------
    // optimize signal/background ratio
    //Double_t ratio = findsig.GetNbg()>0.0 ? 
    //                 findsig.GetNex()/findsig.GetNbg() : 0.0; 
    //f = -ratio;

    //-------------------------------------------
    // final calculations
    //if (iflag == 3)
    //{
    //
    //}    

    //-------------------------------------------------------------
}


// --------------------------------------------------------------------------
//
// Default constructor.
//
MFindSupercutsONOFF::MFindSupercutsONOFF(const char *name, const char *title)
    : fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL)
{
    fName  = name  ? name  : "MFindSupercutsONOFF";
    fTitle = title ? title : "Optimizer of the supercuts";

    //---------------------------
    // camera geometry is needed for conversion mm ==> degree
    //fCam = new MGeomCamCT1Daniel; 
    // fCam = new MGeomCamCT1; 
    fCam = new MGeomCamMagic; 

    // matrices to contain the training/test samples
    fMatrixTrain = new MHMatrix("MatrixTrain");
    fMatrixTest  = new MHMatrix("MatrixTest");
    fMatrixTrainOFF = new MHMatrix("MatrixTrainOFF");
    fMatrixTestOFF  = new MHMatrix("MatrixTestOFF");
    

    // objects of MSupercutsCalcONOFF to which these matrices are attached
    fCalcHadTrain = new MSupercutsCalcONOFF("SupercutsCalcTrain");
    fCalcHadTest  = new MSupercutsCalcONOFF("SupercutsCalcTest");
    fCalcHadTrainOFF = new MSupercutsCalcONOFF("SupercutsCalcTrainOFF");
    fCalcHadTestOFF  = new MSupercutsCalcONOFF("SupercutsCalcTestOFF");

    // Define columns of matrices
    fCalcHadTrain->InitMapping(fMatrixTrain);
    fCalcHadTest->InitMapping(fMatrixTest);
    fCalcHadTrainOFF->InitMapping(fMatrixTrainOFF);
    fCalcHadTestOFF->InitMapping(fMatrixTestOFF);

    // For the time being weights method is not implemented
    //UseWeights = kFALSE; 

    // Normalization factors are initialized to -1
    // Functions determining Number of excess events will not work 
    // with negative norm factors;  
    // This will ensure that Norm factors are computed before running
    // DetExcessONOFF() function.

    fNormFactorTrain = -1.0;
    fNormFactorTest = -1.0;

    // SigmaLiMa and Nex member variables are initialized to 0

    fSigmaLiMaTrain = 0.0;
    fSigmaLiMaTest = 0.0;
    fNexTrain = 0.0; 
    fNexTest = 0.0; 

    
    // Cuts (low and up) in variable MPointingPos.fZd
    // (at some point real theta)
    // The default is not cut, i.e. all values (0-1) are taken
    
    //fThetaMin = 0; // in miliradians  // FIXME: change name
    //fThetaMax = 1570; // in miliradians // FIXME: change name
    
    fThetaMin = 0; // degrees (new in magic)
    fThetaMax = 90; // degrees (new in magic)



    fAlphaSig = 20; // By default, signal is expected in alpha<20 for CT1

    // By default, bkg region is set to alpha range 30-90 for CT1
    fAlphaBkgMin = 30; 
    fAlphaBkgMax = 90;


    // By default, bining for alpha plots is set to 
    
    fNAlphaBins = 54;
    fAlphaBinLow = -12;
    fAlphaBinUp = 96;



    fNormFactorFromAlphaBkg = kTRUE; // By default normalization factor is 
    // computed using bkg region in alpha histograms (after cuts)

//    fActualCosThetaBinCenter = 0;

    fTuneNormFactor = kTRUE; // Norm factors will be corrected using 
    //the total amount of OFF events before cuts and the estimated excess events
    // fNormFactorTrain = fNormFactorTrain - Ngammas/EventsInTrainMatrixOFF


    // use quantities computed from the fits
    // The variable allows the user to NOT use these quantities when there is not 
    // enough statistics and fit not always is possible.
    // Default value is kTRUE
    fUseFittedQuantities = kTRUE;




    // Boolean variable that allows the user to set some limits to the 
    // some of the Minuit parameters. For the time being, only the limits
    // for the parameters which do NOT depend in size, dist and theta are set, 
    // i.e. static limits. The value of this boolean variable is set in the 
    // constructor of the class.
    
    fSetLimitsToSomeMinuitParams = kTRUE;
    
    // Limits for the Minuit parameters. For the time being the values are set in the constructor 
    // of the class. One MUST be very careful to set limits such that the expected final values 
    // (optimized values) are far away from the limits.
    
    // Values in degrees
    

    fMinuitDistUPUpperLimit = 2.0;
    fMinuitDistUPLowerLimit = 0.5;
    fMinuitLengthUPUpperLimit = 0.8;
    fMinuitLengthUPLowerLimit = 0.0;
    fMinuitWidthUPUpperLimit = 0.5;
    fMinuitWidthUPLowerLimit = 0.0;
    fMinuitLeakage1UPUpperLimit = 1.5;
    fMinuitLeakage1UPLowerLimit = -0.5;
    
    fMinuitDistLOWUpperLimit = 1.0;
    fMinuitDistLOWLowerLimit = 0.0;
    fMinuitLengthLOWUpperLimit = 0.5;
    fMinuitLengthLOWLowerLimit = -0.3;
    fMinuitWidthLOWUpperLimit = 0.4;
    fMinuitWidthLOWLowerLimit = -0.3;



    // Boolean variable that controls wether the optimization of the 
    // parameters (MMinuitInterface::CallMinuit(..) in function FindParams(..))
    // takes place or not. kTRUE will skip such optimization.
    // This variable is useful to test the optmized parameters (previously found 
    // and stored in root file) on the TRAIN sample.
    
    fSkipOptimization = kFALSE;

    // Boolean variable that allows the user to write the initial parameters 
    // into the root file that will be used to store the optimum cuts.
    // If fUseInitialSCParams = kTRUE , parameters are written.
    // In this way, the initial SC parameters can be applied on the data (train/test) 
    
    // The initial parameters are ONLY written to the root file if 
    // there is NO SC params optimization, i.e., if variable 
    // fSkipOptimization = kTRUE;

    // The default value is obviously kFALSE.

    fUseInitialSCParams = kTRUE;
    

    // Set wether to use or not hillas dist
    fUseDist = kTRUE;




    fGammaEfficiency = 0.5; // Fraction of gammas that remain after cuts
    // Quantity that will have to be determined with MC, yet for the 
    // time being I set it to 0.5 (standard value)

    fPsFilename = NULL;
    fPsFilename2 = NULL;


    ////////////////////////////////////////////////////
    // TMP

    // There are quite some problems during the data preprocessing. 
    // For the time being, I will add some cuts to the functions 
    // DefineTrainTestMatrixThetaRange and for OFF, so that I can 
    // make a kind of preprocess on my own. This allows me 
    // to make a very silly preprocess with wolfgangs macro, which 
    // might be free of corrupted data, and then I can do on my own.

    fSizeCutLow = 0.1; // To prevent empty events
    fSizeCutUp = 10000000;

    // Angular cuts are converted to mm, which 
    // are the units of the preprocessed data....


    // Angular cuts not yet implemented ...
    fConvMMToDeg = 0.00337034;

    fDistCutLow = 0.4/fConvMMToDeg;
    fDistCutUp = 1.5/fConvMMToDeg;
    
    fLengthCutLow = 0.1/fConvMMToDeg;
    fLengthCutUp = 1/fConvMMToDeg;

    fWidthCutLow = 0.07/fConvMMToDeg;
    fWidthCutUp = 1/fConvMMToDeg;

    // ENDTMP
    ////////////////////////////////////////////////////
    
    
    // Degree of polynomials used to fit the ON OFF data is 
    // set initially to 2.

    fDegreeON = 2;
    fDegreeOFF = 2;


}

// --------------------------------------------------------------------------
//
// Default destructor.
//
MFindSupercutsONOFF::~MFindSupercutsONOFF()
{

    *fLog << "destructor of MFindSupercutsONOFF is called" << endl;

    fPsFilename = NULL;
    fPsFilename2 = NULL;


    delete fCam;
    delete fMatrixTrain;
    delete fMatrixTest;
    delete fCalcHadTrain;
    delete fCalcHadTest;
    delete fMatrixTrainOFF;
    delete fMatrixTestOFF;
    delete fCalcHadTrainOFF;
    delete fCalcHadTestOFF;
    
    *fLog << "destructor of MFindSupercutsONOFF finished successfully" << endl;

}

// --------------------------------------------------------------------------
//
// Function that sets the name of the PostScript file where alpha distributions 
// for the different Theta bins will be stored.
// It also initializes 



void MFindSupercutsONOFF::SetPostScriptFile (TPostScript* PsFile)
{
    fPsFilename = PsFile;

    *fLog << "MFindSupercutsONOFF : Results (alpha distributions with excess and significances) will be stored in PostScript file " 
	  << fPsFilename -> GetName() << endl;

}

void MFindSupercutsONOFF::SetPostScriptFile2 (TPostScript &PsFile)
{
    fPsFilename2 = new TPostScript (PsFile);

    *fLog << "MFindSupercutsONOFF : Results (alpha distributions with excess and significances) will be stored in PostScript file " 
	  << fPsFilename2 -> GetName() << endl;

}

void MFindSupercutsONOFF::SetPsFilenameString (const TString filename)
{
    fPsFilenameString = filename;
}

void MFindSupercutsONOFF::SetSkipOptimization(Bool_t b)
{
    fSkipOptimization = b;
    if (fSkipOptimization)
    {
	*fLog << "MFindSupercutsONOFF :: SetSkipOptimization " << endl
	      << "Variable fSkipOptimization is kTRUE, and therefore "
	      << "the optimization of supercuts is skipped. Hope that's "
	      << "what you want... " << endl;
    }

}


void MFindSupercutsONOFF::SetUseInitialSCParams(Bool_t b)
{
    fUseInitialSCParams = b;
    if (fUseInitialSCParams)
    {
	*fLog << "MFindSupercutsONOFF :: SetUseInitialSCParams " << endl
	      << "Variable fUseInitialSCParams is kTRUE. " << endl;

	if (fSkipOptimization)
	{
	    *fLog << "The Initial SC Parameters will be applied on the selected data."
		  << endl;
	}

	else
	{
	    *fLog << "However, fSkipOptimization = kFALSE, and therefore, the " 
		  << "the supercuts will be optimized. The final cuts that " 
		  << "will be applied to the data will NOT be the initial SC parameters."
		  << endl;

	}



    }
    
}



Bool_t MFindSupercutsONOFF::SetAlphaSig(Double_t alphasig)
{
    // check that alpha is within the limits 0-90
    if (alphasig <= 0 || alphasig > 90)
    {
	*fLog << "MFindSupercutsONOFF ::SetAlphaSig; " 
	      << "value " << alphasig << " is not within the the "
	      << "logical limits of alpha; 0-90" << endl;
	return kFALSE;
    }
    

    fAlphaSig = alphasig;

    return kTRUE;
}

Bool_t MFindSupercutsONOFF::SetAlphaBkgMin(Double_t alpha)
{
    // check that alpha is within the limits 0-90
    if (alpha <= 0 || alpha >= 90)
    {
	*fLog << "MFindSupercutsONOFF ::SetAlphaBkgMin; " 
	      << "value " << alpha << " is not within the the "
	      << "logical limits of alpha; 0-90" << endl;
	return kFALSE;
    }
    

    fAlphaBkgMin = alpha;

    return kTRUE;
}


Bool_t MFindSupercutsONOFF::SetAlphaBkgMax(Double_t alpha)
{
    // check that alpha is within the limits 0-90
    if (alpha <= 0 || alpha > 90.001)
    {
	*fLog << "MFindSupercutsONOFF ::SetAlphaBkgMax; " 
	      << "value " << alpha << " is not within the the "
	      << "logical limits of alpha; 0-90" << endl;
	return kFALSE;
    }
    

    fAlphaBkgMax = alpha;

    return kTRUE;
}


// Function that checks that the values of the member data 
// fAlphaSig, fAlphaBkgMin and fAlphaBkgMax make sense
// (ie, fAlphaSig < fAlphaBkgMin < fAlphaBkgMax)

Bool_t MFindSupercutsONOFF::CheckAlphaSigBkg()
{

    if (fAlphaSig > fAlphaBkgMin)
    {
	*fLog << "MFindSupercutsONOFF ::CheckAlphaSigBkg(); " << endl
	      << "fAlphaSig > fAlphaBkgMin, which should not occur..." << endl
	      << "fAlphaSig = " << fAlphaSig << ", fAlphaBkgMin = " << fAlphaBkgMin
	      << endl;

	return kFALSE;
    }
    
    if (fAlphaBkgMax < fAlphaBkgMin)
    {
	*fLog << "MFindSupercutsONOFF ::CheckAlphaSigBkg(); " << endl
	      << "fAlphaBkgMin > fAlphaBkgMax, which should not occur..." << endl
	      << "fAlphaBkgMin = " << fAlphaBkgMin << ", fAlphaBkgMax = " << fAlphaBkgMax
	      << endl;

	return kFALSE;
    }

    return kTRUE;

}


/*
// Function that computes the normalization factor using COUNTED events ON and OFF
// in alpha region defined by fAlphaBkgMin and fAlphaBkgMax


Double_t MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg(TH1 *histON, TH1 *histOFF)
{

    Double_t NormFactor = 0.0;
    Double_t ONEvents = 0.0; 
    Double_t OFFEvents = 0.0;

    const Double_t SmallQuantity = 0.01;
    
    Double_t xlo = 0.0;
    Double_t xup = 0.0;
    Double_t width = 0.0;

    Int_t BinCounterOFF = 0;
    Int_t BinCounterON = 0;
    

    // I make a copy of the histograms so that nothing happens to the 
    // histograms used in argument

    TH1* HistON = (TH1*) histON->Clone();
    TH1* HistOFF = (TH1*) histOFF->Clone();

    if ( !HistON )
    {
	*fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl
	      << "Clone of ON histogram could not be generated" 
	      << endl;
	return 0.0;
    }


    
    if ( !HistOFF )
    {
	*fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl
	      << " Clone of OFF histogram could not be generated" 
	      << endl;
	return 0.0;
    }

  // Calculate the number of OFF events in the Bkg region 
  // defined by fAlphaBkgMin and fAlphaBkgMax 
  // ___________________________________________________
  

    Int_t nbinsOFF = HistOFF -> GetNbinsX();  
    Double_t binwidthOFF = HistOFF -> GetBinWidth(1);

    for (Int_t i=1; i<=nbinsOFF; i++)
    {
	xlo = HistOFF->GetBinLowEdge(i);
	xup = HistOFF->GetBinLowEdge(i+1);
 
	// bin must be completely contained in the bkg region
	if ( xlo >= (fAlphaBkgMin-SmallQuantity)  &&  xup <= (fAlphaBkgMax+SmallQuantity)    )
	{
	    width = fabs(xup-xlo);
	    if (fabs(width-binwidthOFF) > SmallQuantity)
	    {
		*fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " 
		      << endl << " HistOFF has variable binning, which is not allowed" 
		      << endl;
		return 0.0;
	    }

	    BinCounterOFF++;

	    OFFEvents  += HistOFF->GetBinContent(i);
	    
	}
    }


    // Calculate the number of ON events in the Bkg region 
  // defined by fAlphaBkgMin and fAlphaBkgMax 
  // ___________________________________________________
  

    Int_t nbinsON = HistON -> GetNbinsX();  
    Double_t binwidthON = HistON -> GetBinWidth(1);

    for (Int_t i=1; i<=nbinsON; i++)
    {
	xlo = HistON->GetBinLowEdge(i);
	xup = HistON->GetBinLowEdge(i+1);
 
	// bin must be completely contained in the bkg region
	if ( xlo >= (fAlphaBkgMin-SmallQuantity)  &&  xup <= (fAlphaBkgMax+SmallQuantity)    )
	{
	    width = fabs(xup-xlo);
	    if (fabs(width-binwidthON) > SmallQuantity)
	    {
		*fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " 
		      << endl << " HistON has variable binning, which is not allowed" 
		      << endl;
		return 0.0;
	    }

	    BinCounterON++;
	    ONEvents  += HistON->GetBinContent(i);
	    
	}
    }


  

    // NormFactor is computed

      if (ONEvents < SmallQuantity || OFFEvents < SmallQuantity)
      {
	  *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " 
		<< endl 
		<< "ONEvents or OFFEvents computed in bkg region are < " 
		<< SmallQuantity << endl;
	  return 0.0;
	  
      }

     
	  
    
    NormFactor = ONEvents/OFFEvents;

    Double_t error = 1/ONEvents + 1/OFFEvents;
    error = TMath::Sqrt(error);
    error = error * NormFactor;

     // tmp info
      cout << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg;" << endl
	   << "ON Events in bkg region = " << ONEvents 
	   << " (" << BinCounterON << " bins)"
	   << " , OFF Events in bkg region = " << OFFEvents 
	   << " (" << BinCounterOFF << " bins)" << endl
	   <<"NormFactor computed from bkg region = " << NormFactor 
	   << " +/- " << error << endl;
      // end temp

    return NormFactor;

}


*/


// Function that set the values of fThetaMin and fThetaMax and also 
// fThetaRangeString (with value in miliradians); 
// data members that are used basically to print/plot
// information. 

Bool_t MFindSupercutsONOFF::SetThetaRange(Double_t ThetaMin, Double_t ThetaMax)
{

// Check  that values are reasonable... well ... i guess this was done 
// in previous functions...

  //    fThetaMin = int(ThetaMin*1000.0);
  //  fThetaMax = int(ThetaMax*1000.0);


  // in new magic theta is given in deg
  // 0.5 is added so that the rounding to integer is correct
  fThetaMin = int(ThetaMin + 0.5); 
  fThetaMax = int(ThetaMax + 0.5); 

    fThetaRangeString = ("ThetaRange");
    fThetaRangeString += (fThetaMin);
    fThetaRangeString += ("_");
    fThetaRangeString += (fThetaMax);
    // fThetaRangeString += ("mRad");
    fThetaRangeString += ("_Degrees"); // new in magic


    return kTRUE;

}

// Function that sets Size range
Bool_t MFindSupercutsONOFF::SetSizeRange(Double_t SizeMin, Double_t SizeMax)
{


    fSizeCutLow = SizeMin;
    fSizeCutUp = SizeMax;

    *fLog << "MFindSupercutsONOFF::SetSizeRange" << endl
	  << "Data matrices will be filled with events whose MHillas.fSize " << endl
	  << "is in the range " 
	  << fSizeCutLow <<"-"<<fSizeCutUp << endl;



    return kTRUE;
}


Bool_t MFindSupercutsONOFF::SetFilters(Double_t LeakageMax, Double_t DistMax, Double_t DistMin)
{

    fDistCutLow = DistMin/fConvMMToDeg;
    fDistCutUp  = DistMax/fConvMMToDeg;
    fLeakageMax = LeakageMax;
    
    *fLog << "MFindSupercutsONOFF::SetSizeRange" << endl
	  << "Data matrices will be filled with events whose MHillasSrc.fDist " << endl
	  << "is in the range " 
	  << fDistCutLow <<"-"<< fDistCutUp << " degrees" << endl
	  << "and fLeakage " << "< " << fLeakageMax << endl;

    return kTRUE;
}



// Function that sets the names of all parameter containers 
// used to store the supercuts applied to ON/OFF Train/Test samples

void MFindSupercutsONOFF::SetSupercutsAppliedTreeNames()
{

    char* sc = {"SupercutsApplied"};

    fTrainONSupercutsAppliedName = (sc);  
    fTrainONSupercutsAppliedName += ("TrainON");

    fTrainOFFSupercutsAppliedName = (sc);       
    fTrainOFFSupercutsAppliedName += ("TrainOFF");

    fTestONSupercutsAppliedName = (sc);
    fTestONSupercutsAppliedName += ("TestON");

    fTestOFFSupercutsAppliedName = (sc); 
    fTestOFFSupercutsAppliedName += ("TestOFF");
    

    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::SetSupercutsAppliedTreeNames; "
	      << " fThetaRangeString is NOT defined..." << endl;

	fTrainONSupercutsAppliedName += ("ThetaRangeStringUndefined");
	fTrainOFFSupercutsAppliedName += ("ThetaRangeStringUndefined");
	fTestONSupercutsAppliedName += ("ThetaRangeStringUndefined");
	fTestOFFSupercutsAppliedName += ("ThetaRangeStringUndefined");
    }
    else
    {
	fTrainONSupercutsAppliedName += (fThetaRangeString);
	fTrainOFFSupercutsAppliedName += (fThetaRangeString);
	fTestONSupercutsAppliedName += (fThetaRangeString);
	fTestOFFSupercutsAppliedName += (fThetaRangeString);
    }
     
    // Info about names

    *fLog << "MFindSupercutsONOFF::SetSupercutsAppliedTreeNames; " << endl
	  << "Names of the MTSupercutsApplied Trees for Train (ON/OFF) " 
	  << " and Test (ON/OFF) samples are the following ones: " << endl
	  << fTrainONSupercutsAppliedName << ", " 
	  << fTrainOFFSupercutsAppliedName << ", "
	  << fTestONSupercutsAppliedName << ", "
	  << fTestOFFSupercutsAppliedName << endl;
  
}


void MFindSupercutsONOFF::SetUseOrigDistribution(Bool_t b)
{
  fUseOrigDistribution = b;

  if (fUseOrigDistribution == kTRUE)
  {
 
      *fLog << "MFindSupercutsONOFF : when defining training and test matrices use the original distribution" 
	    << endl;
  }

  else 
  {
      *fLog << "MFindSupercutsONOFF : when defining training and test matrices, events will be selected according to distribution given by a MH3 object" 
	    << endl;
  }
}



// --------------------------------------------------------------------------
//
// Define the matrix 'fMatrixTrain' for the training sample
//
// alltogether 'howmanytrain' events are read from file 'nametrain';
// the events are selected according to a target distribution 'hreftrain'
//
//
Bool_t MFindSupercutsONOFF::DefineTrainMatrix(const TString &nametrain, MH3 &hreftrain,
						 const Int_t howmanytrain, 
						 const TString &filetrain)
{
    if (nametrain.IsNull() || howmanytrain <= 0)
        return kFALSE;

    *fLog << "=============================================" << endl;
    *fLog << "fill training matrix from file '" << nametrain 
          << "',   select " << howmanytrain 
          << " events " << endl;

    if (!fUseOrigDistribution)
    {
      *fLog << "     according to a distribution given by the MH3 object '"
            << hreftrain.GetName() << "'" << endl;
    }
    else
    {
      *fLog << "     randomly" << endl;
    }
    


    MParList  plist;
    MTaskList tlist;

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


    MFEventSelector2 seltrain(hreftrain);
    seltrain.SetNumMax(howmanytrain);
    seltrain.SetName("selectTrain");
    if (fUseOrigDistribution)
    {
	*fLog << "MFindSupercutsONFF; MFEventSelector2::SetUseOrigDistribution(Bool)"
	      << " DOES NOT EXIST NOW..." << endl;
	// seltrain.SetUseOrigDistribution(kTRUE);
    }


    MFillH filltrain(fMatrixTrain);
    filltrain.SetFilter(&seltrain);
    filltrain.SetName("fillMatrixTrain");

    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTrain);

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

    tlist.AddToList(&read);
    tlist.AddToList(&seltrain);
    tlist.AddToList(&filltrain);

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

    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTrain");
    evtloop.SetProgressBar(&bar);

    if (!evtloop.Eventloop())
      return kFALSE;

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTrain->Print("SizeCols");
    Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
    if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
    {
      *fLog << "MFindSupercutsONOFF::DefineTrainMatrix; no.of generated events ("
	    << howmanygenerated 
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }

    *fLog << "training matrix was filled" << endl;
    *fLog << "=============================================" << endl;

    //--------------------------
    // write out training matrix

    if (filetrain != "")
    {
      TFile filetr(filetrain, "RECREATE");
      fMatrixTrain->Write();
      filetr.Close();

      *fLog << "MFindSupercuts::DefineTrainMatrix; Training matrix was written onto file '"
            << filetrain << "'" << endl;
    }


    return kTRUE;
}


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

 

// --------------------------------------------------------------------------
//
// Define the matrix for the test sample
//
// alltogether 'howmanytest' events are read from file 'nametest'
//
// the events are selected according to a target distribution 'hreftest'
//
//
Bool_t MFindSupercutsONOFF::DefineTestMatrix(const TString &nametest, MH3 &hreftest,
						const Int_t howmanytest, const TString &filetest)
{
    if (nametest.IsNull() || howmanytest<=0)
        return kFALSE;

    *fLog << "=============================================" << endl;
    *fLog << "fill test matrix from file '" << nametest 
          << "',   select " << howmanytest << " events " << endl;


    if (!fUseOrigDistribution)
    {
      *fLog << "     according to a distribution given by the MH3 object '"
            << hreftest.GetName() << "'" << endl;
    }
    else
    {
      *fLog << "     randomly" << endl;
    }


    MParList  plist;
    MTaskList tlist;

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

    MFEventSelector2 seltest(hreftest);
    seltest.SetNumMax(howmanytest);
    seltest.SetName("selectTest");
    if (fUseOrigDistribution)
    {
	*fLog << "MFindSupercutsONFF; MFEventSelector2::SetUseOrigDistribution(Bool)"
	      << " DOES NOT EXIST NOW..." << endl;
	// seltest.SetUseOrigDistribution(kTRUE);
    }
      

    MFillH filltest(fMatrixTest);
    filltest.SetFilter(&seltest);
    filltest.SetName("fillMatrixTest");



    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTest);

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

    tlist.AddToList(&read);
    tlist.AddToList(&seltest);
    tlist.AddToList(&filltest);

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

    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTest");
    evtloop.SetProgressBar(&bar);

    if (!evtloop.Eventloop())
      return kFALSE;

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTest->Print("SizeCols");
    const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
    if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest))
    {
      *fLog << "MFindSupercutsONOFF::DefineTestMatrix; no.of generated events ("
	    << howmanygenerated 
            << ") is incompatible with the no.of requested events ("
            << howmanytest << ")" << endl;
    }

    *fLog << "test matrix was filled" << endl;
    *fLog << "=============================================" << endl;

    //--------------------------
    // write out test matrix

    if (filetest != "")
    {
      TFile filete(filetest, "RECREATE", "");
      fMatrixTest->Write();
      filete.Close();

      *fLog << "MFindSupercutsONOFF::DefineTestMatrix; Test matrix was written onto file '"
            << filetest << "'" << endl;
    }


    
    return kTRUE;
}



// --------------------------------------------------------------------------
//
// Define the matrices for the training and test sample respectively
//
//
//
Bool_t MFindSupercutsONOFF::DefineTrainTestMatrix(
    const TString &name, MH3 &href,
    const Int_t howmanytrain, const Int_t howmanytest,
    const TString &filetrain, const TString &filetest)
{
    *fLog << "=============================================" << endl;
    *fLog << "fill training and test matrix from file '" << name 
          << "',   select "   << howmanytrain 
          << " training and " << howmanytest << " test events " << endl;
    if (!fUseOrigDistribution)
    {
      *fLog << "     according to a distribution given by the MH3 object '"
            << href.GetName() << "'" << endl;
    }
    else
    {
      *fLog << "     randomly" << endl;
    }


    MParList  plist;
    MTaskList tlist;

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

     MFEventSelector2 selector(href);
    selector.SetNumMax(howmanytrain+howmanytest);
    selector.SetName("selectTrainTest");
    selector.SetInverted();
    if (fUseOrigDistribution)
    {
	*fLog << "MFindSupercutsONFF; MFEventSelector2::SetUseOrigDistribution(Bool)"
	      << " DOES NOT EXIST NOW..." << endl;
	// selector.SetUseOrigDistribution(kTRUE);
    }

    MContinue cont(&selector);
    cont.SetName("ContTrainTest");

    Double_t prob =  ( (Double_t) howmanytrain )
                   / ( (Double_t)(howmanytrain+howmanytest) );
    MFRandomSplit split(prob);

    MFillH filltrain(fMatrixTrain);
    filltrain.SetFilter(&split);
    filltrain.SetName("fillMatrixTrain");


    // consider this event as candidate for a test event 
    // only if event was not accepted as a training event

    MContinue conttrain(&split);
    conttrain.SetName("ContTrain");

    MFillH filltest(fMatrixTest);
    filltest.SetName("fillMatrixTest");


    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTrain);
    plist.AddToList(fMatrixTest);

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


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

    tlist.AddToList(&read);
    tlist.AddToList(&cont);

    tlist.AddToList(&split);
    tlist.AddToList(&filltrain);
    tlist.AddToList(&conttrain);

    tlist.AddToList(&filltest);

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

    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTrainTest");
    evtloop.SetProgressBar(&bar);

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

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTrain->Print("SizeCols");

    const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
    if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
    {
      *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events ("
	    << generatedtrain 
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }

    fMatrixTest->Print("SizeCols");
    const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
    if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
    {
      *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events ("
	    << generatedtest 
            << ") is incompatible with the no.of requested events ("
            << howmanytest << ")" << endl;
    }


    *fLog << "training and test matrix were filled" << endl;
    *fLog << "=============================================" << endl;


    //--------------------------
    // write out training matrix

    if (filetrain != "")
    {
      TFile filetr(filetrain, "RECREATE");
      fMatrixTrain->Write();
      filetr.Close();

      *fLog << "MFindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '"
            << filetrain << "'" << endl;
    }

    //--------------------------
    // write out test matrix

    if (filetest != "")
    {
      TFile filete(filetest, "RECREATE", "");
      fMatrixTest->Write();
      filete.Close();

      *fLog << "MFindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '"
            << filetest << "'" << endl;
    }

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Define the matrices for the training and test OFF sample respectively
//
//
//

Bool_t MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange(
    const TString &name, 
    const Double_t whichfractiontrain, 
    const Double_t whichfractiontest,  
    Double_t ThetaMin, Double_t ThetaMax, 
    const TString &filetrain, const TString &filetest)
    


{
    *fLog << "=============================================" << endl;
    *fLog << "Fill training and testing OFF matrices from file '" << name 
          << "',   select a fraction of " << whichfractiontrain 
          << " events for the training and a fraction of " << endl
	  << whichfractiontest << " for the testing" << endl;


    MParList  plist;
    MTaskList tlist;

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


    

    
    // Cuts in Theta

    // TString ThetaCutMinString ("ThetaOrig.fVal>");
    TString ThetaCutMinString ("MPointingPos.fZd>"); // new! // for magic
    ThetaCutMinString += ThetaMin;
    MContinue ThetaCutMin(ThetaCutMinString);
    ThetaCutMin.SetInverted();

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


    
    // Cuts in Size,

    TString SizeCutMinString ("MHillas.fSize>"); 
    SizeCutMinString += fSizeCutLow;
    MContinue SizeCutMin(SizeCutMinString);
    SizeCutMin.SetInverted();

    
    TString SizeCutMaxString ("MHillas.fSize<"); 
    SizeCutMaxString += fSizeCutUp;
    MContinue SizeCutMax(SizeCutMaxString);
    SizeCutMax.SetInverted();
    
    // Cuts in Dist
    TString DistCutMinString ("MHillasSrc.fDist>"); 
    DistCutMinString += fDistCutLow;
    MContinue DistCutMin(DistCutMinString);
    DistCutMin.SetInverted();

    TString DistCutMaxString ("MHillasSrc.fDist<"); 
    DistCutMaxString += fDistCutUp; 
    MContinue DistCutMax(DistCutMaxString);
    DistCutMax.SetInverted();


    // Cuts in leakage
    TString LeakCutMaxString ("MNewImagePar.fLeakage1<"); 
    LeakCutMaxString += fLeakageMax; 
    MContinue LeakCutMax(LeakCutMaxString);
    LeakCutMax.SetInverted();




    Double_t prob = whichfractiontrain /(whichfractiontrain+whichfractiontest);
                   
    
    

    MFRandomSplit split(prob);

    MFillH filltrain(fMatrixTrainOFF);
    filltrain.SetName("fillMatrixTrainOFF");
    filltrain.SetFilter(&split);
    
    
    // consider this event as candidate for a test event 
    // only if event was not accepted as a training event

    MContinue conttrain(&split);
    conttrain.SetName("ContTrain");

    MFillH filltest(fMatrixTestOFF);
    filltest.SetName("fillMatrixTestOFF");
    

    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTrainOFF);
    plist.AddToList(fMatrixTestOFF);

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

    tlist.AddToList(&read);
    tlist.AddToList(&ThetaCutMin);
    tlist.AddToList(&ThetaCutMax);
    
    
    tlist.AddToList(&SizeCutMin);
    tlist.AddToList(&SizeCutMax);

    tlist.AddToList(&DistCutMin);
    tlist.AddToList(&DistCutMax);
    tlist.AddToList(&LeakCutMax);

    
    tlist.AddToList(&split);
    tlist.AddToList(&filltrain);

    tlist.AddToList(&conttrain);
    tlist.AddToList(&filltest);

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

    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTrainTestOFF");
    evtloop.SetProgressBar(&bar);
    
    Int_t maxev = -1;
    if (!evtloop.Eventloop(maxev))
      return kFALSE;

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTrainOFF->Print("SizeCols");

    
    fMatrixTestOFF->Print("SizeCols");
    

    *fLog << "train and test matrix OFF were filled" << endl;
    *fLog << "=============================================" << endl;
    
    
    //--------------------------
    // write out training matrix
    
    if (filetrain != "")
    {
	TFile filetr(filetrain, "RECREATE");
	fMatrixTrainOFF->Write();
	filetr.Close();
	
	*fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange; Training matrix was written onto file '"
	      << filetrain << "'" << endl;
    }
    
    //--------------------------
    // write out test matrix

    if (filetest != "")
    {
	TFile filete(filetest, "RECREATE", "");
	fMatrixTestOFF->Write();
	filete.Close();
	
	*fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange; Test matrix was written onto file '"
	      << filetest << "'" << endl;
    }
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Define the matrices for the training and test sample respectively
//
//
//
Bool_t MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange(
    const TString &name, 
    const Double_t whichfractiontrain, 
    const Double_t whichfractiontest,  
    Double_t ThetaMin, Double_t ThetaMax, 
    const TString &filetrain, const TString &filetest)

   
{
    *fLog << "=============================================" << endl;
    *fLog << "Fill training and testing ON matrices from file '" << name 
          << "',   select a fraction of " << whichfractiontrain 
          << " events for the training and a fraction of " << endl
	  << whichfractiontest << " for the testing" << endl;


    MParList  plist;
    MTaskList tlist;

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


     // Cuts in Theta

    //TString ThetaCutMinString ("ThetaOrig.fVal>");
    TString ThetaCutMinString ("MPointingPos.fZd>"); // for magic
    ThetaCutMinString += ThetaMin;
    MContinue ThetaCutMin(ThetaCutMinString);
    ThetaCutMin.SetInverted();

    //TString ThetaCutMaxString ("ThetaOrig.fVal<");
    TString ThetaCutMaxString ("MPointingPos.fZd<"); // for magic
    ThetaCutMaxString += ThetaMax;
    MContinue ThetaCutMax(ThetaCutMaxString);
    ThetaCutMax.SetInverted();


    
    // Cuts in Size,

    TString SizeCutMinString ("MHillas.fSize>"); 
    SizeCutMinString += fSizeCutLow;
    MContinue SizeCutMin(SizeCutMinString);
    SizeCutMin.SetInverted();

    
    TString SizeCutMaxString ("MHillas.fSize<"); 
    SizeCutMaxString += fSizeCutUp;
    MContinue SizeCutMax(SizeCutMaxString);
    SizeCutMax.SetInverted();
    

    // Cuts in Dist
    TString DistCutMinString ("MHillasSrc.fDist>");
    DistCutMinString += fDistCutLow;
    MContinue DistCutMin(DistCutMinString);
    DistCutMin.SetInverted();

    TString DistCutMaxString ("MHillasSrc.fDist<");
    DistCutMaxString += fDistCutUp;
    MContinue DistCutMax(DistCutMaxString);
    DistCutMax.SetInverted();


    // Cuts in leakage
    TString LeakCutMaxString ("MNewImagePar.fLeakage1<");
    LeakCutMaxString += fLeakageMax;
    MContinue LeakCutMax(LeakCutMaxString);
    LeakCutMax.SetInverted();


    Double_t prob = whichfractiontrain/(whichfractiontrain + whichfractiontest);

    

    MFRandomSplit split(prob);

    MFillH filltrain(fMatrixTrain);
    filltrain.SetName("fillMatrixTrain");
    filltrain.SetFilter(&split);
    
    
    // consider this event as candidate for a test event 
    // only if event was not accepted as a training event

    MContinue conttrain(&split);
    conttrain.SetName("ContTrain");

    MFillH filltest(fMatrixTest);
    filltest.SetName("fillMatrixTest");


  

    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTrain);
    plist.AddToList(fMatrixTest);

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

    tlist.AddToList(&read);
    tlist.AddToList(&ThetaCutMin);
    tlist.AddToList(&ThetaCutMax);
    
       
    
    tlist.AddToList(&SizeCutMin);
    tlist.AddToList(&SizeCutMax);

    tlist.AddToList(&DistCutMin);
    tlist.AddToList(&DistCutMax);
    tlist.AddToList(&LeakCutMax);

    

    tlist.AddToList(&split);
    tlist.AddToList(&filltrain);

    tlist.AddToList(&conttrain);
    tlist.AddToList(&filltest);

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

    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTrainTest");
    evtloop.SetProgressBar(&bar);

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

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTrain->Print("SizeCols");
    
    fMatrixTest->Print("SizeCols");
    

    *fLog << "training and test matrix were filled" << endl;
    *fLog << "=============================================" << endl;


    //--------------------------
    // write out training matrix

    if (filetrain != "")
    {
      TFile filetr(filetrain, "RECREATE");
      fMatrixTrain->Write();
      filetr.Close();

      *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange; Train matrix was written onto file '"
            << filetrain << "'" << endl;
    }


     //--------------------------
    // write out test matrix

    if (filetest != "")
    {
      TFile filete(filetest, "RECREATE", "");
      fMatrixTest->Write();
      filete.Close();

      *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange; Test matrix was written onto file '"
            << filetest << "'" << endl;
    }

  return kTRUE;
}





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

// --------------------------------------------------------------------------
//
// Read only training matrices ON and OFF
//
//

Bool_t MFindSupercutsONOFF::ReadMatrixTrain(const TString &filetrainON, const TString &filetrainOFF)
{
  //--------------------------
  // read in training matrix ON sample

  TFile filetr(filetrainON);
  fMatrixTrain->Read("MatrixTrain");
  fMatrixTrain->Print("SizeCols");

  *fLog << "MFindSupercuts::ReadMatrixTrain; Training matrix for ON sample was read in from file '"
        << filetrainON << "'" << endl;
  filetr.Close();


  // read in training matrix OFF sample

  TFile filetrOFF(filetrainOFF);
  fMatrixTrainOFF->Read("MatrixTrainOFF");
  fMatrixTrainOFF->Print("SizeCols");

  *fLog << "MFindSupercutsONOFF::ReadMatrixTrain; Training matrix for OFF sample was read in from file '"
        << filetrainOFF << "'" << endl;
  filetrOFF.Close();

  return kTRUE;  

}

// --------------------------------------------------------------------------
//
// Read only test matrices ON and OFF
//
//

Bool_t MFindSupercutsONOFF::ReadMatrixTest(const TString &filetestON, const TString &filetestOFF)
{
  //--------------------------
  // read in testing matrix ON sample

//--------------------------
  // read in test matrix for ON sample
  
  TFile filete(filetestON);
  fMatrixTest->Read("MatrixTest");
  fMatrixTest->Print("SizeCols");

  *fLog << "MFindSupercuts::ReadMatrixTest; Test matrix for ON sample was read in from file '"
        << filetestON << "'" << endl;
  filete.Close();
  

   //--------------------------
  // read in test matrix for OFF sample
  
  TFile fileteOFF(filetestOFF);
  fMatrixTestOFF->Read("MatrixTestOFF");
  fMatrixTestOFF->Print("SizeCols");

  *fLog << "MFindSupercutsONOFF::ReadMatrixTest; Test matrix for OFF sample was read in from file '"
        << filetestOFF << "'" << endl;
  filete.Close();



  return kTRUE;  

}


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

// --------------------------------------------------------------------------
//
// Read training and test matrices from files
//
//

Bool_t MFindSupercutsONOFF::ReadMatrix(const TString &filetrain, const TString &filetest)
{
  //--------------------------
  // read in training matrix

  TFile filetr(filetrain);
  fMatrixTrain->Read("MatrixTrain");
  fMatrixTrain->Print("SizeCols");

  *fLog << "MFindSupercuts::ReadMatrix; Training matrix was read in from file '"
        << filetrain << "'" << endl;
  filetr.Close();


  //--------------------------
  // read in test matrix
  
  TFile filete(filetest);
  fMatrixTest->Read("MatrixTest");
  fMatrixTest->Print("SizeCols");

  *fLog << "MFindSupercuts::ReadMatrix; Test matrix was read in from file '"
        << filetest << "'" << endl;
  filete.Close();
  
  return kTRUE;  
}


// Read training and test matrices OFF from files
//
//

Bool_t MFindSupercutsONOFF::ReadMatrixOFF(const TString &filetrain, const TString &filetest)
{
  //--------------------------
  // read in training matrix

  TFile filetr(filetrain);
  fMatrixTrainOFF->Read("MatrixTrainOFF");
  fMatrixTrainOFF->Print("SizeCols");

  *fLog << "MFindSupercutsONOFF::ReadMatrixOFF; Training matrix OFF was read in from file '"
        << filetrain << "'" << endl;
  filetr.Close();


  //--------------------------
  // read in test matrix
  
  TFile filete(filetest);
  fMatrixTestOFF->Read("MatrixTestOFF");
  fMatrixTestOFF->Print("SizeCols");

  *fLog << "MFindSupercutsONOFF::ReadMatrixONOFF; Test matrix OFF was read in from file '"
        << filetest << "'" << endl;
  filete.Close();
  
  return kTRUE;  
}



// Function to compute the normalization factor for Train sample.
// The normalization factor is defined as the ratio of OFF/ON events.

Bool_t  MFindSupercutsONOFF::ComputeNormFactorTrain()
{
   Int_t EventsInTrainMatrixON = fMatrixTrain->GetM().GetNrows(); 
   Int_t EventsInTrainMatrixOFF = fMatrixTrainOFF->GetM().GetNrows();

   // Info...
   *fLog << "MFindSupercutsONOFF::ComputeNormFactorTrain; " << endl
	 << "EventsInTrainMatrixON = " << EventsInTrainMatrixON 
	 << " , EventsInTrainMatrixOFF = " << EventsInTrainMatrixOFF << endl;


   fNormFactorTrain = double(EventsInTrainMatrixON)/double(EventsInTrainMatrixOFF);

   
   *fLog << "MFindSupercutsONOFF::ComputeNormFactorTrain; fNormFactorTrain is : " 
	 << fNormFactorTrain << endl;
   
   return kTRUE;
}




// Function to compute the normalization factor for Test sample.
// The normalization factor is defined as the ratio of OFF/ON events.

Bool_t  MFindSupercutsONOFF::ComputeNormFactorTest()
{
   Int_t EventsInTestMatrixON = fMatrixTest->GetM().GetNrows(); 
   Int_t EventsInTestMatrixOFF = fMatrixTestOFF->GetM().GetNrows();

   // Info...
   *fLog << "MFindSupercutsONOFF::ComputeNormFactorTest; " << endl
	 << "EventsInTestMatrixON = " << EventsInTestMatrixON 
	 << " , EventsInTestMatrixOFF = " << EventsInTestMatrixOFF << endl;

   fNormFactorTest = double(EventsInTestMatrixON)/double(EventsInTestMatrixOFF);
   
   *fLog << "MFindSupercutsONOFF::ComputeNormFactorTest; fNormFactorTest is : " 
	 << fNormFactorTest << endl;

   return kTRUE;
}


Bool_t  MFindSupercutsONOFF::SetGammaEfficiency(Double_t gammaeff)
{

    if (gammaeff < 0.0 || gammaeff >1.0)
    {
	*fLog << "MFindSupercutsONOFF::SetGammaEfficiency; " << endl
	      << "The argument of SetGammaEfficiency ("
	      << gammaeff << ") is outside the range 0-1" << endl
	      << "This is not allowed..." << endl;
	return kFALSE;
	    
    }

    fGammaEfficiency = gammaeff;
    
    // TEMP INFO
    cout << "fGammaEfficiency has been set to : " << fGammaEfficiency << endl;
    // END TEMP
    return kTRUE;
	

}

//------------------------------------------------------------------------
//
// Steering program for optimizing the supercuts
// ---------------------------------------------
//
//      the criterion for the 'optimum' is 
//
//          - a maximum significance of the gamma signal in the alpha plot, 
//            in which the supercuts have been applied
//
// The various steps are :
//
// - setup the 2 event loops (ON and OFF) to be executed for each call to fcnSupercuts 
// - call TMinuit to do the minimization of (-significance) :
//        the fcnSupercuts function calculates the significance 
//                                             for the current cut values
//        for this - the 2 alpha plots (ON and OFF) are produced in the 2 event loops, 
//                   in which the SAME cuts have been applied
//                 - the significance of the gamma signal is extracted from the 
//                   the 2 alpha plots (ON OFF) using MHFindSignificanceONOFF::FindSigmaONOFF.
//                   
//
// Needed as input : (to be set by the Set functions)
//
// - fFilenameParam      name of file to which optimum values of the 
//                       parameters are written
//
// - fHadronnessName     name of container where MSupercutsCalcONOFF will
//                       put the supercuts hadronness for ON data

// - fHadronnessNameOFF  name of container where MSupercutsCalcONOFF will
//                       put the supercuts hadronness for OFF data
//
// - fAlphaSig, fAlphaBkgMin, fAlphaBkgMax ; Signal and Bkg region in alpha histo. 

// - for the minimization, the starting values of the parameters are taken  
//     - from file parSCinit (if it is != "")
//     - or from the arrays params and/or steps 
//
//----------------------------------------------------------------------
Bool_t MFindSupercutsONOFF::FindParams(TString parSCinit,
					  TArrayD &params, TArrayD &steps)
{
    // Setup the event loops which will be executed in the 
    //                 fcnSupercuts function  of MINUIT
    //
    // parSCinit is the name of the file containing the initial values 
    // of the parameters; 
    // if parSCinit = ""   'params' and 'steps' are taken as initial values

    *fLog << "=============================================" << endl;
    *fLog << "Setup event loop for fcnSupercuts" << endl;

    if (fHadronnessName.IsNull())
    {
      *fLog << "MFindSupercutsONOFF::FindParams; hadronness name for ON data is not defined... aborting"
            << endl;
      return kFALSE;
    }


    if (fHadronnessNameOFF.IsNull())
    {
      *fLog << "MFindSupercutsONOFF::FindParams; hadronness name for OFF data is not defined... aborting"
            << endl;
      return kFALSE;
    }


    if (fMatrixTrain == NULL)
    {
      *fLog << "MFindSupercuts::FindParams; training matrix is not defined... aborting"
            << endl;
      return kFALSE;
    }

    if (fMatrixTrain->GetM().GetNrows() <= 0)
    {
      *fLog << "MFindSupercuts::FindParams; training matrix has no entries"
            << endl;
      return kFALSE;
    }

    

    if (fMatrixTrainOFF == NULL)
    {
      *fLog << "MFindSupercutsONOFF::FindParams; training matrix OFF is not defined... aborting"
            << endl;
      return kFALSE;
    }

    if (fMatrixTrainOFF->GetM().GetNrows() <= 0)
    {
      *fLog << "MFindSupercutsONOFF::FindParams; training matrix OFF has no entries"
            << endl;
      return kFALSE;
    }

    if (fAlphaDistributionsRootFilename.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::FindParams; fAlphaDistributionsRootFilename is not defined... program aborting..." << endl;

	return kFALSE;
    }

   

    //--------------------------------
    // create container for the supercut parameters
    // and set them to their initial values
    MSupercuts super;

    // take initial values from file parSCinit
    if (parSCinit != "")
    {
      TFile inparam(parSCinit);
      super.Read("MSupercuts");
      inparam.Close();
      *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from file "
            << parSCinit << endl;
    }

    // take initial values from 'params' and/or 'steps'
    else if (params.GetSize() != 0  || steps.GetSize()  != 0 )
    {
	if (params.GetSize()  != 0)
	{
	    *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from 'params'"
		  << endl;
	    super.SetParameters(params);
	}
	if (steps.GetSize()  != 0)
	{
	    *fLog << "MFindSupercutsONOFF::FindParams; initial step sizes are taken from 'steps'"
		  << endl;
	    super.SetStepsizes(steps);
	}
	
	/*
	// TMP
	// Print parameters
	if (params.GetSize() == steps.GetSize())
	{
	    *fLog << "MFindSupercutsONOFF; SC parameters and Setps are: " << endl;
	    for (Int_t z = 0; z < params.GetSize(); z++)
	    {
		cout << params[z] << setw(20) << steps[z] << endl;
	    }
	}

	// ENDTMP

	*/

	/*
	// TMP2

	TArrayD paramsBis = super.GetParameters();
	TArrayD stepsBis = super.GetStepsizes();
	if (paramsBis.GetSize() == stepsBis.GetSize())
	{
	    *fLog << "MFindSupercutsONOFF; SC parametersBis and SetpsBis are: " << endl;
	    for (Int_t z = 0; z < paramsBis.GetSize(); z++)
	    {
		cout << paramsBis[z] << setw(20) << stepsBis[z] << endl;
	    }
	}
	
	// ENDTMP2
	*/

	else
	{
	    *fLog << "MFindSupercutsONOFF; ERROR: params.GetSize() != steps.GetSize() " 
		  << endl;
	    
	}



    }
    else
    {
        *fLog << "MFindSupercutsONOFF::FindParams; initial values and step sizes are taken from the MSupercuts constructor"
              << endl;
    }


    // Computation of fNormFactorTrain and creation of a container for 
    // storing this value and include it  in the MParList for being 
    // used by function fcnSupercuts

    if (!ComputeNormFactorTrain())
    {
	*fLog << "Normalization factor for train sample (ON-OFF) could not be computed. Aborting..." << endl;
	return kFALSE;
    }

    
    *fLog << "MFindSupercutsONOFF::FindParams; fNormFactorTrain = " << fNormFactorTrain << endl;

    

    MDataValue NormFactorCont(fNormFactorTrain);
    NormFactorCont.SetName("NormFactorTrain");
                                 

    // Value of fAlphaSig, fAlphaBkgMin and fAlphaBkgMax 
    // are stored in MDataValue containers,
    // then they will be passed to the MParList and will be 
    // accessible to fcnsupercuts

    MDataValue AlphaSigCont(fAlphaSig);
    AlphaSigCont.SetName("AlphaSigValue");

    MDataValue AlphaBkgMinCont(fAlphaBkgMin);
    AlphaBkgMinCont.SetName("AlphaBkgMinValue");

    MDataValue AlphaBkgMaxCont(fAlphaBkgMax);
    AlphaBkgMaxCont.SetName("AlphaBkgMaxValue");
    

    MDataValue DegreeONCont (fDegreeON);
    DegreeONCont.SetName("DegreeON");

    MDataValue DegreeOFFCont (fDegreeOFF);
    DegreeOFFCont.SetName("DegreeOFF");




    
    // Value of fUseFittedQuantities is stored in container 
    // and passed to the MParList and will be 
    // accessible to fcnsupercuts

    

    MDataValue UseFittedQuantitiesCont(fUseFittedQuantities);
    UseFittedQuantitiesCont.SetName("UseFittedQuantitiesValue");


     // Value of fNormFactorFromAlphaBkg is stored in container 
    // and passed to the MParList and will be 
    // accessible to fcnsupercuts

    MDataValue UseNormFactorFromAlphaBkgCont(fNormFactorFromAlphaBkg);
    UseNormFactorFromAlphaBkgCont.SetName("UseNormFactorFromAlphaBkgValue");



    // A vector of pointers to objects of the class MEvtLoop is defined.
    // One of the pointed loops will be used to compute ON data alpha plot
    // and the other for computing OFF data alpha plot.

    MEvtLoop evtloopfcn[2] = {MEvtLoop("ONDataEvtLoopFCN"), 
			      MEvtLoop("OFFDataEvtLoopFCN")};




    // ******************************************************************
    // EVENT LOOP FOR COMPUTING ALPHA HISTOGRAM FOR ON DATA IS SET
    // ******************************************************************

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

    MParList  parlistfcn;
    MTaskList tasklistfcn;

    // loop over rows of matrix
    MMatrixLoop loop(fMatrixTrain);


    //--------------------------------
    // calculate supercuts hadronness
    fCalcHadTrain->SetHadronnessName(fHadronnessName);

    // Boolean variable that controls in class MSupercutsCalcONOFF
    // wether the supercuts are stored or not
    // is set to the default value, kFALSE. (no storage of supercuts)

    fCalcHadTrain -> SetStoreAppliedSupercuts(kFALSE);


    
    // Set boolean variable that controls wether cuts are 
    // dynamic or static

    fCalcHadTrain -> SetVariableUseStaticCuts(fUseStaticCuts);

    if (!fUseStaticCuts)
    {
	// Set boolean variable that controls wether the theta variable 
	// is used or not in the computation of the dynamical cuts
	
	fCalcHadTrain -> SetVariableNotUseTheta(fNotUseTheta);
    }

    // apply the supercuts
    MF scfilter(fHadronnessName+".fHadronness>0.5");
    MContinue supercuts(&scfilter);

    // plot |alpha|
    const TString  mh3Name = "AlphaFcn";
    MBinning binsalpha("Binning"+mh3Name);
    //binsalpha.SetEdges(54, -12.0, 96.0);
    binsalpha.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
     
     
    *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl;

    // |alpha| is assumed to be in column 7 of the matrix
    MH3 alpha("MatrixTrain[7]");
    alpha.SetName(mh3Name);

    MFillH fillalpha(&alpha);

    // book histograms to be filled during the optimization
    //                              ! not in the event loop !
    MHCT1Supercuts plotsuper;
    plotsuper.SetupFill(&parlistfcn);

    //******************************
    // entries in MParList (extension of old MParList)
    
    parlistfcn.AddToList(&tasklistfcn);
    parlistfcn.AddToList(&super);
    parlistfcn.AddToList(&NormFactorCont);
    parlistfcn.AddToList(&AlphaSigCont);
    parlistfcn.AddToList(&AlphaBkgMinCont);
    parlistfcn.AddToList(&AlphaBkgMaxCont);
    parlistfcn.AddToList(&DegreeONCont);
    parlistfcn.AddToList(&UseFittedQuantitiesCont);
    parlistfcn.AddToList(&UseNormFactorFromAlphaBkgCont);

    parlistfcn.AddToList(fCam);
    parlistfcn.AddToList(fMatrixTrain);

    parlistfcn.AddToList(&binsalpha);
    parlistfcn.AddToList(&alpha);

    parlistfcn.AddToList(&plotsuper);

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

    tasklistfcn.AddToList(&loop);
    tasklistfcn.AddToList(fCalcHadTrain);
    tasklistfcn.AddToList(&supercuts);
    tasklistfcn.AddToList(&fillalpha);

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


    // &evtloopfcn[0] = new MEvtLoop("ONDataEvtLoopFCN");
    evtloopfcn[0].SetParList(&parlistfcn);
    // MEvtLoop evtloopfcn("EvtLoopFCN");
    // evtloopfcn.SetParList(&parlistfcn);
    *fLog << "Event loop for computing ON data alpha histo in fcnSupercuts has been setup" << endl;


//---- End of setting of loop for compuing ON data alpha plot ---------------

    

     // ******************************************************************
    // EVENT LOOP FOR COMPUTING ALPHA HISTOGRAM FOR OFF DATA IS SET
    // ******************************************************************

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

    MParList  parlistfcnOFF;
    MTaskList tasklistfcnOFF;

    // loop over rows of matrix
    MMatrixLoop loopOFF(fMatrixTrainOFF);


    //--------------------------------
    // calculate supercuts hadronness
    fCalcHadTrainOFF->SetHadronnessName(fHadronnessNameOFF);


    // Boolean variable that controls in class MSupercutsCalcONOFF
    // wether the supercuts are stored or not
    // is set to the default value, kFALSE. (no storage of supercuts)

    fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kFALSE);


    
    // Set boolean variable that controls wether cuts are 
    // dynamic or static

    fCalcHadTrainOFF -> SetVariableUseStaticCuts(fUseStaticCuts);


    if (!fUseStaticCuts)
    {
	// Set boolean variable that controls wether the theta variable 
	// is used or not in the computation of the dynamica cuts
	
	fCalcHadTrainOFF -> SetVariableNotUseTheta(fNotUseTheta);
    }

    // apply the supercuts
    MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5");
    MContinue supercutsOFF(&scfilterOFF);

    // plot |alpha|
    const TString  mh3NameOFF = "AlphaOFFFcn";
    MBinning binsalphaOFF("Binning"+mh3NameOFF);
    //binsalphaOFF.SetEdges(54, -12.0, 96.0);
    binsalphaOFF.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);


    *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl;

    // |alpha| is assumed to be in column 7 of the matrix
    MH3 alphaOFF("MatrixTrainOFF[7]");
    alphaOFF.SetName(mh3NameOFF);

    MFillH fillalphaOFF(&alphaOFF);

    

    //******************************
    // entries in MParList (extension of old MParList)
    
    parlistfcnOFF.AddToList(&tasklistfcnOFF);
    parlistfcnOFF.AddToList(&super);
    parlistfcnOFF.AddToList(fCam);
    parlistfcnOFF.AddToList(fMatrixTrainOFF);

    parlistfcnOFF.AddToList(&binsalphaOFF);
    parlistfcnOFF.AddToList(&alphaOFF);

    parlistfcnOFF.AddToList(&DegreeOFFCont);
    

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

    tasklistfcnOFF.AddToList(&loopOFF);
    tasklistfcnOFF.AddToList(fCalcHadTrainOFF);
    tasklistfcnOFF.AddToList(&supercutsOFF);
    tasklistfcnOFF.AddToList(&fillalphaOFF);

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


    // &evtloopfcn[1] = new MEvtLoop("OFFDataEvtLoopFCN");
    evtloopfcn[1].SetParList(&parlistfcnOFF);
    // MEvtLoop evtloopfcn("EvtLoopFCN");
    // evtloopfcn.SetParList(&parlistfcn);
    *fLog << "Event loop for computing OFF data alpha histo in fcnSupercuts has been setup" << endl;


//---- End of setting of loop for compuing OFF data alpha plot ---------------


    // address of evtloopfcn is used in CallMinuit()


    //-----------------------------------------------------------------------
    //
    //----------   Start of minimization part   --------------------
    //
    // Do the minimization with MINUIT
    //
    // Be careful: This is not thread safe
    //
    *fLog << "========================================================" << endl;
    *fLog << "Start minimization for supercuts" << endl;


    // -------------------------------------------
    // prepare call to MINUIT
    //

    // get initial values of parameters 
    fVinit = super.GetParameters();
    fStep  = super.GetStepsizes();

    TString name[fVinit.GetSize()];
    fStep.Set(fVinit.GetSize());
    fLimlo.Set(fVinit.GetSize());
    fLimup.Set(fVinit.GetSize());
    fFix.Set(fVinit.GetSize());

    fNpar = fVinit.GetSize();

    for (UInt_t i=0; i<fNpar; i++)
    {
        name[i]   = "p";
        name[i]  += i+1;
        //fStep[i]  = TMath::Abs(fVinit[i]/10.0);
        fLimlo[i] = -100.0;
        fLimup[i] =  100.0;
        fFix[i]   =     0;
    }


    


    if(fSetLimitsToSomeMinuitParams)
      {

	// Limits for some Minuit parameters are set. For the time being the values are set in the constructor 
	// of the class. One MUST be very careful to set limits such that the expected final values 
	// (optimized values) are far away from the limits.
	// Limits are set to those Minuit parameters not depending in Size,dist and theta (static values). 
	// Limits are set only to Dist, Length and Width.

	
	*fLog << endl 
	      <<"Limits set to hillas parameters Length, Width, Dist and Leakage1 for the optimization with Minuit." 
	      << endl;
	

	
	fLimup[0] = fMinuitLengthUPUpperLimit;
	fLimlo[0] = fMinuitLengthUPLowerLimit;
	fLimup[8] = fMinuitLengthLOWUpperLimit;
	fLimlo[8] = fMinuitLengthLOWLowerLimit;
	fLimup[16] = fMinuitWidthUPUpperLimit;
	fLimlo[16] = fMinuitWidthUPLowerLimit;
	fLimup[24] = fMinuitWidthLOWUpperLimit;
	fLimlo[24] = fMinuitWidthLOWLowerLimit;
	fLimup[32] = fMinuitDistUPUpperLimit;
	fLimlo[32] = fMinuitDistUPLowerLimit;
	fLimup[40] = fMinuitDistLOWUpperLimit;
	fLimlo[40] = fMinuitDistLOWLowerLimit;

	fLimup[80] = fMinuitLeakage1UPUpperLimit;
	fLimlo[80] = fMinuitLeakage1UPLowerLimit;

	

      }
    
   

    // Wolfgang sets manyally some parameters.
    // That might give an error if parameters are changed...
    // To be changed in future...

    // these parameters make no sense, fix them at 0.0
    fVinit[33] = 0.0;
    fStep[33]  = 0.0;
    fFix[33]   = 1;

    fVinit[36] = 0.0;
    fStep[36]  = 0.0;
    fFix[36]   = 1;

    fVinit[39] = 0.0;
    fStep[39]  = 0.0;
    fFix[39]   = 1;

    fVinit[41] = 0.0;
    fStep[41]  = 0.0;
    fFix[41]   = 1;

    fVinit[44] = 0.0;
    fStep[44]  = 0.0;
    fFix[44]   = 1;

    fVinit[47] = 0.0;
    fStep[47]  = 0.0;
    fFix[47]   = 1;


   
    // ADDITIONAL PARAMETERS ARE FIXED ACCORDING TO THE 
    // VALUES OF THE BOOLEAN VARIABLES fUseStaticCuts, AND
    // fNotUseTheta


    if (fUseStaticCuts)
    {    
        // Static cuts will be used; all parameters with index 
	// 1-7 are fixed for ALL variables (width, length, ...)
	
       for (UInt_t i=0; i<fNpar; i = i+8)
       {
	   // At some point, the number of paraemters in the 
	   // dynamical cuts parameterezattion (currently 7) 
           // will become a variable that can be set from 
	   // outside the class

	   for (UInt_t j=i+1; j<=i+7; j++)
	   {
	       fVinit[j] = 0.0;
	       fStep[j]  = 0.0;
	       fFix[j]   = 1;
	   }
       }
    }
    else
    {
	if(fNotUseTheta)
	{
	    // Theta is not used in the parameterization of the dynamical cut
	    // Parameters with index 2 and 5 are fixed for all variables 
	    // (width, length, ...)

	    for (UInt_t i=0; i<fNpar; i = i+8)
	    {
		// At some point, the number of paraemters in the 
		// dynamical cuts parameterezattion (currently 7) 
		// will become a variable that can be set from 
		// outside the class
		
		for (UInt_t j=i+2; j<=i+7; j = j+3)
		{
		    fVinit[j] = 0.0;
		    fStep[j]  = 0.0;
		    fFix[j]   = 1;
		}
	    }


	}

	if(!fUseDist)
	  {
	    // parameterization using dist is removed
	    // Parameters with index 1 and 4 and 7 are fixed for all variables 
	    // (width, length, ...)
	     for (UInt_t i=0; i<fNpar; i = i+8)
	    {
		// At some point, the number of paraemters in the 
		// dynamical cuts parameterezattion (currently 7) 
		// will become a variable that can be set from 
		// outside the class
		
		for (UInt_t j=i+1; j<=i+7; j = j+3)
		{
		    fVinit[j] = 0.0;
		    fStep[j]  = 0.0;
		    fFix[j]   = 1;
		}
	    }



	  }




    }






    // vary only first 48 parameters
    //for (UInt_t i=0; i<fNpar; i++)
    //{
    //    if (i >= 48)
    //	{
    //      fStep[i] = 0.0;
    //      fFix[i]  =   1;
    //	}
    //}
 
    // -------------------------------------------
    // call MINUIT
    
    Bool_t rc;

    if(fSkipOptimization)
    {
	*fLog << "Parameter optimization is skipped. Using previously optimized parameters."
	      << endl;
	

	rc = kTRUE;
	

	if(fUseInitialSCParams)
	{

	    // write Initial parameter values onto root file filenameParam
	    // so that later they can be applied on the data

	     *fLog << "Initial SC parameter values are written onto file '"
		  << fFilenameParam << "' :" << endl;

	    
	    TFile outparam(fFilenameParam, "RECREATE"); 
	    super.Write();
	    outparam.Close();
	    

	    
	     const TArrayD &check = super.GetParameters();
	    for (Int_t i=0; i<check.GetSize(); i++)
		*fLog << check[i] << ",  ";
	    *fLog << endl;

	
	}


    }
    else
    {
	TStopwatch clock;
	clock.Start();
	
	*fLog << "before calling CallMinuit" << endl;
    
	MMinuitInterface inter;               
	rc = inter.CallMinuit(fcnSupercuts, name,
			      fVinit, fStep, fLimlo, fLimup, fFix,
			      evtloopfcn, "SIMPLEX", kFALSE);
	
	*fLog << "after calling CallMinuit" << endl;
	
	*fLog << "Time spent for the minimization in MINUIT :   " << endl;;
	clock.Stop();
	clock.Print();

	// plotsuper.DrawClone();

	if (!rc)
	{
	    *fLog << "*************  ERROR !!! *********************" << endl
		  << "Minimization could not finish properly..." << endl
		  << "MMinuitInterface.CallMinuit() returned kFALSE" << endl;
	    return kFALSE;
	}
	
	
	
	
	*fLog << "Minimization for supercuts finished" << endl;
	*fLog << "========================================================" << endl;
	
	
	// -----------------------------------------------------------------
	// in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put 
	//                    into MSupercuts
	
	// write optimum parameter values onto file filenameParam
	
	TFile outparam(fFilenameParam, "RECREATE"); 
	super.Write();
	outparam.Close();
	
	*fLog << "Optimum parameter values for supercuts were written onto file '"
	      << fFilenameParam << "' :" << endl;
    
	const TArrayD &check = super.GetParameters();
	for (Int_t i=0; i<check.GetSize(); i++)
	    *fLog << check[i] << ",  ";
	*fLog << endl;
	
	
	
	
    }
    
	
    
    
    
    *fLog << "End of  supercuts optimization part" << endl;
    *fLog << "======================================================" << endl;
    
    
    *fLog << "Applying the optimized supercuts on the TRAIN sample and " 
	  << "producing plots related to the output results" << endl;

    *fLog << "======================================================" << endl;
    
    if (!TestParamsOnTrainSample())
    {
	*fLog << "MFindSupercutsONOFF::FindParams;" 
	      << "TestParamsOnTrainSample failed"
	      << endl;
	
	return kFALSE;
    }
    
    
    
    return kTRUE;
}







// -----------------------------------------------------------------------
//
// Test the supercuts on the test sample using ON and OFF data
// 2 loops are used to fill the 2 alpha distributions
// Supercuts applied (LengthUp, LengthLow...) 
// to all the individual events, as well as 
// the features of the shower images (Length,Width...) are 
// stored in 2 TTree objects (for ON and OFF data) in the 
// root file specified by variable "fAlphaDistributionsRootFilename"



Bool_t MFindSupercutsONOFF::TestParamsOnTestSample()
{
    if (fMatrixTest->GetM().GetNrows() <= 0)
    {
        *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; test matrix ON has no entries" 
              << endl;
        return kFALSE;
    }


    if (fMatrixTestOFF->GetM().GetNrows() <= 0)
    {
        *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; test matrix OFF has no entries" 
              << endl;
        return kFALSE;
    }

    // -------------   TEST the supercuts    ------------------------------
    //
    *fLog << "Test the supercuts on the test sample ON and OFF" << endl;

    // -----------------------------------------------------------------
    // read optimum parameter values from file filenameParam
    // into array 'supercutsPar'

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

    *fLog << "Optimum parameter values for supercuts were read from file '";
    *fLog << fFilenameParam << "' :" << endl;

    const TArrayD &supercutsPar = scin.GetParameters();
    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
        *fLog << supercutsPar[i] << ",  ";
    *fLog << endl;
    //---------------------------


    // -----------------------------------------------------------------
    if (fHadronnessName.IsNull())
    {
        *fLog << "MFindSupercutsONOFF::TestParamsONOFF; hadronness name for ON data is not defined... aborting";
        *fLog << endl;
        return kFALSE;
    }

    if (fHadronnessNameOFF.IsNull())
    {
        *fLog << "MFindSupercutsONOFF::TestParamsONOFF; hadronness name for OFF data is not defined... aborting";
        *fLog << endl;
        return kFALSE;
    }

    if(fTestONSupercutsAppliedName.IsNull())
    {
        *fLog << "MFindSupercutsONOFF::TestParamsONOFF; " 
	      << " MTSupercutsApplied tree name for ON Test data is not defined... aborting"
	      << endl;
        return kFALSE;
    }

    if(fTestOFFSupercutsAppliedName.IsNull())
    {
        *fLog << "MFindSupercutsONOFF::TestParamsONOFF; " 
	      << " MTSupercutsApplied container name for OFF Test data is not defined... aborting"
	      << endl;
        return kFALSE;
    }
    

    MSupercuts supercuts;
    supercuts.SetParameters(supercutsPar);

    // Objects containing the trees used to store supercuts applied
    // TTree Branches are also created
    MTSupercutsApplied supapptestON(fTestONSupercutsAppliedName.Data(), "");
    supapptestON.CreateTreeBranches();

    MTSupercutsApplied supapptestOFF(fTestOFFSupercutsAppliedName.Data(), "");
    supapptestOFF.CreateTreeBranches();

    

    
    fCalcHadTestOFF -> SetHadronnessName(fHadronnessNameOFF);
    fCalcHadTest -> SetHadronnessName(fHadronnessName);


    // Settings related to the storage of teh supercuts applied 
    fCalcHadTest -> SetSupercutsAppliedName(fTestONSupercutsAppliedName);
    fCalcHadTest -> SetStoreAppliedSupercuts(kTRUE);

    fCalcHadTestOFF -> SetSupercutsAppliedName(fTestOFFSupercutsAppliedName);
    fCalcHadTestOFF -> SetStoreAppliedSupercuts(kTRUE);
    


    // Set boolean variable that controls wether cuts are 
    // dynamic or static

    fCalcHadTest -> SetVariableUseStaticCuts(fUseStaticCuts);
    fCalcHadTestOFF -> SetVariableUseStaticCuts(fUseStaticCuts);

    if (!fUseStaticCuts)
    {
	// Set boolean variable that controls wether the theta variable 
	// is used or not in the computation of the dynamical cuts
	
	fCalcHadTest -> SetVariableNotUseTheta(fNotUseTheta);
	fCalcHadTestOFF -> SetVariableNotUseTheta(fNotUseTheta);
    }


     // apply the supercuts to OFF data

    MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5");
    MContinue applysupercutsOFF(&scfilterOFF);


     // apply the supercuts to ON data
    MF scfilter(fHadronnessName+".fHadronness>0.5");
    MContinue applysupercuts(&scfilter);

   

    // ****************************************************
    // Filling OFF alpha distribution 

    MParList parlistOFF;
    MTaskList tasklistOFF;

    MMatrixLoop loopOFF(fMatrixTestOFF);

    // plot alpha before applying the supercuts
    const TString  mh3NameOFFB = "AlphaOFFBefore";
    MBinning binsalphaOFFbef("Binning"+mh3NameOFFB);
    // binsalphaOFFbef.SetEdges(54, -12.0, 96.0);
    binsalphaOFFbef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);

    // |alpha| is assumed to be in column 7 of the matrix
    MH3 alphaOFFbefore("MatrixTestOFF[7]");
    alphaOFFbefore.SetName(mh3NameOFFB);


    MFillH fillalphaOFFbefore(&alphaOFFbefore);
    fillalphaOFFbefore.SetName("FillAlphaOFFBefore");

   

    // plot alpha OFF after applying the supercuts
    const TString  mh3NameOFFA = "AlphaOFFAfter";
    MBinning binsalphaOFFaft("Binning"+mh3NameOFFA);
    // binsalphaOFFaft.SetEdges(54, -12.0, 96.0);
    binsalphaOFFaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
    
    MH3 alphaOFFafter("MatrixTestOFF[7]");
    alphaOFFafter.SetName(mh3NameOFFA);


    MFillH fillalphaOFFafter(&alphaOFFafter);
    fillalphaOFFafter.SetName("FillAlphaOFFAfter");

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

    parlistOFF.AddToList(&tasklistOFF);

    parlistOFF.AddToList(&supercuts);

    parlistOFF.AddToList(&supapptestOFF);

    parlistOFF.AddToList(fCam);
    parlistOFF.AddToList(fMatrixTestOFF);

    parlistOFF.AddToList(&binsalphaOFFbef);
    parlistOFF.AddToList(&alphaOFFbefore);

    parlistOFF.AddToList(&binsalphaOFFaft);
    parlistOFF.AddToList(&alphaOFFafter);


    //******************************
    // Entries in MtaskList
    tasklistOFF.AddToList(&loopOFF);
    tasklistOFF.AddToList(&fillalphaOFFbefore);

    tasklistOFF.AddToList(fCalcHadTestOFF);
    tasklistOFF.AddToList(&applysupercutsOFF);

    tasklistOFF.AddToList(&fillalphaOFFafter);

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

    MProgressBar barOFF;
    MEvtLoop evtloopOFF;
    evtloopOFF.SetParList(&parlistOFF);
    evtloopOFF.SetName("EvtLoopTestParamsOFF");
    evtloopOFF.SetProgressBar(&barOFF);
    Int_t maxeventsOFF = -1;
    if (!evtloopOFF.Eventloop(maxeventsOFF))
        return kFALSE;

    tasklistOFF.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // draw the alpha plots

    MH3* alphaOFFBefore = (MH3*)parlistOFF.FindObject(mh3NameOFFB, "MH3");
    TH1  &alphaOFFHistb = alphaOFFBefore->GetHist();
    TString alphaOFFHistbName ("TestAlphaOFFBeforeCuts");
    
    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
	      << " fThetaRangeString is NOT defined..." << endl;
	alphaOFFHistbName += ("ThetaRangeStringUndefined");
    }
    else
    {
	alphaOFFHistbName += (fThetaRangeString);
    }


    alphaOFFHistb.SetXTitle("|\\alpha|  [\\circ]");
    alphaOFFHistb.SetName(alphaOFFHistbName);


    MH3* alphaOFFAfter = (MH3*)parlistOFF.FindObject(mh3NameOFFA, "MH3");
    TH1  &alphaOFFHista = alphaOFFAfter->GetHist();
    TString alphaOFFHistaName ("TestAlphaOFFAfterCuts");

    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
	      << " fThetaRangeString is NOT defined..." << endl;
	alphaOFFHistaName += ("ThetaRangeStringUndefined");
    }
    else
    {
	alphaOFFHistaName += (fThetaRangeString);
    }


    alphaOFFHista.SetXTitle("|\\alpha|  [\\circ]");
    alphaOFFHista.SetName(alphaOFFHistaName);

    



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



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


    MParList  parlist2;
    MTaskList tasklist2;

    MMatrixLoop loop(fMatrixTest);

    // plot alpha before applying the supercuts
    const TString  mh3NameB = "AlphaBefore";
    MBinning binsalphabef("Binning"+mh3NameB);
    //binsalphabef.SetEdges(54, -12.0, 96.0);
    binsalphabef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
    
    // |alpha| is assumed to be in column 7 of the matrix
    MH3 alphabefore("MatrixTest[7]");
    alphabefore.SetName(mh3NameB);


    MFillH fillalphabefore(&alphabefore);
    fillalphabefore.SetName("FillAlphaBefore");


    // plot alpha after applying the supercuts
    const TString  mh3NameA = "AlphaAfter";
    MBinning binsalphaaft("Binning"+mh3NameA);
    // binsalphaaft.SetEdges(54, -12.0, 96.0);
    binsalphaaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);

    MH3 alphaafter("MatrixTest[7]");
    alphaafter.SetName(mh3NameA);


    MFillH fillalphaafter(&alphaafter);
    fillalphaafter.SetName("FillAlphaAfter");

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

    parlist2.AddToList(&tasklist2);

    parlist2.AddToList(&supercuts);

    parlist2.AddToList(&supapptestON);

    parlist2.AddToList(fCam);
    parlist2.AddToList(fMatrixTest);

    parlist2.AddToList(&binsalphabef);
    parlist2.AddToList(&alphabefore);

    parlist2.AddToList(&binsalphaaft);
    parlist2.AddToList(&alphaafter);

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

    tasklist2.AddToList(&loop);
    tasklist2.AddToList(&fillalphabefore);

    tasklist2.AddToList(fCalcHadTest);
    tasklist2.AddToList(&applysupercuts);

    tasklist2.AddToList(&fillalphaafter);

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

    MProgressBar bar2;
    MEvtLoop evtloop2;
    evtloop2.SetParList(&parlist2);
    evtloop2.SetName("EvtLoopTestParams");
    evtloop2.SetProgressBar(&bar2);
    Int_t maxevents2 = -1;
    if (!evtloop2.Eventloop(maxevents2))
        return kFALSE;

    tasklist2.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // draw the alpha plots

    MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3");
    TH1  &alphaHist1 = alphaBefore->GetHist();
    TString alphaHist1Name ("TestAlphaBeforeCuts");
    
    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
	      << " fThetaRangeString is NOT defined..." << endl;
	alphaHist1Name += ("ThetaRangeStringUndefined");
    }
    else
    {
	alphaHist1Name += (fThetaRangeString);
    }




    alphaHist1.SetName(alphaHist1Name);
    alphaHist1.SetXTitle("|\\alpha|  [\\circ]");
    


    MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3");
    TH1  &alphaHist2 = alphaAfter->GetHist();
    TString alphaHist2Name ("TestAlphaAfterCuts");
    
    
    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
	      << " fThetaRangeString is NOT defined..." << endl;
	alphaHist2Name += ("ThetaRangeStringUndefined");
    }
    else
    {
	alphaHist2Name += (fThetaRangeString);
    }

    alphaHist2.SetXTitle("|\\alpha|  [\\circ]");
    alphaHist2.SetName(alphaHist2Name);


       
    //fPsFilename  -> NewPage();


    // Canvas is deleted after saving the plot in ps file
    // This is to allow for GRID analysis

    TCanvas *c = new TCanvas("TestAlphaBeforeAfterSC", 
			     "TestAlphaBeforeAfterSC", 
			     600, 600);
    c->Divide(2,2);
    
    

    c->cd(1);
    alphaOFFHistb.DrawCopy();

    c->cd(2);
    alphaOFFHista.DrawCopy();
    
  

    c->cd(3);
    alphaHist1.DrawCopy();

    c->cd(4);
    alphaHist2.DrawCopy();

    c->Modified();
    c-> Update();


    // ********************************************
    // TMP solutions while the TPostScript thing is not working 
    // PsFileName for storing these plots is derived 
    // from PsFilenameString. 
    

    cout << "Alpha distributions will be saved in PostScript file  " ;


    if (!fPsFilenameString.IsNull())
    {
	TString filename = (fPsFilenameString);
	filename += ("AlphaTESTSampleONOFFBeforeAfterCuts.ps");
	cout << filename << endl;	
	c -> SaveAs(filename);
    }

    // END OF TEMPORAL SOLUTION
    // ********************************************
    
    
    // Canvas is deleted after saving the plot in ps file
    // This is to allow for GRID analysis
    delete c;

    //-------------------------------------------
    // fit alpha distribution to get the number of excess events and
    // calculate significance of gamma signal in the alpha plot
  
    const Double_t alphasig = fAlphaSig;
    const Double_t alphamin = fAlphaBkgMin; 
   // alpha min for bkg region in ON data
    
    const Double_t alphamax = fAlphaBkgMax;
//    const Int_t    degree   =    2;
//    const Int_t    degreeOFF   =    2;
    const Bool_t   drawpoly  = kTRUE;
    const Bool_t   fitgauss  = kTRUE;
    const Bool_t   print     = kTRUE;
    
    const Bool_t saveplots   = kTRUE; // Save plots in Psfile
   
     

    if (!ComputeNormFactorTest())
    {
	*fLog << "Normalization factor for test sample (ON-OFF) couold not be computed. Aborting..." << endl;
	return kFALSE;
    }



     // Normalization factor factor is corrected using the estimated 
    // number of excess events and the total number of ON events


 
    if (fTuneNormFactor)
    {


	// Run MHFindSignificanceONOFF::FindSigmaONOFF
	// to get estimation of number of excess gamma events. 
	// This number will be used in the correction of the 
	// normalization factor 
	
	// No plots will be produced this time...

	const Bool_t   drawpolyTest  = kFALSE;
	const Bool_t   fitgaussTest  = kFALSE;
	const Bool_t   printTest     = kFALSE;
	
	const Bool_t saveplotsTest   = kFALSE; // Save plots in Psfile
	
	// TPostScript* DummyPs  = new TPostScript("dummy.ps");

	TString DummyPs = ("Dummy");

	MHFindSignificanceONOFF findsigTest;
	findsigTest.SetRebin(kTRUE);
	findsigTest.SetReduceDegree(kFALSE);
	findsigTest.SetUseFittedQuantities(fUseFittedQuantities);
	

	if (findsigTest.FindSigmaONOFF(&alphaHist2, &alphaOFFHista, 
				       fNormFactorTest, 
        			       alphamin, alphamax, 
				       fDegreeON, fDegreeOFF,
				       alphasig, drawpolyTest, 
				       fitgaussTest, printTest, 
				       saveplotsTest, DummyPs))
	{
	    
	    // Update values of Nex and SigmaLiMa for Test sammple
	
	    fSigmaLiMaTest = double(findsigTest.GetSignificance());

	    if(fUseFittedQuantities)
	    {
		fNexTest = double(findsigTest.GetNexONOFFFitted());
	    }
	    else
	    {
		fNexTest = double(findsigTest.GetNexONOFF());
	    }


	}
	else
	{

	    *fLog << "Normalization Factor tuning for Train sample is not possible"
		  << endl;
	    fSigmaLiMaTest = 0.0;
	    fNexTest = 0.0;
	}
    

//	DummyPs -> Close();
	//DummyPs = NULL;
	//      delete DummyPs;

	
   
	
	
	Int_t EventsInTestMatrixOFF = fMatrixTestOFF->GetM().GetNrows();
	Double_t Ngammas = fNexTest/fGammaEfficiency;
	Double_t GammaFraction = Ngammas/EventsInTestMatrixOFF;
	
	
	*fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
	      << "fNormFactorTest is corrected with fraction of gammas in ON "
	      << "Data Test sample and the number of OFF events before cuts." 
	      << endl
	      << "EventsInTestMatrixOFF = " << EventsInTestMatrixOFF
	      << " Excess events in Test ON sample = " << fNexTest
	      << " fGammaEfficiency = " << fGammaEfficiency << endl
	      << "fNormFactorTest Value before correction is " 
	      << fNormFactorTest << endl;
	
	fNormFactorTest = fNormFactorTest - GammaFraction;

	*fLog << "fNormFactorTest Value after correction is "
	      << fNormFactorTest << endl;
	    

    }
    

    
   
    if (fNormFactorFromAlphaBkg)
    {
	// Normalization factor computed using alpha bkg region
	Double_t NewNormFactor = 
	    ComputeNormFactorFromAlphaBkg(&alphaHist2, &alphaOFFHista, 
					  fAlphaBkgMin, fAlphaBkgMax);
	
	*fLog << "Normalization factor computed from alpha plot (after cuts) " << endl
	      << "using counted number of ON and OFF events in alpha region " << endl
	      << "defined by range " << fAlphaBkgMin << "-" << fAlphaBkgMax << endl
	      << "Normalization factor = " << NewNormFactor << endl;

	*fLog << "Normalization factor used is the one computed in bkg region; " << endl
	      << "i.e. " << NewNormFactor << " instead of " << fNormFactorTest << endl;

	fNormFactorTest = NewNormFactor;


    }



    // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF


    MHFindSignificanceONOFF findsig;
    findsig.SetRebin(kTRUE);
    findsig.SetReduceDegree(kFALSE);
    findsig.SetUseFittedQuantities(fUseFittedQuantities);

    TString psfilename; // Name of psfile where Alpha plot will be stored
    if (!fPsFilenameString.IsNull())
    {
	psfilename += (fPsFilenameString);
	psfilename += ("TESTSample"); 
    }
    else
    {
	psfilename += ("TESTSample");
    }


    findsig.FindSigmaONOFF(&alphaHist2, &alphaOFFHista, 
			   fNormFactorTest, 
			   alphamin, alphamax, 
			   fDegreeON, fDegreeOFF,
			   alphasig, drawpoly, fitgauss, print, 
			   saveplots, psfilename);
    
    const Double_t significance = findsig.GetSignificance();
    const Double_t alphasi = findsig.GetAlphasi();


    // Update values of Nex and SigmaLiMa for Test sammple

    fSigmaLiMaTest = double(findsig.GetSignificance());

    if(fUseFittedQuantities)
    {
	fNexTest = double(findsig.GetNexONOFFFitted());
    }
    else
    {
	fNexTest = double(findsig.GetNexONOFF());
    }


    // Alpha distributions and containers with teh supercuts applied 
    // are written into root file 

    TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE",
		   "Alpha Distributions for several Theta bins");

    alphaHist2.Write();
    alphaOFFHista.Write();
    
    (supapptestON.GetTreePointer())->Write();
    (supapptestOFF.GetTreePointer())->Write();

    *fLog  << "TTree objects containing supercuts applied to TEST sample are "
	   << "written into root file " << fAlphaDistributionsRootFilename << endl; 
 
    (supapptestON.GetTreePointer())->Print();
    (supapptestOFF.GetTreePointer())->Print();

    //rootfile.Close();
    

    // Boolean variable that controls in class MSupercutsCalcONOFF
    // wether the supercuts are stored or not
    // is set to the default value, kFALSE.

    fCalcHadTest -> SetStoreAppliedSupercuts(kFALSE);
    fCalcHadTestOFF -> SetStoreAppliedSupercuts(kFALSE);



    *fLog << "Significance of gamma signal after supercuts in Test sample : "
         << significance << " (for |alpha| < " << alphasi << " degrees)" 
         << endl;
    //-------------------------------------------


    *fLog << "Test of supercuts on TEST sample finished" << endl;
    *fLog << "======================================================" << endl;

    return kTRUE;
}





// Function that applies the optimized supercuts on the TEST sample.

// Supercuts applied (LengthUp, LengthLow...) 
// to all the individual events, as well as 
// the features of the shower images (Length,Width...) are 
// stored in 2 TTree objects (for ON and OFF data) in the 
// root file specified by variable "fAlphaDistributionsRootFilename"


Bool_t MFindSupercutsONOFF::TestParamsOnTrainSample()
{
    if (fMatrixTrain->GetM().GetNrows() <= 0)
    {
        *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; train matrix ON has no entries" 
              << endl;
        return kFALSE;
    }


    if (fMatrixTrainOFF->GetM().GetNrows() <= 0)
    {
        *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; train matrix OFF has no entries" 
              << endl;
        return kFALSE;
    }

    // -------------   TEST the supercuts on TRAIN sample  --------------------------
    //
    *fLog << "Producing output of the SUPERCUTS ON-OFF optimization " << endl;

    // -----------------------------------------------------------------
    // read optimum parameter values from file filenameParam
    // into array 'supercutsPar'

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

    *fLog << "Optimum parameter values for supercuts were read from file '";
    *fLog << fFilenameParam << "' :" << endl;

    const TArrayD &supercutsPar = scin.GetParameters();
    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
        *fLog << supercutsPar[i] << ",  ";
    *fLog << endl;
    //---------------------------


    // -----------------------------------------------------------------
    if (fHadronnessName.IsNull())
    {
        *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; hadronness name for ON data is not defined... aborting";
        *fLog << endl;
        return kFALSE;
    }

    if (fHadronnessNameOFF.IsNull())
    {
        *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; hadronness name for OFF data is not defined... aborting";
        *fLog << endl;
        return kFALSE;
    }


    if(fTrainONSupercutsAppliedName.IsNull())
    {
        *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; " 
	      << " MTSupercutsApplied Tree name for ON Train data is not defined... aborting"
	      << endl;
        return kFALSE;
    }

    if(fTrainOFFSupercutsAppliedName.IsNull())
    {
        *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; " 
	      << " MTSupercutsApplied Tree name for OFF Train data is not defined... aborting"
	      << endl;
        return kFALSE;
    }
    
   
    MSupercuts supercuts;
    supercuts.SetParameters(supercutsPar);

    // Objects containing the trees used to store supercuts applied
    // TTree Branches are also created
    MTSupercutsApplied supapptrainON(fTrainONSupercutsAppliedName,"");
    supapptrainON.CreateTreeBranches();

    MTSupercutsApplied supapptrainOFF(fTrainOFFSupercutsAppliedName,"");
    supapptrainOFF.CreateTreeBranches();

    fCalcHadTrainOFF -> SetHadronnessName(fHadronnessNameOFF);
    fCalcHadTrain -> SetHadronnessName(fHadronnessName);


    // Settings related to the storage of teh supercuts applied 
    fCalcHadTrain -> SetSupercutsAppliedName(fTrainONSupercutsAppliedName);
    fCalcHadTrain -> SetStoreAppliedSupercuts(kTRUE);

    fCalcHadTrainOFF -> SetSupercutsAppliedName(fTrainOFFSupercutsAppliedName);
    fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kTRUE);
    

    // Set boolean variable that controls wether cuts are 
    // dynamic or static

    fCalcHadTrain -> SetVariableUseStaticCuts(fUseStaticCuts);
    fCalcHadTrainOFF -> SetVariableUseStaticCuts(fUseStaticCuts);


    if (!fUseStaticCuts)
    {
	// Set boolean variable that controls wether the theta variable 
	// is used or not in the computation of the dynamical cuts
	
	fCalcHadTrain -> SetVariableNotUseTheta(fNotUseTheta);
	fCalcHadTrainOFF -> SetVariableNotUseTheta(fNotUseTheta);
    }


    // apply the supercuts to OFF data

    MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5");
    MContinue applysupercutsOFF(&scfilterOFF);

     // apply the supercuts to ON data
    MF scfilter(fHadronnessName+".fHadronness>0.5");
    MContinue applysupercuts(&scfilter);

    


    // ****************************************************
    // Filling OFF alpha distribution 

    MParList parlistOFF;
    MTaskList tasklistOFF;

    MMatrixLoop loopOFF(fMatrixTrainOFF);

    // plot alpha before applying the supercuts
    const TString  mh3NameOFFB = "AlphaOFFBefore";
    MBinning binsalphaOFFbef("Binning"+mh3NameOFFB);
    //binsalphaOFFbef.SetEdges(54, -12.0, 96.0);
    binsalphaOFFbef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);


    // |alpha| is assumed to be in column 7 of the matrix
    MH3 alphaOFFbefore("MatrixTrainOFF[7]");
    alphaOFFbefore.SetName(mh3NameOFFB);

    
    MFillH fillalphaOFFbefore(&alphaOFFbefore);
    fillalphaOFFbefore.SetName("FillAlphaOFFBefore");

   

    // fill OFF alpha after applying the supercuts
    const TString  mh3NameOFFA = "AlphaOFFAfter";
    MBinning binsalphaOFFaft("Binning"+mh3NameOFFA);
    //binsalphaOFFaft.SetEdges(54, -12.0, 96.0);
    binsalphaOFFaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);

    MH3 alphaOFFafter("MatrixTrainOFF[7]");
    alphaOFFafter.SetName(mh3NameOFFA);

    //TH1 &alphaOFFhista = alphaOFFafter.GetHist();
    //alphaOFFhista.SetName("alphaOFFAfter-OptimizationParamOutput");

    MFillH fillalphaOFFafter(&alphaOFFafter);
    fillalphaOFFafter.SetName("FillAlphaOFFAfter");

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

    parlistOFF.AddToList(&tasklistOFF);

    parlistOFF.AddToList(&supercuts);

    parlistOFF.AddToList(&supapptrainOFF);

    parlistOFF.AddToList(fCam);
    parlistOFF.AddToList(fMatrixTrainOFF);

    parlistOFF.AddToList(&binsalphaOFFbef);
    parlistOFF.AddToList(&alphaOFFbefore);

    parlistOFF.AddToList(&binsalphaOFFaft);
    parlistOFF.AddToList(&alphaOFFafter);


    //******************************
    // Entries in MtaskList
    tasklistOFF.AddToList(&loopOFF);
    tasklistOFF.AddToList(&fillalphaOFFbefore);

    tasklistOFF.AddToList(fCalcHadTrainOFF);
    tasklistOFF.AddToList(&applysupercutsOFF);

    tasklistOFF.AddToList(&fillalphaOFFafter);

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

    MProgressBar barOFF;
    MEvtLoop evtloopOFF;
    evtloopOFF.SetParList(&parlistOFF);
    evtloopOFF.SetName("EvtLoopOptimizationParamOutputOFF");
    evtloopOFF.SetProgressBar(&barOFF);
    Int_t maxeventsOFF = -1;
    if (!evtloopOFF.Eventloop(maxeventsOFF))
        return kFALSE;

    tasklistOFF.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // draw the alpha plots

    MH3* alphaOFFBefore = (MH3*)parlistOFF.FindObject(mh3NameOFFB, "MH3");
    TH1  &alphaOFFHistb = alphaOFFBefore->GetHist();
    alphaOFFHistb.SetXTitle("|\\alpha|  [\\circ]");

    TString alphaOFFHistbName ("TrainAlphaOFFBeforeCuts");
    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
	      << " fThetaRangeString is NOT defined..." << endl;
	alphaOFFHistbName += ("ThetaRangeStringUndefined");
    }
    else
    {
	alphaOFFHistbName += (fThetaRangeString);
    }

    alphaOFFHistb.SetName(alphaOFFHistbName);

    MH3* alphaOFFAfter = (MH3*)parlistOFF.FindObject(mh3NameOFFA, "MH3");
    TH1  &alphaOFFHista = alphaOFFAfter->GetHist();
    TString alphaOFFHistaName ("TrainAlphaOFFAfterCuts");
    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
	      << " fThetaRangeString is NOT defined..." << endl;
	alphaOFFHistaName += ("ThetaRangeStringUndefined");
    }
    else
    {
	alphaOFFHistaName += (fThetaRangeString);
    }


    alphaOFFHista.SetXTitle("|\\alpha|  [\\circ]");
    alphaOFFHista.SetName(alphaOFFHistaName);


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



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


    MParList  parlist2;
    MTaskList tasklist2;

    MMatrixLoop loop(fMatrixTrain);

    // plot alpha before applying the supercuts
    const TString  mh3NameB = "AlphaBefore";
    MBinning binsalphabef("Binning"+mh3NameB);
    //binsalphabef.SetEdges(54, -12.0, 96.0);
    binsalphabef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
    
    // |alpha| is assumed to be in column 7 of the matrix
    MH3 alphabefore("MatrixTrain[7]");
    alphabefore.SetName(mh3NameB);

    MFillH fillalphabefore(&alphabefore);
    fillalphabefore.SetName("FillAlphaBefore");


    // plot alpha after applying the supercuts
    const TString  mh3NameA = "AlphaAfter";
    MBinning binsalphaaft("Binning"+mh3NameA);
    //binsalphaaft.SetEdges(54, -12.0, 96.0);
    binsalphaaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);

    MH3 alphaafter("MatrixTrain[7]");
    alphaafter.SetName(mh3NameA);

    
    MFillH fillalphaafter(&alphaafter);
    fillalphaafter.SetName("FillAlphaAfter");

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

    parlist2.AddToList(&tasklist2);

    parlist2.AddToList(&supercuts);

    parlist2.AddToList(&supapptrainON);

    parlist2.AddToList(fCam);
    parlist2.AddToList(fMatrixTrain);

    parlist2.AddToList(&binsalphabef);
    parlist2.AddToList(&alphabefore);

    parlist2.AddToList(&binsalphaaft);
    parlist2.AddToList(&alphaafter);

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

    tasklist2.AddToList(&loop);
    tasklist2.AddToList(&fillalphabefore);

    tasklist2.AddToList(fCalcHadTrain);
    tasklist2.AddToList(&applysupercuts);

    tasklist2.AddToList(&fillalphaafter);

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

    MProgressBar bar2;
    MEvtLoop evtloop2;
    evtloop2.SetParList(&parlist2);
    evtloop2.SetName("EvtLoopOptimizationParamOutput");
    evtloop2.SetProgressBar(&bar2);
    Int_t maxevents2 = -1;
    if (!evtloop2.Eventloop(maxevents2))
        return kFALSE;

    tasklist2.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // draw the alpha plots

    MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3");
    TH1  &alphaHist1 = alphaBefore->GetHist();
    TString alphaHist1Name ("TrainAlphaBeforeCuts");
    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
	      << " fThetaRangeString is NOT defined..." << endl;
	alphaHist1Name += ("ThetaRangeStringUndefined");
    }
    else
    {
	alphaHist1Name += (fThetaRangeString);
    }


    alphaHist1.SetXTitle("|\\alpha|  [\\circ]");
    alphaHist1.SetName(alphaHist1Name);

    MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3");
    TH1  &alphaHist2 = alphaAfter->GetHist();
    TString alphaHist2Name ("TrainAlphaAfterCuts");
    if (fThetaRangeString.IsNull())
    {
	*fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
	      << " fThetaRangeString is NOT defined..." << endl;
	alphaHist2Name += ("ThetaRangeStringUndefined");
    }
    else
    {
	alphaHist2Name += (fThetaRangeString);
    }

    alphaHist2.SetXTitle("|\\alpha|  [\\circ]");
    alphaHist2.SetName(alphaHist2Name);


    // fPsFilename -> NewPage();

    
    // Canvas is deleted after saving the plot in ps file
    // This is to allow for GRID analysis

    TCanvas *c = new TCanvas("TrainAlphaBeforeAfterSC", 
			     "TrainAlphaBeforeAfterSC", 
			     600, 600);

    c->Divide(2,2);



    c->cd(1);
    alphaOFFHistb.DrawCopy();

    c->cd(2);
    alphaOFFHista.DrawCopy();


    //c->Modified();
    //c-> Update();
    

    // fPsFilename -> NewPage();

    c->cd(3);
    alphaHist1.DrawCopy();

    c->cd(4);
    alphaHist2.DrawCopy();


    c -> Modified();
    c -> Update();
    
     // ********************************************
    // TMP solutions while the TPostScript thing is not working 
    // PsFileName for storing these plots is derived 
    // from PsFilenameString. 

    cout << "Alpha distributions will be saved in PostScript file  " ;

    if (!fPsFilenameString.IsNull())
    {
	TString filename = (fPsFilenameString);
	filename += ("AlphaTRAINSampleONOFFBeforeAfterCuts.ps");
	cout << filename << endl;
	c -> SaveAs(filename);
    }

    // END OF TEMPORAL SOLUTION
    // ********************************************
	    



    // fPsFilename -> NewPage();


    // Canvas is deleted after saving the plot in ps file
    // This is to allow for GRID analysis

    delete c;




    //-------------------------------------------
    // fit alpha distribution to get the number of excess events and
    // calculate significance of gamma signal in the alpha plot
  
    const Double_t alphasig = fAlphaSig;
    const Double_t alphamin = fAlphaBkgMin; 
   // alpha min for bkg region in ON data
    
    const Double_t alphamax = fAlphaBkgMax;
//    const Int_t    degree   =    2;
//    const Int_t    degreeOFF   =    2;
    const Bool_t   drawpoly  = kTRUE;
    const Bool_t   fitgauss  = kTRUE;
    const Bool_t   print     = kTRUE;
    
    Bool_t saveplots   = kTRUE; // Save plots in Psfile
  
    

    
    

    if (!ComputeNormFactorTrain())
    {
	*fLog << "Normalization factor for train sample (ON-OFF) could not be computed. Aborting..." << endl;
	return kFALSE;
    }
	    

    // Normalization factor factor is corrected using the estimated 
    // number of excess events and the total number of OFF events
    // fNormFactor = fNormFactor - Ngammas/EventsInMatrixOFF

 
    if (fTuneNormFactor)
    {


	// Run MHFindSignificanceONOFF::FindSigmaONOFF
	// to get estimation of number of excess gamma events. 
	// This number will be used in the correction of the 
	// normalization factor 
	
	// No plots will be produced this time...

	const Bool_t   drawpolyTest  = kFALSE;
	const Bool_t   fitgaussTest  = kFALSE;
	const Bool_t   printTest     = kFALSE;
	
	const Bool_t saveplotsTest   = kFALSE; // Save plots in Psfile
	
	//TPostScript* DummyPs  = new TPostScript("dummy.ps");

	TString DummyPs = ("Dummy");

	MHFindSignificanceONOFF findsigTest;
	findsigTest.SetRebin(kTRUE);
	findsigTest.SetReduceDegree(kFALSE);
	findsigTest.SetUseFittedQuantities(fUseFittedQuantities);

	if(findsigTest.FindSigmaONOFF(&alphaHist2, &alphaOFFHista, 
				      fNormFactorTrain, 
				      alphamin, alphamax, 
				      fDegreeON, fDegreeOFF,
				      alphasig, drawpolyTest, 
				      fitgaussTest, printTest, 
				      saveplotsTest, DummyPs))

	{

	    // Update values of Nex and SigmaLiMa for Test sammple
	
	    fSigmaLiMaTrain = double(findsigTest.GetSignificance());

	    if(fUseFittedQuantities)
	    {
		fNexTrain = double(findsigTest.GetNexONOFFFitted());
	    }
	    else
	    {
		fNexTrain = double(findsigTest.GetNexONOFF());
	    }
	    
	}

	else
	{
	    *fLog << "Normalization Factor tuning for Train sample is not possible"
		  << endl;

	    fSigmaLiMaTrain = 0.0;
	    fNexTrain = 0.0;

	}

	//DummyPs -> Close();
	//DummyPs = NULL;
	//delete DummyPs;

	


	Int_t EventsInTrainMatrixOFF = fMatrixTrainOFF->GetM().GetNrows();
	Double_t Ngammas = fNexTrain/fGammaEfficiency;
	Double_t GammaFraction = Ngammas/EventsInTrainMatrixOFF;
	

	*fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
	      << "fNormFactorTrain is corrected with fraction of gammas in ON "
	      << "Data Train sample and the number of OFF events before cuts." 
	      << endl
	      << "EventsInTrainMatrixOFF = " << EventsInTrainMatrixOFF
	      << " Excess events in Train ON sample = " << fNexTrain
	      << " fGammaEfficiency = " << fGammaEfficiency << endl
	      << "fNormFactorTrain Value before correction is " 
	      << fNormFactorTrain << endl;

	fNormFactorTrain = fNormFactorTrain - GammaFraction;

	*fLog << "fNormFactorTrain Value after correction is "
	      << fNormFactorTrain << endl;
       
	    

    }

    if (fNormFactorFromAlphaBkg)
    {
	Double_t NewNormFactor = 
	    ComputeNormFactorFromAlphaBkg(&alphaHist2, &alphaOFFHista, 
					  fAlphaBkgMin, fAlphaBkgMax);
	
	*fLog << "Normalization factor computed from alpha plot (after cuts) " << endl
	      << "using counted number of ON and OFF events in alpha region " << endl
	      << "defined by range " << fAlphaBkgMin << "-" << fAlphaBkgMax << endl
	      << "Normalization factor = " << NewNormFactor << endl;

	*fLog << "Normalization factor used is the one computed in bkg region; " << endl
	      << "i.e. " << NewNormFactor << " instead of " << fNormFactorTrain << endl;

	fNormFactorTrain = NewNormFactor;


    }


    MHFindSignificanceONOFF findsig;
    findsig.SetRebin(kTRUE);
    findsig.SetReduceDegree(kFALSE);
    findsig.SetUseFittedQuantities(fUseFittedQuantities);

    TString psfilename; // Name of psfile where Alpha plot will be stored
    if (!fPsFilenameString.IsNull())
    {
	psfilename += (fPsFilenameString);
	psfilename += ("TRAINSample"); 
    }
    else
    {
	psfilename += ("TRAINSample");
    }


    findsig.FindSigmaONOFF(&alphaHist2, &alphaOFFHista, 
			   fNormFactorTrain, 
			   alphamin, alphamax, 
			   fDegreeON, fDegreeOFF,
			   alphasig, drawpoly, fitgauss, print, 
			   saveplots, psfilename);
    
    const Double_t significance = findsig.GetSignificance();
    const Double_t alphasi = findsig.GetAlphasi();

    // Update Train values Nex, SigmaLiMa

  fSigmaLiMaTrain = double(findsig.GetSignificance());

  if(fUseFittedQuantities)
  {
      fNexTrain = double(findsig.GetNexONOFFFitted()); 
  }
  else
  {
      fNexTrain = double(findsig.GetNexONOFF()); 
  }
  

    // Alpha distributions and containers with the supercuts applied 
    // are written into root file 

    TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE",
		   "Alpha Distributions for several Theta bins");

    alphaHist2.Write();
    alphaOFFHista.Write();
    
    (supapptrainON.GetTreePointer())->Write();
    (supapptrainOFF.GetTreePointer())->Write();

    *fLog  << "TTree objects containing supercuts applied to TEST sample are "
	   << "written into root file " << fAlphaDistributionsRootFilename << endl; 
 
    (supapptrainON.GetTreePointer())->Print();
    (supapptrainOFF.GetTreePointer())->Print();


    //rootfile.Close();
    

    // Boolean variable that controls in class MSupercutsCalcONOFF
    // wether the supercuts are stored or not
    // is set to the default value, kFALSE.

    fCalcHadTrain -> SetStoreAppliedSupercuts(kFALSE);
    fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kFALSE);
    




    *fLog << "Significance of gamma signal after supercuts in train sample : "
         << significance << " (for |alpha| < " << alphasi << " degrees)" 
         << endl;
    //-------------------------------------------


    *fLog << "Test of supercuts on TRAIN sample finished" << endl;
    *fLog << "======================================================" << endl;

    return kTRUE;
}





















