/* ======================================================================== *\
!
! *
! * 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): David Paneque   12/2003 <mailto:dpaneque@mppmu.mpg.de>
!              
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MFindSupercutsONOFFThetaLoop                                             //
//                                                                         //
// Class for optimizing the parameters of the supercuts                     //
// Using ON and OFF data runing over data which are subdivided in 
// several Theta ranges.                                                  //
//                                                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////


#include "MFindSupercutsONOFFThetaLoop.h"


#include <math.h>            // fabs 

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

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include "MFindSupercutsONOFF.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" // only for magic
#include "MFRandomSplit.h"
#include "MH3.h"
#include "MHCT1Supercuts.h"
#include "MHFindSignificance.h" // To be removed at some point...
#include "MHFindSignificanceONOFF.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;

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


    // Cuts (low and up) in variable ThetaOrig.fVal
    // The default is not cut, i.e. all values (0-1) are taken
    // fThetaOrig.fVal is measured in Radians; thus 1 = 57 degrees.

    fThetaMin = 0.0;
    fThetaMax = 1.0;

    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;

    fActualCosThetaBinCenter = 0;
    
    fNormFactorTrainHist = NULL;
    //fNormFactorTrainHist -> SetDirectory(0);
    fNormFactorTestHist = NULL;
    //fNormFactorTestHist -> SetDirectory(0);
    fSigmaLiMaTrainHist = NULL;
    //fSigmaLiMaTrainHist -> SetDirectory(0);
    fSigmaLiMaTestHist = NULL;
    //fSigmaLiMaTestHist -> SetDirectory(0);
    fNexTrainHist = NULL;
    //fSigmaLiMaTestHist -> SetDirectory(0);
    fNexTestHist = NULL;
    //fNexTestHist -> SetDirectory(0);


    fSuccessfulThetaBinsHist = NULL;

    /*
    fNexVSAlphaSigTrain = NULL;
    fNexVSAlphaSigTest = NULL;
    
    fSignificanceVSAlphaSigTrain = NULL;
    fSignificanceVSAlphaSigTest = NULL;
    */


    fNEvtsInTrainMatrixONHist = NULL; 
    fNEvtsInTestMatrixONHist = NULL; 
    fNEvtsInTrainMatrixOFFHist = NULL; 
    fNEvtsInTestMatrixOFFHist = NULL; 


    fTrainMatrixONFilenameVector = NULL;
    fTestMatrixONFilenameVector= NULL;
    fTrainMatrixOFFFilenameVector= NULL;
    fTestMatrixOFFFilenameVector= NULL;

    
    fOptSCParamFilenameVector = NULL;

    fThetaRangeStringVector = NULL;

    fPsFilename = NULL;



    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
    
    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

    // boolean variables are set to kFALSE as default 

    fReadMatricesFromFile = kFALSE;
    fOptimizeParameters = kFALSE;
    fTestParameters = kFALSE;

    // Fraction of Train/Test ON/OFF samples is set to 0.5 as default

   fWhichFractionTrain = 0.5;  // number <= 1; specifying fraction of ON Train events
   fWhichFractionTest = 0.5;  // number <= 1; specifying fraction of ON Test events

   fWhichFractionTrainOFF = 0.5;  // number <= 1; specifying fraction of OFF  Train events
   fWhichFractionTestOFF = 0.5;  // number <= 1; specifying fraction of OFF Test events



    // 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
    fUseFittedQuantities = kTRUE;



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



    // fInitSCPar  and fInitSCParSteps initialized to zero
    
    fInitSCPar.Set(0);

    fInitSCParSteps.Set(0);
        
}

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

    *fLog << "Destructor of MFindSupercutsONOFFThetaLoop is called" 
	  << endl;

    fPsFilename = NULL;

    delete fNormFactorTrainHist;
    delete fNormFactorTestHist;
    delete fSigmaLiMaTrainHist;
    delete fSigmaLiMaTestHist;
    delete fNexTrainHist;
    delete fNexTestHist;


    
    delete fSuccessfulThetaBinsHist;
    

    delete fNEvtsInTrainMatrixONHist; 
    delete fNEvtsInTestMatrixONHist; 
    delete fNEvtsInTrainMatrixOFFHist; 
    delete fNEvtsInTestMatrixOFFHist; 

    delete [] fTrainMatrixONFilenameVector;
    delete [] fTestMatrixONFilenameVector;
    delete [] fTrainMatrixOFFFilenameVector;
    delete [] fTestMatrixOFFFilenameVector;
    
    delete [] fOptSCParamFilenameVector;
    
    delete [] fThetaRangeStringVector;


    *fLog << "Destructor of MFindSupercutsONOFFThetaLoop 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 MFindSupercutsONOFFThetaLoop::SetPostScriptFile (TPostScript* PsFile)
{
    fPsFilename = PsFile;

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

}


Bool_t MFindSupercutsONOFFThetaLoop::SetGammaEfficiency (Double_t gammaeff)
{
    if (gammaeff <= 0.0 || gammaeff > 1)
    {
	*fLog << " ******   ERROR  *********"  << endl
 	      << "MFindSupercutsONOFFThetaLoop :: SetGammaEfficiency; " << endl
	      << "The requested Gamma efficiency " << gammaeff << " can not occur" << endl
	      << "Member data fGammaEfficiency can not be set properly." << endl
	      << "Default value for fGammaEfficiency, which is " 
	      << fGammaEfficiency << " will remain." << endl;


	return kFALSE;
    }

    fGammaEfficiency = gammaeff;

    return kTRUE;

}


Bool_t MFindSupercutsONOFFThetaLoop::ReadSCParamsFromAsciiFile(const char* filename,
								  Int_t Nparams)
{
    const Int_t NColumns = 2; 
    
    Double_t tmp_vector[NColumns];
    Int_t counter = 0;
    Double_t tmpdouble = 0.0;

    
    // File is read and data stored in data_vector
    
    ifstream infile (filename, ios::in);
    
    if ( !infile )
    {
	cout << "Input file could not be opened" << endl; 	
	return kFALSE;
      //	exit (1);
    }
    
    counter = 0;
    while (infile >> tmp_vector [0])
    {
	for (Int_t i = 1; i < NColumns; i++)
	{
	    infile >> tmp_vector[i];
	}
	counter++;
    }
    
    infile.close();
    
    
    // Length of TArray vectors is set to the length of the 
    // vectors defined in ascii file.
    // If such length does not match with the specified 
    // by variable Nparams, the function complains and 
    // aborts.
    
    if (counter != Nparams)
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::ReadSCParamsFromAsciiFile" << endl
	      << "Length of vectors in ascii file " << filename << "  is  " << counter 
	      << " , which does not match with the expected length " << Nparams << endl  
	      << " fInitSCPar and fInitSCParSteps are NOT retrieved from ascii file" << endl;

	return kFALSE;
    }
    
    fInitSCPar.Set(counter);
    fInitSCParSteps.Set(counter);
    
    
  // Vectors are filled

  ifstream infile2 (filename, ios::in);
  
  if ( !infile2 )
    {
      cout << "Input file could not be opened" << endl; 
      return kFALSE;
      // exit (1);
    }
  
  for (Int_t i = 0; i < fInitSCPar.GetSize(); i++)
    {

	infile2 >> tmpdouble;
	fInitSCPar[i] = tmpdouble;
	
	infile2 >> tmpdouble;
	fInitSCParSteps[i] = tmpdouble;
	
    }


  // Print the vectors read 

  *fLog << "*****  INITIAL SC PARAMETERS AND STEPS TAKEN FROM ASCII FILE ********* " << endl;
  for (Int_t i = 0; i < fInitSCPar.GetSize(); i++)
    {
	cout << fInitSCPar[i] << setw(20) <<  fInitSCParSteps[i] << endl;
    }

  // 
  
 // File close 
  infile2.close();
  
  return kTRUE;
}


Bool_t MFindSupercutsONOFFThetaLoop::SetInitSCPar (TArrayD &d)
{
    *fLog << "MFindSupercutsONOFFThetaLoop; Initial SC parameters are "
	  << "set using function SetInitSCPar (TArrayD &d)" << endl;

    fInitSCPar = d;

    return kTRUE;
}


Bool_t MFindSupercutsONOFFThetaLoop::SetInitSCParSteps (TArrayD &d)
{

    *fLog << "MFindSupercutsONOFFThetaLoop; Initial SC parameters STEPS are "
	  << "set using function SetInitSCParSteps (TArrayD &d)" << endl;

    fInitSCParSteps = d;

    return kTRUE;
}



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

    fAlphaSig = alphasig;

    return kTRUE;
}

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

    fAlphaBkgMin = alpha;

    return kTRUE;
}


Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaBkgMax(Double_t alpha)
{
    // check that alpha is within the limits 0-90
    if (alpha <= 0 || alpha > 90.001)
    {
	*fLog << "MFindSupercutsONOFFThetaLoop ::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 MFindSupercutsONOFFThetaLoop::CheckAlphaSigBkg()
{

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

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

	return kFALSE;
    }

    return kTRUE;

}




// Function that sets the value of fCosThetaRangeVector and 
// fCosThetaBinCenterVector

Bool_t MFindSupercutsONOFFThetaLoop::SetCosThetaRangeVector(
    const TArrayD &d)
{
    if (d.GetSize() < 2)
    {
        *fLog << err 
	      << "MFindSupercutsONOFFThetaLoop::SetCosThetaRangeVector: " 
	      << endl
	      << "Size of d is smaller than 2... "
	      << "fCosThetaRangeVector can not be defined properly"
	      << endl;
        return kFALSE;
    }

    fCosThetaRangeVector.Set(d.GetSize());
    fCosThetaRangeVector = d;


    // Give a "soft warning" if Vector has increasing values 
    // of CosThetas; which means decreasing values of Thetas
    /*
    if (fCosThetaRangeVector[0] <= fCosThetaRangeVector[1])
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::SetCosThetaRangeVector; " 
	      << endl
	      << "Values in vector fCosThetaRangeVector are in increasing "
	      << "order, i.e., values in theta will be in decreasing order"
	      << endl
	      << "DO NOT forget it !!" << endl;
    }

    */

    // fCosThetaBinCenterVector is defined

    Int_t dim = fCosThetaRangeVector.GetSize() - 1;
    fCosThetaBinCenterVector.Set(dim);

    for (Int_t i = 0; i < dim; i++)
    {
	fCosThetaBinCenterVector[i] = 
	    (fCosThetaRangeVector[i+1]+fCosThetaRangeVector[i])/2;
    }

    return kTRUE;
}

// Function that sets the range of Theta (fThetaMin, fThetaMax) 
// where data will be used in the iteration 
// thetabin. This theta range is used when geting data from 
// file and filling matrices. If matrices are already built 
// and just read from root files, this theta range is useless. 


// **** IMPORTANT *********
// thetabin is an integer  value that goes from 1 to 
// fCosThetaBinCenterVector.GetSize()
// ************************

Bool_t MFindSupercutsONOFFThetaLoop::SetThetaRange (Int_t thetabin)
{
    // Check that iteration thetabin is smaller than 
    // components of fCosThetaBinCenterVector
   
    if (fCosThetaRangeVector.GetSize() < 2 ||
	fCosThetaBinCenterVector.GetSize() < 1)
	
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl
	      << "Dimension of vectors fCosThetaRangeVector or/and " 
	      << "fCosThetaBinCenterVector are smaller than 2 and 1 " 
	      << "respectively. Range for theta can not be defined." << endl;
	
	return kFALSE;
    }
    
    if (thetabin < 1 || 
	thetabin > fCosThetaBinCenterVector.GetSize())
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl
	      << "thetabin value is " << thetabin 
	      << " Is is outside the possible values given by " 
	      << "vector fCosThetaRangeVector. "
	      << "Range for theta can not be defined." << endl;
	
	return kFALSE;
    }
    
    
    fThetaMin = fCosThetaRangeVector[thetabin-1];
    fThetaMax = fCosThetaRangeVector[thetabin];
    

    // Now costheta values have to be converted to Radians
    // Function TMath::ACos(x) gives angle in Radians

    
    fThetaMin = TMath::ACos(fThetaMin);
    fThetaMax = TMath::ACos(fThetaMax);

    // If fThetaMin > fThetaMax means that fCosThetaRangeVector
    // is defined with increasing values of CosTheta, which 
    // means decreasing values of Theta

    // In such case, values of fThetaMin and fThetaMax are swapped

    if (fThetaMin > fThetaMax)
    {
	Double_t tmp = fThetaMin;
	fThetaMin = fThetaMax;
	fThetaMax = tmp;
    }

    // Info for user
    *fLog << "-------------------------------------------------"<< endl
	  << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl
	  << "Theta range for this iteration is set to "
	  << fThetaMin << "-" << fThetaMax << " Radians." << endl 
	  << "-------------------------------------------------"<< endl;
    
    return kTRUE;


}


Bool_t MFindSupercutsONOFFThetaLoop::SetNormFactorTrainHist()
{
    if (fCosThetaRangeVector.GetSize() <= 1)
    {
        *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNormFactorTrainHist: " 
	      << endl
	      << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
	      << endl;
        return kFALSE;
    }

    
    Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
    Double_t* BinsVector = fCosThetaRangeVector.GetArray();

    fNormFactorTrainHist = new TH1F ("NormFactorTrainHist",
				     "NormFactorTrainHist",
				     HistoNbins,BinsVector );

    return kTRUE;

}


Bool_t MFindSupercutsONOFFThetaLoop::SetNormFactorTestHist()
{
    if (fCosThetaRangeVector.GetSize() <= 1)
    {
        *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNormFactorTestHist: " 
	      << endl
	      << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
	      << endl;
        return kFALSE;
    }


    Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
    Double_t* BinsVector = fCosThetaRangeVector.GetArray();

    fNormFactorTestHist = new TH1F ("NormFactorTestHist",
				     "NormFactorTestHist",
				     HistoNbins, BinsVector);

    return kTRUE;
}


Bool_t MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTrainHist()
{
    if (fCosThetaRangeVector.GetSize() <= 1)
    {
        *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTrainHist: " 
	      << endl
	      << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
	      << endl;
        return kFALSE;
    }


    // At some point, definition of bins should be able to 
    // contain variable bins...

    Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
    
    Double_t* BinsVector = fCosThetaRangeVector.GetArray();

    fSigmaLiMaTrainHist = new TH1F ("SigmaLiMaTrainHist",
				     "SigmaLiMaTrainHist",
				     HistoNbins, BinsVector);

    return kTRUE;
}



Bool_t MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTestHist()
{
    if (fCosThetaRangeVector.GetSize() <= 1)
    {
        *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSigmaLiMaTestHist: " 
	      << endl
	      << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
	      << endl;
        return kFALSE;
    }


   
    Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
    
    Double_t* BinsVector = fCosThetaRangeVector.GetArray();
    fSigmaLiMaTestHist = new TH1F ("SigmaLiMaTestHist",
				     "SigmaLiMaTestHist",
				     HistoNbins, BinsVector);

    return kTRUE;

}


Bool_t MFindSupercutsONOFFThetaLoop::SetNexTrainHist()
{
    if (fCosThetaRangeVector.GetSize() <= 1)
    {
        *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNexTrainHist: " 
	      << endl
	      << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
	      << endl;
        return kFALSE;
    }


    Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
    
    Double_t* BinsVector = fCosThetaRangeVector.GetArray();

    fNexTrainHist = new TH1F ("NexTrainHist",
			      "NexTrainHist",
			      HistoNbins, BinsVector);

    return kTRUE;

}


Bool_t MFindSupercutsONOFFThetaLoop::SetNexTestHist()
{
    if (fCosThetaRangeVector.GetSize() <= 1)
    {
        *fLog << err << "MFindSupercutsONOFFThetaLoop::SetNexTestHist: " 
	      << endl
	      << "Size of fCosThetaRangeVector is <=1 ... fCosThetaRangeVector must exist and have, at least, 2 components defining a single bin"
	      << endl;
        return kFALSE;
    }


    Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
    Double_t* BinsVector = fCosThetaRangeVector.GetArray();

    fNexTestHist = new TH1F ("NexTestHist",
			     "NexTestHist",
			     HistoNbins, BinsVector);

    return kTRUE;

}

Bool_t MFindSupercutsONOFFThetaLoop::SetSuccessfulThetaBinsHist()
{
    if (fCosThetaRangeVector.GetSize() <= 1)
    {
        *fLog << err << "MFindSupercutsONOFFThetaLoop::SetSuccessfulThetaBinsHist(): " 
	      << endl
	      << "Size of fCosThetaRangeVector is <=1 ... "
	      << "fCosThetaRangeVector must exist and have >= 2 components (defining a single bin)"
	      << endl;
        return kFALSE;
    }


    Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
    Double_t* BinsVector = fCosThetaRangeVector.GetArray();

    fSuccessfulThetaBinsHist = new TH1F ("SuccessfulThetaBinsHist",
					 "SuccessfulThetaBinsHist",
					 HistoNbins, BinsVector);

    return kTRUE;

}




Bool_t MFindSupercutsONOFFThetaLoop::SetNexSigmaLiMaNormFactorNEvtsTrainTestHist()
{
    if (fCosThetaRangeVector.GetSize() <= 1)
    {
        *fLog << err 
	      << "MFindSupercutsONOFFThetaLoop::SetNexSigmaLiMaNormFactorTrainTestHist: " 
	      << endl
	      << "Size of fCosThetaRangeVector is <=1 ... " 
	      << "fCosThetaRangeVector must exist and have, at least, "
	      << "2 components defining a single bin"
	      << endl;
        return kFALSE;
    }


    Int_t HistoNbins = fCosThetaRangeVector.GetSize() - 1;
    
    Double_t* BinsVector = fCosThetaRangeVector.GetArray();

    char* x_axis_title = {"Cos(Theta)"};

    fNormFactorTrainHist = new TH1F ("NormFactorTrainHist",
				     "NormFactorTrainHist",
				     HistoNbins, BinsVector);

    
    fNormFactorTrainHist -> SetXTitle(x_axis_title);
    fNormFactorTrainHist -> SetYTitle("Normalization Factor (Non/Noff)");
    

    fNormFactorTestHist = new TH1F ("NormFactorTestHist",
				     "NormFactorTestHist",
				     HistoNbins, BinsVector);

    fNormFactorTestHist -> SetXTitle(x_axis_title);
    fNormFactorTestHist -> SetYTitle("Normalization Factor (Non/Noff)");

    fSigmaLiMaTrainHist = new TH1F ("SigmaLiMaTrainHist",
				     "SigmaLiMaTrainHist",
				     HistoNbins, BinsVector);

    fSigmaLiMaTrainHist -> SetXTitle(x_axis_title);
    fSigmaLiMaTrainHist -> SetYTitle("Significance");


    fSigmaLiMaTestHist = new TH1F ("SigmaLiMaTestHist",
				   "SigmaLiMaTestHist",
				   HistoNbins, BinsVector);

    fSigmaLiMaTestHist -> SetXTitle(x_axis_title);
    fSigmaLiMaTestHist -> SetYTitle("Significance");



    fNexTrainHist = new TH1F ("NexTrainHist",
			      "NexTrainHist",
			      HistoNbins, BinsVector);

    fNexTrainHist -> SetXTitle(x_axis_title);
    fNexTrainHist -> SetYTitle("Excess Events");


    fNexTestHist = new TH1F ("NexTestHist",
			     "NexTestHist",
			     HistoNbins, BinsVector);

    fNexTestHist -> SetXTitle(x_axis_title);
    fNexTestHist -> SetYTitle("Excess Events");



    fNEvtsInTrainMatrixONHist = new TH1F("NEvtsInTrainMatrixONHist",
					 "NEvtsInTrainMatrixONHist", 
					 HistoNbins, BinsVector);


    fNEvtsInTrainMatrixONHist -> SetXTitle(x_axis_title);
    fNEvtsInTrainMatrixONHist -> SetYTitle("Number of ON Events");
    

    fNEvtsInTestMatrixONHist = new TH1F("NEvtsInTestMatrixONHist",
					 "NEvtsInTestMatrixONHist", 
					 HistoNbins, BinsVector);

    fNEvtsInTestMatrixONHist -> SetXTitle(x_axis_title);
    fNEvtsInTestMatrixONHist -> SetYTitle("Number of ON Events");
    

    fNEvtsInTrainMatrixOFFHist = new TH1F("NEvtsInTrainMatrixOFFHist",
					  "NEvtsInTrainMatrixOFFHist", 
					  HistoNbins, BinsVector);

    fNEvtsInTrainMatrixOFFHist -> SetXTitle(x_axis_title);
    fNEvtsInTrainMatrixOFFHist -> SetYTitle("Number of OFF Events");


    fNEvtsInTestMatrixOFFHist = new TH1F("NEvtsInTestMatrixOFFHist",
					 "NEvtsInTestMatrixOFFHist", 
					 HistoNbins, BinsVector);

    fNEvtsInTestMatrixOFFHist -> SetXTitle(x_axis_title);
    fNEvtsInTestMatrixOFFHist -> SetYTitle("Number of OFF Events");
    

    // Histograms references are removed from the current directory

    /* IT DOES NOT WORK !!! THESE STATEMENTS PRODUCE SEGMENTATION FAULT.
    fNormFactorTrainHist -> SetDirectory(NULL);
    
    fNormFactorTestHist -> SetDirectory(NULL);
    
    fSigmaLiMaTrainHist -> SetDirectory(NULL);
    
    fSigmaLiMaTestHist -> SetDirectory(NULL);
    
    fSigmaLiMaTestHist -> SetDirectory(NULL);
    
    fNexTestHist -> SetDirectory(NULL);
    */


    return kTRUE;

}

void MFindSupercutsONOFFThetaLoop::WriteSuccessfulThetaBinsHistToFile()
{
    *fLog << "MFindSupercutsONOFFThetaLoop : Writing histogram "
	  << " 'SuccessfulThetaBinsHist'"
	  << " (if contents != 0) into root file " 
	  << fAlphaDistributionsRootFilename << endl;

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


    if (fSuccessfulThetaBinsHist)
    {
	if (fSuccessfulThetaBinsHist->GetEntries() > 0.5)
	    fSuccessfulThetaBinsHist -> Write();
    }


    rootfile.Close();

}

void MFindSupercutsONOFFThetaLoop::WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile()
{

    *fLog << "MFindSupercutsONOFFThetaLoop : Writing histograms 'NexSigmaLiMaNormFactorNEvtsTrainTest'"
	  << " (if they have contents != 0) into root file " 
	  << fAlphaDistributionsRootFilename << endl;

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


    if (fNormFactorTrainHist)
    {
	if (fNormFactorTrainHist->GetEntries() > 0.5)
	    fNormFactorTrainHist -> Write();
    }

    if (fNormFactorTestHist)
    {
	if (fNormFactorTestHist->GetEntries() > 0.5)
	    fNormFactorTestHist -> Write();
    }

    if (fSigmaLiMaTrainHist)
    {
	if (fSigmaLiMaTrainHist->GetEntries() > 0.5)
	    fSigmaLiMaTrainHist -> Write();
    }

    if (fSigmaLiMaTestHist)
    {
	if (fSigmaLiMaTestHist->GetEntries() > 0.5)
	    fSigmaLiMaTestHist -> Write();
    }

    if (fNexTrainHist)
    {
	if (fNexTrainHist->GetEntries() > 0.5)
	    fNexTrainHist -> Write();
    }

    if (fNexTestHist)
    {
	if (fNexTestHist->GetEntries() > 0.5)
	    fNexTestHist -> Write();
    }
   

    if (fNEvtsInTrainMatrixONHist)
    {
	if (fNEvtsInTrainMatrixONHist->GetEntries() > 0.5)
	    fNEvtsInTrainMatrixONHist -> Write();
    }

    if (fNEvtsInTestMatrixONHist)
    {
	if (fNEvtsInTestMatrixONHist->GetEntries() > 0.5)
	    fNEvtsInTestMatrixONHist -> Write();
    }
    
    if (fNEvtsInTrainMatrixOFFHist)
    {
	if (fNEvtsInTrainMatrixOFFHist->GetEntries() > 0.5)
	    fNEvtsInTrainMatrixOFFHist -> Write();
    }


    if (fNEvtsInTestMatrixOFFHist)
    {
	if (fNEvtsInTestMatrixOFFHist->GetEntries() > 0.5)
	    fNEvtsInTestMatrixOFFHist -> Write();
    }


    rootfile.Close();

}

void  MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents(
    Double_t fontrain, 
    Double_t fontest, 
    Double_t fofftrain,
    Double_t fofftest)
{

    fWhichFractionTrain = fontrain;
    fWhichFractionTest = fontest;
    fWhichFractionTrainOFF = fofftrain;
    fWhichFractionTestOFF = fofftest;


}

void MFindSupercutsONOFFThetaLoop::SetReadMatricesFromFile (Bool_t b)
{

    *fLog << "---------------------------------------------------" << endl
	  << "MFindSupercutsONOFFThetaLoop: " << endl;

    if (b)
    {	*fLog  << "Matrices are read from root files." << endl;}
    else
    {   *fLog  << "Matrices are produced and stored into root files." 
	       << endl;
    }


    fReadMatricesFromFile = b;
}

Bool_t MFindSupercutsONOFFThetaLoop::SetAlphaDistributionsRootFilename(
    const TString &name)
{

    if (fPathForFiles.IsNull())
    {
	*fLog << err 
	      << "MFindSupercutsONOFFThetaLoop::SetAlphaDistributionsRootFilename; " 
	      << endl
	      << "fPathForFiles is not defined yet, hence "
	      << "name for rootfile containing alpha distributions " 
	      << "can not be defined properly" << endl;

	return kFALSE;
    }

    
    fAlphaDistributionsRootFilename = (fPathForFiles);
    fAlphaDistributionsRootFilename += (name);


    *fLog << "MFindSupercutsONOFFThetaLoop: fAlphaDistributionsRootFilename = "
	  << fAlphaDistributionsRootFilename << endl;


    return kTRUE;


}

// Function to set Names of root files containing matrices and  
// optimizedparameters and also fThetaRangeStringVector vector

Bool_t MFindSupercutsONOFFThetaLoop::SetALLNames()
{
    // Names for files only possible if CosThetaRangeVector 
    // is defined properly
    if (fCosThetaRangeVector.GetSize() < 2)
    {
        *fLog << err 
	      << "MFindSupercutsONOFFThetaLoop::SetALLNames; " 
	      << endl
	      << "Size of fCosThetaRangeVector is smaller than 2... "
	      << "fCosThetaRangeVector is not defined properly, "
	      << "and thus, names for rootfiles containing data matrices "
	      << "and optimized parameters for the different theta "
	      << "ranges can not defined properly." << endl;
 
        return kFALSE;
    }

    if (fPathForFiles.IsNull())
    {
	*fLog << err 
	      << "MFindSupercutsONOFFThetaLoop::Setnames; " 
	      << endl
	      << "fPathForFiles is not defined yet, hence "
	      << "names for rootfiles containing data matrices "
	      << "and optimized parameters for the different theta "
	      << "ranges can not defined properly." << endl;

	return kFALSE;
    }


     // Name of rootfiles for Supercuts parameters follow the 
    // follwing structure

    // Path + "OptSCParametersONOFF" + "ThetaRange" 
    // + int {(fThetaMin_fThetaMax)*1000} + ".root"


    // Name of rootfiles for Matrices follow the follwing structure

    // Path + "ON(OFF)Data" + "Train(Test)" + "Matrix" + Fraction 
    // + "ThetaRange" + int {(fThetaMin_fThetaMax)*1000} + ".root"


    Int_t tmpdim = fCosThetaRangeVector.GetSize() - 1;
    const Int_t VectorDimension = tmpdim;

    fThetaRangeStringVector = new TString[VectorDimension];

    fOptSCParamFilenameVector = new TString[VectorDimension];

    fTrainMatrixONFilenameVector = new TString[VectorDimension];
    fTestMatrixONFilenameVector = new TString[VectorDimension];
    fTrainMatrixOFFFilenameVector = new TString[VectorDimension];
    fTestMatrixOFFFilenameVector = new TString[VectorDimension];


    Int_t ThetaMinTimes1000 = 0;
    Int_t ThetaMaxTimes1000 = 0;
    Int_t FractionONTrainTimes100 = 0;
    Int_t FractionONTestTimes100 = 0;
    Int_t FractionOFFTrainTimes100 = 0;
    Int_t FractionOFFTestTimes100 = 0;

    char ThetaRangeString[100];

    for (Int_t i = 0; i < VectorDimension; i++)
    {
	if (!SetThetaRange(i+1))
	{
	    *fLog << "MFindSupercutsONOFFThetaLoop::SetALLNames; " 
		  << endl
		  << "Values for fThetaMin and fThetaMax could NOT "
		  << "be computed for theta bin " 
		  << (i+1) << "; SetNames can not go on successfully."
		  << endl;

	    return kFALSE;
	}
	
	
	ThetaMinTimes1000 = int(fThetaMin*1000);
	ThetaMaxTimes1000 = int(fThetaMax*1000);
	sprintf(ThetaRangeString, "ThetaRange%d_%dmRad", 
		ThetaMinTimes1000, ThetaMaxTimes1000);

	fThetaRangeStringVector[i] = (ThetaRangeString);
	
	// tmp
	*fLog << "ThetaRangeString = " << fThetaRangeStringVector[i] << endl;
	// endtmp

	Double_t tmpdouble = 0.0;
	tmpdouble = fWhichFractionTrain*100 - int(fWhichFractionTrain*100); 
	
	FractionONTrainTimes100 = tmpdouble < 0.5 ? int(fWhichFractionTrain*100):
	    int(fWhichFractionTrain*100) + 1;
	

	tmpdouble = fWhichFractionTest*100 - int(fWhichFractionTest*100);

	FractionONTestTimes100 = tmpdouble < 0.5 ? int(fWhichFractionTest*100):
	    int(fWhichFractionTest*100) + 1;
	
	tmpdouble = fWhichFractionTrainOFF*100 - int(fWhichFractionTrainOFF*100);

	FractionOFFTrainTimes100 = tmpdouble < 0.5 ? int(fWhichFractionTrainOFF*100):
	    int(fWhichFractionTrainOFF*100) + 1;

	tmpdouble = fWhichFractionTestOFF*100 - int(fWhichFractionTestOFF*100);
	FractionOFFTestTimes100 = tmpdouble < 0.5 ? int(fWhichFractionTestOFF*100):
	    int(fWhichFractionTestOFF*100) + 1;
	
	

	// File for SC parameters 
	fOptSCParamFilenameVector[i] = (fPathForFiles);
	fOptSCParamFilenameVector[i] += ("OptSCParametersONOFF");
	fOptSCParamFilenameVector[i] += (fThetaRangeStringVector[i]);
	fOptSCParamFilenameVector[i] += (".root");
	
	
	
	// File for Data Matrices

	fTrainMatrixONFilenameVector[i] = (fPathForFiles);
	fTrainMatrixONFilenameVector[i] += ("ONDataTrainMatrixFraction");
	fTrainMatrixONFilenameVector[i] += (FractionONTrainTimes100);
	fTrainMatrixONFilenameVector[i] += (fThetaRangeStringVector[i]);
	fTrainMatrixONFilenameVector[i] += (".root");

	fTestMatrixONFilenameVector[i] = (fPathForFiles);
	fTestMatrixONFilenameVector[i] += ("ONDataTestMatrixFraction");
	fTestMatrixONFilenameVector[i] += (FractionONTestTimes100);
	fTestMatrixONFilenameVector[i] += (fThetaRangeStringVector[i]);
	fTestMatrixONFilenameVector[i] += (".root");


	fTrainMatrixOFFFilenameVector[i] = (fPathForFiles);
	fTrainMatrixOFFFilenameVector[i] += ("OFFDataTrainMatrixFraction");
	fTrainMatrixOFFFilenameVector[i] += (FractionOFFTrainTimes100);
	fTrainMatrixOFFFilenameVector[i] += (fThetaRangeStringVector[i]);
	fTrainMatrixOFFFilenameVector[i] += (".root");

	fTestMatrixOFFFilenameVector[i] = (fPathForFiles);
	fTestMatrixOFFFilenameVector[i] += ("OFFDataTestMatrixFraction");
	fTestMatrixOFFFilenameVector[i] += (FractionOFFTestTimes100);
	fTestMatrixOFFFilenameVector[i] += (fThetaRangeStringVector[i]);
	fTestMatrixOFFFilenameVector[i] += (".root");

    }
    // Info concerning names is printed
    
    *fLog << "--------------------------------------------------" << endl
	  << "MFindSupercutsONOFFThetaLoop::Setnames; " 
	  << endl
	  << "Names for root files containing Opt SC Param. and Data "
	  << "matrices for the different theta bins are the "
	  << "following ones: " << endl
	  << "MeanCosTheta, OptSCParm, ONDataTrainMatrix, ..." << endl;
    
    for (Int_t i = 0; i < VectorDimension; i++)
    {
	*fLog << fCosThetaBinCenterVector[i] << ",  " << endl
	      << fOptSCParamFilenameVector[i] << ",  "<< endl
	      << fTrainMatrixONFilenameVector[i] << ",  " << endl
	      << fTestMatrixONFilenameVector[i] << ",  " << endl
	      << fTrainMatrixOFFFilenameVector[i] << ",  " << endl
	      << fTestMatrixOFFFilenameVector[i] << ",  " << endl
	      << endl << endl;
    }
    
    *fLog << "--------------------------------------------------" << endl;
    
    return kTRUE;
    
}



Bool_t MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges()
{


    *fLog << "--------------------------------------------------------------" << endl
	  << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges function starts " 
	  << endl;
	  


    if (fDataONRootFilename.IsNull())
    {
      *fLog << "MFindSupercutsONOFF::LoopOverThetaRanges; root file for ON data is not defined... " 
	    << "aborting..."
            << endl;
      return kFALSE;
    }


    if (fDataOFFRootFilename.IsNull())
    {
      *fLog << "MFindSupercutsONOFF:: LoopOverThetaRanges; root file for OFF data is not defined... aborting"
            << endl;
      return kFALSE;
    }

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


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


    if (fCosThetaRangeVector.GetSize() < 2)
	
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
	      << "Dimension of vector fCosThetaRangeVector " 
	      << "is smaller than 2" 
	      << "Theta ranges are not well defined... aborting ..." << endl;
	
	return kFALSE;
    }


    // Check if pointers to root file names (opt SC parameters and matrices) 
    // are empty


    if (!fTrainMatrixONFilenameVector || 
	!fTestMatrixONFilenameVector ||
	!fTrainMatrixOFFFilenameVector ||
	!fTestMatrixOFFFilenameVector ||
	!fOptSCParamFilenameVector)
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
	      << "Pointers to root files names for Opt SC parameters and data "
	      << "matrices are not well defined. Aborting... " << endl;

	return kFALSE;
    }


    if (fAlphaDistributionsRootFilename.IsNull())
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
	      << "Root filename where to store alpha distributions "
	      << "(fAlphaDistributionsRootFilename) is not defined." << endl
	      << "Loop over theta angles can not go on..." << endl;

	return kTRUE;
    }


    
	
    if (!fPsFilename)
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
	      << "Object from Class TPostScript is not defined. "
	      << "Plots will be stored in a a file named 'PsFilesTmp.ps' located "
	      << "in current directory" << endl;

	// fPsFilename = new TPostScript("PsFilesTmp.ps" ,111);


    }


    



    // All drinks are in house... the party starts !!

    Int_t ThetaLoopIterations = fCosThetaRangeVector.GetSize() - 1;
                                
    *fLog << "looping over the following Cos Theta Angle Ranges : "<< endl;
    for (Int_t i = 0; i < ThetaLoopIterations; i++)
    {
	*fLog << i+1 << " -> " 
	      << fCosThetaRangeVector[i] << "-" << fCosThetaRangeVector[i+1] 
	      << endl;
    }
	  
    *fLog << "--------------------------------------------------------------------" 
	  << endl;



    // Let's create histograms where SigmaLiMa, Nex and Normfactors will be stored

    if (!SetNexSigmaLiMaNormFactorNEvtsTrainTestHist())
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
	      << "Histograms where to store Significances, Nex and Normalization "
	      << "factors, NEvts for all theta bins could not be created... Aborting."
	      << endl;

	return kFALSE;
    }


    // Let's create histogram where successfully optimized theta bins will be stored

    if (!SetSuccessfulThetaBinsHist())
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
	      << "Histogram where to store successfully optimized theta bins "
	      << " could not be created... Aborting."
	      << endl;

	return kFALSE;
    }


    


    Double_t NormFactor = 0.0;
    Double_t Nex = 0.0;
    Double_t SigmaLiMa = 0.0;
    Int_t NEvtsON = 0;
    Int_t NEvtsOFF = 0;
    
    Bool_t MinimizationSuccessful = kTRUE;
    Double_t MinimizationValue = 1.0; // Value quantifying the minimization, and stored in 
    // histogram fSuccessfulThetaBinsHist. For the time being it is 0.0 (minimization failed). 
    // and 1.0 (minimization succeeded). Yet in future it may take other values quantifying 
    // other aspects...


    for (Int_t i = 0; i < ThetaLoopIterations; i++)
    {

	// Get the Theta range that will be used in this iteration
	// Remember that argument of the function is iteration number, 
	// and it starts at 1...

	if (!SetThetaRange(i+1))
	{
	    *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " 
		  << endl
		  << "Values for fThetaMin and fThetaMax could NOT "
		  << "be computed for theta bin " 
		  << (i+1) << ";  LoopOverThetaRanges is aborting."
		  << endl;

	    return kFALSE;
	}
	


	// Object of class MFindSupercutsONOFF will be created and used 


	MFindSupercutsONOFF FindSupercuts("MFindSupercutsONOFF", 
					     "Optimizer for the supercuts");

	
	// Set Theta range that will be used
	FindSupercuts.SetThetaRange(fThetaMin, fThetaMax);


	// Alphamax where signal is expected is set
	FindSupercuts.SetAlphaSig(fAlphaSig);

	// Bkg alpha region is set 
	FindSupercuts.SetAlphaBkgMin(fAlphaBkgMin);
	FindSupercuts.SetAlphaBkgMax(fAlphaBkgMax);

	// alpha bkg and signal region set in object FindSupercuts
	// are re-checked in order to be sure that make sense
	
	FindSupercuts.CheckAlphaSigBkg();


	// bining for alpha plots is set 
	
	FindSupercuts.SetAlphaPlotBinining(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);



	
	// Bool variable NormFactorFromAlphaBkg is set in FindSupercuts
	// in order to decide the way in which mormalization factors
	// for TRAIN and TEST samples are computed


	FindSupercuts.SetNormFactorFromAlphaBkg(fNormFactorFromAlphaBkg);


	// Define size range of events used to fill the data matrices

	FindSupercuts.SetSizeRange(fSizeCutLow, fSizeCutUp);



	
	// 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.

	FindSupercuts.SetSkipOptimization(fSkipOptimization);

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

	FindSupercuts.SetUseInitialSCParams(fUseInitialSCParams);


        // Set boolean variable that controls wether the theta variable 
	// is used or not in the computation of the dynamica cuts

	FindSupercuts.SetVariableNotUseTheta(fNotUseTheta);

	// Set boolean variable that controls wether cuts are dynamical 
	// or static.

	FindSupercuts.SetVariableUseStaticCuts(fUseStaticCuts);

	
	FindSupercuts.SetUseFittedQuantities(fUseFittedQuantities);

	                                        
	// fGammaEfficiency is set in object FindSupercuts

	FindSupercuts.SetGammaEfficiency (fGammaEfficiency);


	// Filename with SuperCuts parameters is set

	FindSupercuts.SetFilenameParam(fOptSCParamFilenameVector[i]);

	// Hadroness containers are set
	
	FindSupercuts.SetHadronnessName(fHadronnessName);
	FindSupercuts.SetHadronnessNameOFF(fHadronnessNameOFF);
    
	// Rootfilename where to store alpha distribution histograms is set

	FindSupercuts.SetAlphaDistributionsRootFilename(fAlphaDistributionsRootFilename);


	// Set the names of all trees 
        // used to store the supercuts applied to ON/OFF Train/Test samples

	// In future, names might be controlled by user, but for the time 
	// being I'll make things simple, and names are set to some defaults 
	// taking into account the theta range used
	
	// Containers will be stored in rootfile specified by variable 
	// fAlphaDistributionsRootFilename

	FindSupercuts.SetSupercutsAppliedTreeNames();





	// Object of the class TPostScritp is set...
	// BUT THAT'S NOT WORKING PROPERLY !!!!!

	// FindSupercuts.SetPostScriptFile (fPsFilename);

	// ********************************************************
	// Due to the failure of the use of object TPostScript
	// to make a Ps document with all plots, I decided to use the 
	// standard way (SaveAs(filename.ps)) to store plots related 
	// to alpha distributions 
	// for ON and OFF and BEFORE and AFTER cuts (VERY IMPORTANT
	// TO COMPUTE EFFICIENCIES IN CUTS) 

	// Psfilename is set inside function LoopOverThetaRanges()
	// and given to the object MFindSupercutsONOFF created 
	// within this loop.

	// This will have to be removed as soon as the TPostScript
	// solutions works...
	// ********************************************************

	TString PsFilenameString = (fPathForFiles);
	PsFilenameString += ("SC");
	PsFilenameString += (fThetaRangeStringVector[i]);
	// Name to be completed inside object MFindSupercutsONOFF

	FindSupercuts.SetPsFilenameString(PsFilenameString);


	if (fReadMatricesFromFile)
	{
	    *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
		  << "Reading Data matrices from root files... " << endl;

	    FindSupercuts.ReadMatrix(fTrainMatrixONFilenameVector[i], 
				     fTestMatrixONFilenameVector[i]);
	    FindSupercuts.ReadMatrixOFF(fTrainMatrixOFFFilenameVector[i], 
					fTestMatrixOFFFilenameVector[i]);

	}

	else
	{
	    *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
		  << "Reading Data from files " << fDataONRootFilename 
		  << " and " << fDataOFFRootFilename << "(for ON and OFF data "
		  << "respectively " 
		  << " and producing data matrices that will be stored "
		  << "in root files. " << endl;

	    // ON data 
	    
	    FindSupercuts.DefineTrainTestMatrixThetaRange(
		fDataONRootFilename,
		fWhichFractionTrain,
		fWhichFractionTest, 
		fThetaMin, fThetaMax, 
		fTrainMatrixONFilenameVector[i], 
		fTestMatrixONFilenameVector[i]);

	    // OFF data
	    
	    FindSupercuts.DefineTrainTestMatrixOFFThetaRange(
		fDataOFFRootFilename,
		fWhichFractionTrainOFF,
		fWhichFractionTestOFF, 
		fThetaMin, fThetaMax, 
		fTrainMatrixOFFFilenameVector[i], 
		fTestMatrixOFFFilenameVector[i]);
	    
	}


	MinimizationSuccessful = kTRUE;
	
	if (fOptimizeParameters)
	{
	    MinimizationSuccessful = 
		FindSupercuts.FindParams("",fInitSCPar, fInitSCParSteps); 
	    
	    if (MinimizationSuccessful)
	    {
		// Minimization was successful	
		// NormFactor, Nex, SigmaLiMa and NEvts Train histograms are updated 

		*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " 
		      << endl
		      << "Minimization was successful" << endl
		      << "Updating values for NormFactor, Nex, SigmaLiMa and NEvts " 
		      << "Train histograms for " <<  fThetaRangeStringVector[i] << endl;
		
		
		NormFactor = FindSupercuts.GetNormFactorTrain();
		SigmaLiMa = FindSupercuts.GetSigmaLiMaTrain();
		Nex = FindSupercuts.GetNexTrain();
		
		NEvtsON = FindSupercuts.GetMatrixTrain()->GetM().GetNrows();
		NEvtsOFF = FindSupercuts.GetMatrixTrainOFF()->GetM().GetNrows();

		// MinimizationValue is updated...

		MinimizationValue = 1.0;		

	    }
	    else
	    {
		// Minimization was NOT successful	
		// NormFactor, Nex, SigmaLiMa and NEvts Train histograms are updated 

		*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " 
		      << endl
		      << "Minimization was NOT successful" << endl
		      << "Updating values for NormFactor, and NEvts " 
		      << "Train histograms for " <<  fThetaRangeStringVector[i] << endl
		      << "Nex and SigmaLiMa are set to ZERO" << endl;
		
		
		NormFactor = FindSupercuts.GetNormFactorTrain();
		
		NEvtsON = FindSupercuts.GetMatrixTrain()->GetM().GetNrows();
		NEvtsOFF = FindSupercuts.GetMatrixTrainOFF()->GetM().GetNrows();

		SigmaLiMa = 0.0;
		Nex = 0.0;
		
		// MinimizationValue is updated...

		MinimizationValue = 0.0;

	    }


	    
	    cout << "---- Train Histograms will be updated with the following numbers: ------ " 
		 << endl
		 << "Mean CosTheta for this Thetabin = " 
		 << fCosThetaBinCenterVector[i] << endl
		 << "Minimization status (1.0 is successful, 0.0 is UNsuccessful) = "
		 << MinimizationValue << endl
		 << "Number of ON Events in Train sample used = "
		 << NEvtsON << endl
		 << "Number of OFF Events in Train sample used = "
		 << NEvtsOFF << endl
		 << "Normalization factor (Non/Noff) for Train sample = "
		 << NormFactor << endl 
		 << "Excess events and SigmaLiMa: "
		 << Nex << "; " << SigmaLiMa << endl << endl;
		 
	    

	    fNormFactorTrainHist -> Fill(fCosThetaBinCenterVector[i], 
					  NormFactor);
	    
	    fNexTrainHist -> Fill(fCosThetaBinCenterVector[i], 
				   Nex);

	    fSigmaLiMaTrainHist -> Fill(fCosThetaBinCenterVector[i], 
					 SigmaLiMa);


	    fNEvtsInTrainMatrixONHist -> Fill(fCosThetaBinCenterVector[i], 
					 NEvtsON);

	    fNEvtsInTrainMatrixOFFHist -> Fill(fCosThetaBinCenterVector[i], 
						NEvtsOFF);

	    
	    fSuccessfulThetaBinsHist -> Fill(fCosThetaBinCenterVector[i],
					     MinimizationValue);
	    
	    cout << "Histograms (TRAIN) updated successfully " << endl;
	    

	    
	    
	}
	
	
	if (fTestParameters && MinimizationSuccessful)
	{
	    
	    // Cuts are applied to test sample
	    
	    FindSupercuts.TestParamsOnTestSample();

	    // NormFactor, Nex and SigmaLiMa Test histograms are updated 

	    *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
		  << "Updating values for NormFactor, Nex, SigmaLiMa and NEvts " 		  
		  << "Test histograms for " <<  fThetaRangeStringVector[i] << endl;


	    NormFactor = FindSupercuts.GetNormFactorTest();
	    SigmaLiMa = FindSupercuts.GetSigmaLiMaTest();
	    Nex = FindSupercuts.GetNexTest();

	    NEvtsON = FindSupercuts.GetMatrixTest()->GetM().GetNrows();
	    NEvtsOFF = FindSupercuts.GetMatrixTestOFF()->GetM().GetNrows();


       cout << "---- Test Histograms will be updated with the following numbers: ------ " 
	    << endl
		 << "Mean CosTheta for this Thetabin = " 
		 << fCosThetaBinCenterVector[i] << endl
		 << "Number of ON Events in Test sample used = "
		 << NEvtsON << endl
		 << "Number of OFF Events in Test sample used = "
		 << NEvtsOFF << endl
		 << "Normalization factor (Non/Noff) for Test sample = "
		 << NormFactor << endl 
		 << "Excess events and SigmaLiMa: "
		 << Nex << "; " << SigmaLiMa << endl << endl;
		 



	    fNormFactorTestHist -> Fill(fCosThetaBinCenterVector[i], 
					 NormFactor);
	    
	    fNexTestHist -> Fill(fCosThetaBinCenterVector[i], 
				  Nex);

	    fSigmaLiMaTestHist -> Fill(fCosThetaBinCenterVector[i], 
				       SigmaLiMa);


	    fNEvtsInTestMatrixONHist -> Fill(fCosThetaBinCenterVector[i], 
					    NEvtsON);

	    fNEvtsInTestMatrixOFFHist -> Fill(fCosThetaBinCenterVector[i], 
					      NEvtsOFF);

	    
	    cout << "Histograms (TEST) updated successfully " << endl;
	    



	}

    }


    // PsFile is closed


    //*fLog << "Closing PsFile ";
    //*fLog << fPsFilename -> GetName() << endl;
    
    //fPsFilename -> Close();

    

    if (fOptimizeParameters || fTestParameters)
    {
	// Histograms cotanining SigmaLiMa, Nex and NormFactors and NEvts are written 
	// in the same root file where alpha distributions where stored
	// i.e. fAlphaDistributionsRootFilename
	
	*fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
	      << "Executing WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile()"
	      << endl;

	WriteNexSigmaLiMaNormFactorNEvtsTrainTestHistToFile();


	if (fOptimizeParameters)
	{   // Histogram containing minimization status for all theta bins defined 
	    // in vector fCosThetaVector is stored in file fAlphaDistributionsRootFilename

	    *fLog << "MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges; " << endl
	      << "Executing WriteSuccessfulThetaBinsHistToFile()"
	      << endl;

	    WriteSuccessfulThetaBinsHistToFile();
	    
	}

    }


    return kTRUE;
}





 // Function that loops over the alpha distributions (ON-OFF)
   // stored in root file defined by fAlphaDistributionsRootFilename
   // and computes the significance and Nex (using MHFindSignificanceONOFF::FindSigma)
   // for several cuts in alpha (0-fAlphaSig; in bins defined for alpha distributions 
   // by user).

   // It initializes the histograms (memeber data), fills them and store them 
   // in root file defined by fAlphaDistributionsRootFilename. A single histogram 
   // for each theta bin.



Bool_t MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()
{


  *fLog << endl << endl 
	<< "****************************************************************************" << endl
	<< "Computing number of excess events and significance vs " << endl
	<< " the cut in alpha parameter (fAlphaSig) using function " << endl
	<< " MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig() " << endl
	<< "****************************************************************************" << endl
	<< endl;
  

  // Check that all elements needed by the function are available
 

  if (fAlphaDistributionsRootFilename.IsNull())
    {
	*fLog << err 
	      << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig; " 
	      << endl
	      << "fAlphaDistributionsRootFilename is not defined yet, hence "
	      << "histograms  containing alpha distributions for the different theta "
	      << "ranges can not be retrieved. " << endl;

	return kFALSE;
    }



    if (fCosThetaRangeVector.GetSize() < 2)
    {
        *fLog << err 
	      << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig; " 
	      << endl
	      << "Size of fCosThetaRangeVector is smaller than 2... "
	      << "fCosThetaRangeVector is not defined properly, "
	      << "and thus, number of theta loops and names "
	      << " of the histograms  containing alpha distributions "
	      << "for the different theta "
	      << "ranges can not defined properly." << endl;
 
        return kFALSE;
    }

   

    // It seems that i have all I need for the time being. Let's start the party...


    Int_t tmpint =  0; // variable that will be used to store temporally integers 
    Double_t tmpdouble = 0.0; // variable that will be used to store temporally doubles

    tmpint = fCosThetaRangeVector.GetSize() - 1;
    const Int_t NThetaBins = tmpint;
 

    // Definition of the bins for Nex and Significance vs Alpha histograms
    // The bins in Alpha will start with ZERO  and will end 
    // with 5 bins further after fAlphaSig.

    Double_t AlphaBinWidth = double((fAlphaBinUp-fAlphaBinLow)/fNAlphaBins);
    Double_t FirstBinLow = 0.0;
    tmpint = int((fAlphaSig/AlphaBinWidth) + 5.0);

    while ((FirstBinLow + AlphaBinWidth*tmpint) > fAlphaBkgMin)
      {
	tmpint--;
      }
	
    *fLog << "Information about the binning of the histograms Nex and Significance vs alphasig " 
	  << endl
	  << "FirstBinLow, AlphaBinWidth, NBins : " 
	  << FirstBinLow << ", " << AlphaBinWidth << ", " << tmpint-1 << endl;



    const Int_t NexSignVSAlphaXBinsVectorDimension = tmpint;
                
 
    Double_t* XBinsVector = new Double_t [NexSignVSAlphaXBinsVectorDimension];
    
    for (Int_t i = 0; i < NexSignVSAlphaXBinsVectorDimension; i++)
      {
	XBinsVector[i] = FirstBinLow + i*AlphaBinWidth;
	// cout << XBinsVector[i] << endl;
      }


   
    // Axis labels for Nex and Significance vs Alpha histograms

    TString x_axis_label = ("Cut in alpha [\\circ]");
    TString y_axis_label_nex = ("Excess Events");
    TString y_axis_label_significance = ("Significance (LiMa formula 17)");


    // Arguments for function MHFindSignificanceONOFF::FindSigmaONOFF()
    // are defined in here. The same arguments (obviously) will be used 
    // in all theta bins and TRAIN and TEST samples.

    // variable alphasig will be defined and varied withing the theta loop

    Double_t alphasig = 0.0;

    // Minimum value of alpha for bkg region in ON data
    const Double_t alphabkgmin = fAlphaBkgMin;
    
    // Maximum value of alpha for bkg region in ON data
    const Double_t alphabkgmax = fAlphaBkgMax;
    
    // degree of polynomial function used to fit ON data in background region
    const Int_t    degree   =    2;
    
    // degree of polynomial function used to fit OFF data in all alpha region (0-90)
    const Int_t    degreeOFF   =    2;
    
    const Bool_t   drawpoly  = kFALSE;
    const Bool_t   fitgauss  = kFALSE;
    const Bool_t   print     = kFALSE;
    
    const Bool_t saveplots   = kFALSE; // No plots are saved
   
    const Bool_t RebinHistograms = kTRUE;
    const Bool_t ReduceDegree = kFALSE;



    Double_t Nex = 0.0;
    Double_t Significance = 0.0;


    // Names of histograms containing Nex and Significance vs alphasig

    TString NexVSAlphaName = ("NexVSAlpha"); // Names to be completed in theta loop
    TString SignificanceVSAlphaName = ("SignificanceVSAlpha");// Names to be completed in theta loop
    


    // Names of histograms to be retrieved from root file 
    // fAlphaDistributionsRootFilename

    TString NormFactorTrainHistoName = ("NormFactorTrainHist");
    TString NormFactorTestHistoName = ("NormFactorTestHist");
    
    TString SuccessfulThetaBinsHistoName = ("SuccessfulThetaBinsHist");

 

    // Histogram containing information about the success of the optimization
    // for the several theta bins. The content of those bins where optimization was not 
    // successful, are set to ZERO.

    
    TH1F* SuccessfulThetaBinsHisto;
 


    // Histograms to store the normalization factors

    TH1F* NormFactorTrainHisto;
    TH1F* NormFactorTestHisto;



    // Histograms that will be used to store the alpha distributions 
    // retrieved from   fAlphaDistributionsRootFilename, and that will be use as arguments in 
    // function MHFindSignificanceONOFF::FindSigma

    // The same histograms (pointer to histograms actually) will be used for TRAIN and TEST 
    // and for the different ThetaNBinsUsed.

    TH1F* AlphaONAfterCuts;
    TH1F* AlphaOFFAfterCuts;



   
      // Vector containing the histogram index of the successfully 
    // optimized theta bins. 
    // It is filled with the contents of the histogram 
    // fSuccessfulThetaBinsHist retrieved from the root 
    // file fAlphaDistributionsRootFilename
    // Remember that histogram index START AT 1 !!!
    TArrayI SuccessfulThetaBinsVector; 
   


    



     // Open root file for retrieving information
 
    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl
	  << "Opening root file " << fAlphaDistributionsRootFilename << endl;

    TFile rootfile (fAlphaDistributionsRootFilename, "READ");

    // Retrieving SuccessfulThetaBinsHisto from root file. 
    
    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl
	  << "Retrieving SuccessfulThetaBinsHisto from root file and filling vector "
	  << " SuccessfulThetaBinsVector" << endl;

                                    
    cout << "Getting histogram with name " << SuccessfulThetaBinsHistoName << endl;

    SuccessfulThetaBinsHisto = (TH1F*) rootfile.Get(SuccessfulThetaBinsHistoName);
    


    // Check that bins in this histo correspond with the 
    // CosTheta intervals defined by vector fCosThetaRangeVector

	
    if (SuccessfulThetaBinsHisto -> GetNbinsX() != NThetaBins)
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" 
	      << endl
	      << "Number of theta bins defined by vector fCosThetaRangeVector (" 
	      << NThetaBins << ") does not match with the bins in histogram "
	      << SuccessfulThetaBinsHistoName << "(" 
	      << SuccessfulThetaBinsHisto -> GetNbinsX() << ")" << endl
	      << "Function execution will be ABORTED." << endl;
	return kFALSE;
    }


    // Filling vector SuccessfulThetaBinsVector

    tmpint = 0;
    for (Int_t i = 0; i < NThetaBins; i++)
    {
	if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5)
	{
	    tmpint++;
	}
    }

    SuccessfulThetaBinsVector.Set(tmpint);
    
    tmpint = 0;
    for (Int_t i = 0; i < NThetaBins; i++)
    {
	if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5)
	{
	    if(tmpint < SuccessfulThetaBinsVector.GetSize())
	    {
		SuccessfulThetaBinsVector[tmpint] = i+1;
		tmpint++;
	    }
	    else
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" 
		      << endl
		      << "Problem when filling vector 'SuccessfulThetaBinsVector'" 
		      << endl
		      << "Function execution will be ABORTED." << endl;
		return kFALSE;
	    }
	    

	}
    }

    // variable defining the theta bins where optimization was successful, and 
    // hence, alpha distributions stored in root file are meaningful (up to some extent)
   
    tmpint = SuccessfulThetaBinsVector.GetSize();
    const Int_t SuccessfulNThetaBins = tmpint;
    
    // TEMP
    //  cout << "const Int_t SuccessfulNThetaBins = " << SuccessfulNThetaBins << endl;
    // ENDTEMP

    // Vectors of pointers to histograms are defined to store the histograms 
    // containing Nex and significance for TRAIN and TEST sample for all the 
    // successful theta bins

    // Histograms will be initialized in the loop over theta bins
          
    TH1F* NexVSAlphaSigTrainVector[SuccessfulNThetaBins];
    TH1F* NexVSAlphaSigTestVector[SuccessfulNThetaBins];

    
    TH1F* SignificanceVSAlphaSigTrainVector[SuccessfulNThetaBins];
    TH1F* SignificanceVSAlphaSigTestVector[SuccessfulNThetaBins];
    
    
    // ***************************************************
    // HISTOGRAMS FOR TRAIN SAMPLE ARE COMPUTED
    // ***************************************************

    

    if (fOptimizeParameters)
      { // Computation for TRAIN sample


	*fLog << endl << endl
	      << "****************************************************************************" << endl
	      << "Computing Nex and significance vs cut in alpha (alphasig)  for TRAIN sample" << endl
	      << "****************************************************************************" << endl
	      << endl << endl;


	// Retrieving NormFactorTrainHisto from root file. 

	 cout << "Getting histogram with name " << NormFactorTrainHistoName << endl;

	NormFactorTrainHisto = (TH1F*) rootfile.Get(NormFactorTrainHistoName);
	

	// Check that bins in this histo correspond with the 
	// CosTheta intervals defined by vector fCosThetaRangeVector

	////////////////// START CHECK  //////////////////////////////////////
	
	if (NormFactorTrainHisto -> GetNbinsX() != NThetaBins)
	{
	    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" 
		  << endl
		  << "Number of theta bins defined by vector fCosThetaRangeVector (" 
		  << NThetaBins << ") does not match with the bins in histogram "
		  << NormFactorTrainHistoName << "(" 
		  << NormFactorTrainHisto -> GetNbinsX() << ")" << endl
		  << "Function execution will be ABORTED." << endl;
	    return kFALSE;
	}

	// histo test 
	/*
	cout << "NORMFACTORTRAINHIST  Test: BinCenter and Value" << endl;
	for (Int_t k = 0; k < NThetaBins; k++)
	{
	    
	    cout << NormFactorTrainHisto -> GetBinCenter(k+1) 
		 << ";  " << NormFactorTrainHisto -> GetBinContent(k+1)<< endl;
	}
	*/

	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	{
	    tmpdouble = NormFactorTrainHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]);
	    if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble || 
		fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble)
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" 
		      << endl
		      << "Bins defined by vector fCosThetaRangeVector " 
		      << "do not match with the bins of histogram "
		      << NormFactorTrainHistoName << endl
		      << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1]
		      << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl
		      << "NormFactorHist Bin: " <<  tmpdouble << endl
		      << "Function execution will be ABORTED." << endl;

		return kFALSE;
	    }
	    
	}

	////////////////// END CHECK  //////////////////////////////////////



	
	// Check that the theta range string (needed to compute names for alpha histograms)
	// is well defined
	
	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	  {
	    if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
	      {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" 
		      << endl
		      << "Component " << SuccessfulThetaBinsVector[i]-1 
		      << " of fThetaRangeStringVector is EMPTY"
		      << endl
		      << "Function execution will be ABORTED." << endl;

		return kFALSE;

		
	    }
	  }


	// *************************************************************
	// *************************************************************
	
	// Normfactors and successfultheta bins are ok. 
	// I can now compute the Nex and Sigma for train sample (ON-OFF) 
	
	// *************************************************************
	// *************************************************************


	

	// Loop over the several successful theta bins doing the following 

	// 0) Initializing histograms Nex and Significance (with names!!)
	// 1) Retrieving alpha histograms (ON-OFF), and normalization factor
	// 2) Computing Nex and Signficance vs alpha and filling histograms


	*fLog << endl
	      << "Starting loop over the successfully optimized theta ranges " << endl
	      << endl;

	Double_t NormalizationFactor = 0.0;
	
	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	  {

	    *fLog << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] << endl << endl;
	
	    // 0) Initializing histograms Nex and Significance (with names!!)

	    TString NexVSAlphaNameTmp = (NexVSAlphaName);
	    TString SignificanceVSAlphaNameTmp = (SignificanceVSAlphaName);

	    NexVSAlphaNameTmp += "Train";
	    NexVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    SignificanceVSAlphaNameTmp += "Train";
	    SignificanceVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	
	    // TEMP    
	    // cout << NexVSAlphaNameTmp << endl;
	    // cout << SignificanceVSAlphaNameTmp << endl;
	    // ENDTEMP
	      


    
	    NexVSAlphaSigTrainVector[i] = new TH1F (NexVSAlphaNameTmp.Data(), 
						    NexVSAlphaNameTmp.Data(), 
						    NexSignVSAlphaXBinsVectorDimension - 1,
						    XBinsVector);
	    
	    // It is important to unlink the the created histograms from 
	    // the current directory (defined now by "rootfile"), 
	    // so that we can do whatever we want to 
	    // with them regardless of such directory is opened or closed 
	    // The reference of the histogram is removed from the current 
	    // directory by means of the member function  "SetDirectory(0)"


	    NexVSAlphaSigTrainVector[i] -> SetDirectory(0);


	    SignificanceVSAlphaSigTrainVector[i] = new TH1F (SignificanceVSAlphaNameTmp.Data(), 
							     SignificanceVSAlphaNameTmp.Data(),
							     NexSignVSAlphaXBinsVectorDimension - 1,
							     XBinsVector);
	    
	    SignificanceVSAlphaSigTrainVector[i] -> SetDirectory(0);
	    
	    

	    // axis are defined

	    NexVSAlphaSigTrainVector[i] -> SetXTitle(x_axis_label.Data());
	    NexVSAlphaSigTrainVector[i] -> SetYTitle(y_axis_label_nex.Data());
	    
	    SignificanceVSAlphaSigTrainVector[i] ->SetXTitle(x_axis_label.Data());
	    SignificanceVSAlphaSigTrainVector[i] ->SetYTitle(y_axis_label_significance.Data());

	
	    // 1) Retrieving alpha histograms (ON-OFF) and normalization factor

	    NormalizationFactor = 
	      NormFactorTrainHisto -> GetBinContent(SuccessfulThetaBinsVector[i]);


	    TString TmpAlphaHistoName = ("TrainAlphaAfterCuts");
	    TmpAlphaHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    *fLog << "Getting Histogram with name " << endl
		  << TmpAlphaHistoName << endl;

	    AlphaONAfterCuts = (TH1F*) rootfile.Get(TmpAlphaHistoName);

	    
	    TString TmpAlphaOFFHistoName = ("TrainAlphaOFFAfterCuts");
	    TmpAlphaOFFHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    *fLog << "Getting Histogram with name " << endl
		  << TmpAlphaOFFHistoName << endl;

	    AlphaOFFAfterCuts = (TH1F*) rootfile.Get(TmpAlphaOFFHistoName);
	    

	 // 2) Computing Nex and Signficance vs alpha and filling histograms
	    // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF  


	    *fLog  << endl << endl
		   << "Starting loop to compute Nex and Significance vs alphasig " << endl
		   << "for TRAIN sample and " 
		   << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] 
		   << " :" << endl;
	    
	   	    
	    for (Int_t j = 1; j < NexSignVSAlphaXBinsVectorDimension; j++)
	      {
		

		alphasig = XBinsVector[j];
		
		MHFindSignificanceONOFF findsig;
		findsig.SetRebin(RebinHistograms);
		findsig.SetReduceDegree(ReduceDegree);	
		findsig.SetUseFittedQuantities(fUseFittedQuantities);

		// Dummy name for psfile 
		TString psfilename = ("DummyPs"); 
	    
		findsig.FindSigmaONOFF(AlphaONAfterCuts, AlphaOFFAfterCuts, 
				       NormalizationFactor, 
				       alphabkgmin, alphabkgmax, 
				       degree, degreeOFF,
				       alphasig, 
				       drawpoly, fitgauss, print, 
				       saveplots, psfilename);

		
		Significance = double(findsig.GetSignificance());

		// so far there is not a variable that contains 
		// NexONOFF or NexONOFFFitted (in class MHFindSignificanceONOFF) 
		// according to the value of the variable fUseFittedQuantities.
		// As soon as I have some time I will implement it (with the 
		// corresponging GET member function) and will change 
		// MFindSignificanceONOFF accordingly to make use of it.
		
		// For the time being, I will survive with an "if statement".

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


				
		*fLog << endl << "Cut in alpha, Nex and significance : " 
		      << alphasig << ", " << Nex << ", " << Significance << endl << endl;
		
	      
		// updating histograms

		tmpdouble = alphasig - AlphaBinWidth/2.0;

		NexVSAlphaSigTrainVector[i] -> Fill(tmpdouble, Nex);
		SignificanceVSAlphaSigTrainVector[i] -> Fill(tmpdouble, Significance);
		
	  

		



	      }	  
	

	    // Dynamic memory allocated in this loop is released
	    
    
	    *fLog  << "Memory allocated by temporal alpha histograms "
		   << "will be released" << endl;
	    
	    delete AlphaONAfterCuts;
	    delete AlphaOFFAfterCuts;
	    
	  }


      }

    
    
     // ***************************************************
    // HISTOGRAMS FOR TEST SAMPLE ARE COMPUTED
    // ***************************************************


 

    if(fTestParameters)
      {// computation for TEST sample
	
	*fLog << endl << endl
	      << "**************************************************************************" << endl
	      << "Computing Nex and significance vs cut in alpha (alphasig) for TEST sample" << endl
	      << "**************************************************************************" << endl
	      << endl << endl;

	// Retrieving NormFactorTestHisto from root file. 

	 *fLog << "Getting histogram with name " << NormFactorTestHistoName << endl;

	NormFactorTestHisto = (TH1F*) rootfile.Get(NormFactorTestHistoName);
	

	// Check that bins in this histo correspond with the 
	// CosTheta intervals defined by vector fCosThetaRangeVector

	////////////////// START CHECK  //////////////////////////////////////
	
	if (NormFactorTestHisto -> GetNbinsX() != NThetaBins)
	{
	    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" 
		  << endl
		  << "Number of theta bins defined by vector fCosThetaRangeVector (" 
		  << NThetaBins << ") does not match with the bins in histogram "
		  << NormFactorTestHistoName << "(" 
		  << NormFactorTestHisto -> GetNbinsX() << ")" << endl
		  << "Function execution will be ABORTED." << endl;
	    return kFALSE;
	}

	// histo test 
	/*
	cout << "NORMFACTORTRAINHIST  Test: BinCenter and Value" << endl;
	for (Int_t k = 0; k < NThetaBins; k++)
	{
	    
	    cout << NormFactorTestHisto -> GetBinCenter(k+1) 
		 << ";  " << NormFactorTestHisto -> GetBinContent(k+1)<< endl;
	}
	*/

	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	{
	    tmpdouble = NormFactorTestHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]);
	    if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble || 
		fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble)
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" 
		      << endl
		      << "Bins defined by vector fCosThetaRangeVector " 
		      << "do not match with the bins of histogram "
		      << NormFactorTestHistoName << endl
		      << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1]
		      << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl
		      << "NormFactorHist Bin: " <<  tmpdouble << endl
		      << "Function execution will be ABORTED." << endl;

		return kFALSE;
	    }
	    
	}

	////////////////// END CHECK  //////////////////////////////////////



	
	// Check that the theta range string (needed to compute names for alpha histograms)
	// is well defined
	
	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	  {
	    if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
	      {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" 
		      << endl
		      << "Component " << SuccessfulThetaBinsVector[i]-1 
		      << " of fThetaRangeStringVector is EMPTY"
		      << endl
		      << "Function execution will be ABORTED." << endl;

		return kFALSE;

		
	    }
	  }


	// *************************************************************
	// *************************************************************
	
	// Normfactors and successfultheta bins are ok. 
	// I can now compute the Nex and Sigma for test sample (ON-OFF) 
	
	// *************************************************************
	// *************************************************************

	// Loop over the several successful theta bins doing the following 

	// 0) Initializing histograms Nex and Significance (with names!!)
	// 1) Retrieving alpha histograms (ON-OFF), and normalization factor
	// 2) Computing Nex and Signficance vs alpha and filling histograms


	*fLog << endl
	      << "Starting loop over the successfully optimized theta ranges " << endl
	      << endl;

	
	Double_t NormalizationFactor = 0.0;
	
	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	  {

	    *fLog << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] << endl << endl;

	    // 0) Initializing histograms Nex and Significance (with names!!)

	    TString NexVSAlphaNameTmp = (NexVSAlphaName);
	    TString SignificanceVSAlphaNameTmp = (SignificanceVSAlphaName);

	    NexVSAlphaNameTmp += "Test";
	    NexVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    SignificanceVSAlphaNameTmp += "Test";
	    SignificanceVSAlphaNameTmp += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	
	    // TEMP    
	    // cout << NexVSAlphaNameTmp << endl;
	    // cout << SignificanceVSAlphaNameTmp << endl;
	    // ENDTEMP
	      


    
	    NexVSAlphaSigTestVector[i] = new TH1F (NexVSAlphaNameTmp.Data(), 
						    NexVSAlphaNameTmp.Data(), 
						    NexSignVSAlphaXBinsVectorDimension - 1,
						    XBinsVector);
	    
	    // It is important to unlink the the created histograms from 
	    // the current directory (defined now by "rootfile"), 
	    // so that we can do whatever we want to 
	    // with them regardless of such directory is opened or closed 
	    // The reference of the histogram is removed from the current 
	    // directory by means of the member function  "SetDirectory(0)"


	    NexVSAlphaSigTestVector[i] -> SetDirectory(0);


	    SignificanceVSAlphaSigTestVector[i] = new TH1F (SignificanceVSAlphaNameTmp.Data(), 
							     SignificanceVSAlphaNameTmp.Data(),
							     NexSignVSAlphaXBinsVectorDimension - 1,
							     XBinsVector);
	    
	    SignificanceVSAlphaSigTestVector[i] -> SetDirectory(0);
	    
	    

	    // axis are defined

	    NexVSAlphaSigTestVector[i] -> SetXTitle(x_axis_label.Data());
	    NexVSAlphaSigTestVector[i] -> SetYTitle(y_axis_label_nex.Data());
	    
	    SignificanceVSAlphaSigTestVector[i] ->SetXTitle(x_axis_label.Data());
	    SignificanceVSAlphaSigTestVector[i] ->SetYTitle(y_axis_label_significance.Data());

	
	    // 1) Retrieving alpha histograms (ON-OFF) and normalization factor

	    NormalizationFactor = 
	      NormFactorTestHisto -> GetBinContent(SuccessfulThetaBinsVector[i]);


	    TString TmpAlphaHistoName = ("TestAlphaAfterCuts");
	    TmpAlphaHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    *fLog << "Getting Histogram with name " << endl
		  << TmpAlphaHistoName << endl;

	    AlphaONAfterCuts = (TH1F*) rootfile.Get(TmpAlphaHistoName);

	    
	    TString TmpAlphaOFFHistoName = ("TestAlphaOFFAfterCuts");
	    TmpAlphaOFFHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    *fLog << "Getting Histogram with name " << endl
		  << TmpAlphaOFFHistoName << endl;

	    AlphaOFFAfterCuts = (TH1F*) rootfile.Get(TmpAlphaOFFHistoName);
	    

	 // 2) Computing Nex and Signficance vs alpha and filling histograms
	    // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF  


	    *fLog  << endl << endl
		   << "Starting loop to compute Nex and Significance vs alphasig " << endl
		   << "for TRAIN sample and " 
		   << fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1] 
		   << " :" << endl;
	    
	   	    
	    for (Int_t j = 1; j < NexSignVSAlphaXBinsVectorDimension; j++)
	      {
		

		alphasig = XBinsVector[j];
		
		MHFindSignificanceONOFF findsig;
		findsig.SetRebin(RebinHistograms);
		findsig.SetReduceDegree(ReduceDegree);	
		findsig.SetUseFittedQuantities(fUseFittedQuantities);

		// Dummy name for psfile 
		TString psfilename = ("DummyPs"); 
	    
		findsig.FindSigmaONOFF(AlphaONAfterCuts, AlphaOFFAfterCuts, 
				       NormalizationFactor, 
				       alphabkgmin, alphabkgmax, 
				       degree, degreeOFF,
				       alphasig, 
				       drawpoly, fitgauss, print, 
				       saveplots, psfilename);



		Significance = double(findsig.GetSignificance());

		// so far there is not a variable that contains 
		// NexONOFF or NexONOFFFitted (in class MHFindSignificanceONOFF) 
		// according to the value of the variable fUseFittedQuantities.
		// As soon as I have some time I will implement it (with the 
		// corresponging GET member function) and will change 
		// MFindSignificanceONOFF accordingly to make use of it.
		
		// For the time being, I will survive with an "if statement".

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


		*fLog << endl << "Cut in alpha, Nex and significance : " 
		      << alphasig << ", " << Nex << ", " << Significance << endl << endl;
		
	      
		// updating histograms

		tmpdouble = alphasig - AlphaBinWidth/2.0;

		NexVSAlphaSigTestVector[i] -> Fill(tmpdouble, Nex);
		SignificanceVSAlphaSigTestVector[i] -> Fill(tmpdouble, Significance);
		
	  

		



	      }	  
	

	    // Dynamic memory allocated in this loop is released
	    
    
	    *fLog  << "Memory allocated by temporal alpha histograms "
		   << "will be released" << endl;
	    
	    delete AlphaONAfterCuts;
	    delete AlphaOFFAfterCuts;
	    
	  }








      }
    
    rootfile.Close();


   
     // *************************************************************
    // Histograms Nex and Significance vs alpha are written into 
    // root file fAlphaDistributionsRootFilename

    


     *fLog <<"Histograms containing Nex and significance vs alphasig " 
	  << " are saved into root file " << endl 
	  << fAlphaDistributionsRootFilename << endl;
    

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

 
   for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
    {

      if (fOptimizeParameters)
	{// Histograms for train sample are written to root file

	  if(NexVSAlphaSigTrainVector[i])
	    {
	      if(NexVSAlphaSigTrainVector[i]-> GetEntries() > 0.5)
		NexVSAlphaSigTrainVector[i]-> Write();

	    }

	   if(SignificanceVSAlphaSigTrainVector[i])
	    {
	      if(SignificanceVSAlphaSigTrainVector[i]-> GetEntries() > 0.5)
		SignificanceVSAlphaSigTrainVector[i]-> Write();
	      
	    }
	}



      if (fTestParameters)
	{// Histograms for test sample are written to root file

	  if(NexVSAlphaSigTestVector[i])
	    {
	      if(NexVSAlphaSigTestVector[i]-> GetEntries() > 0.5)
		NexVSAlphaSigTestVector[i]-> Write();

	    }

	   if(SignificanceVSAlphaSigTestVector[i])
	    {
	      if(SignificanceVSAlphaSigTestVector[i]-> GetEntries() > 0.5)
		SignificanceVSAlphaSigTestVector[i]-> Write();
	      
	    }
	}
	  


    }
    

    rootfiletowrite.Close();

 


    // Dynamic memory allocated is released 
    
    delete  XBinsVector; 

    // function finished successfully... hopefully...
    return kTRUE;

}



// Function that gets the histograms with the alpha distributions
// (for all the theta bins specified by fCosThetaRangeVector)
// stored in fAlphaDistributionsRootFilename, and combine them
// (correcting OFF histograms with the normalization factors stored 
// in NormFactorTrain or NormFactorTest) to get one single 
// Alpha distribution for ON and another one for OFF. 
// Then these histograms are given as arguments to 
// the function MHFindSignificanceONOFF::FindSigmaONOFF, 
// (Object of this class is created) to compute the 
// Overall Excess events and  significance, that will be 
// stored in variables fOverallNexTrain/Test and fOverallSigmaLiMaTrain/Test




Bool_t MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance(Bool_t CombineTrainData,
								   Bool_t CombineTestData)
{

    // First of all let's check that we have the proper the "drinks" for the party

    if (fAlphaDistributionsRootFilename.IsNull())
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl
	      << "Root file containing all alpha distributions "
	      <<" (fAlphaDistributionsRootFilename) is not well defined." << endl
	      << "Execution of the function is ABORTED" << endl;
	return kFALSE;
    }


    if (fCosThetaRangeVector.GetSize() < 2)
	
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl
	      << "Dimension of vector fCosThetaRangeVector is smaller than 2." 
	      << "So, vector fCosThetaRangeVector is NOT properly defined" << endl
	      << "Function execution will be ABORTED" << endl;

	return kFALSE;
    }

    
    TArrayI SuccessfulThetaBinsVector; 
    // Vector containing the histogram index of the successfully 
    // optimized theta bins. 
    // It is filled with the contents of the histogram 
    // fSuccessfulThetaBinsHist retrieved from the root 
    // file fAlphaDistributionsRootFilename
    // Remember that histogram index START AT 1 !!!

    TH1F* CombinedAlphaTrainON = NULL;
    TString CombinedAlphaTrainONName = ("CombinedAlphaTrainON");
    TH1F* CombinedAlphaTrainOFF = NULL;
    TString CombinedAlphaTrainOFFName = ("CombinedAlphaTrainOFF");
   
    TH1F* CombinedAlphaTestON = NULL;
    TString CombinedAlphaTestONName = ("CombinedAlphaTestON");
   
    TH1F* CombinedAlphaTestOFF = NULL;
    TString CombinedAlphaTestOFFName = ("CombinedAlphaTestOFF");
   

    

    
    Int_t NAlphaBins = 0; // Amount of bins of alpha histograms is expected to 
    // be equal for all alpha histograms

    // Vectors storing the computed error for each of the NAlphaBins 
    // Only 2 vectors are needed for OFF Train and OFF Test sample, whose 
    // alpha distributions are corrected with normalization factors. 

    // The ON sample is just added (no normalization factors are applied) 
    // and therefore, the error of the bin is the assumed by default 
    // (i.e Sqrt(BinContent)

    TArrayD ErrorInCombinedAlphaTrainOFF;
    TArrayD ErrorInCombinedAlphaTestOFF;

    TArrayD CombinedAlphaTrainOFFWithoutNormalization;
    TArrayD CombinedAlphaTestOFFWithoutNormalization;
    


    const Double_t SmallQuantity = 0.01; // This quantity will be used as the 
    // maximum difference between quantities (double) that are expected to be equal, 
    // as the center of the bins in the several alpha distributions 

    TH1F* NormFactorTrainHisto;
    TString NormFactorTrainHistoName = ("NormFactorTrainHist");
    TH1F* NormFactorTestHisto;
    TString NormFactorTestHistoName = ("NormFactorTestHist");
    TString SuccessfulThetaBinsHistoName = ("SuccessfulThetaBinsHist");
    

    TH1F* SuccessfulThetaBinsHisto;

    TH1F* TmpTH1FHisto;

    
    // Retrieve number of theta bins that have to be combined
    // from TArrayD fCosThetaRangeVector

    Int_t tmpdim = fCosThetaRangeVector.GetSize() - 1;
    const Int_t NThetaBins = tmpdim;

    // Variables to store the mean normalization factor of the 
    // combined OFF sample. 
    // MeanNormFactor = Sum (NormFactor_i*Noff_i)/Sum (Noff_i)
    

    Double_t MeanNormFactorTrain = 0.0;
    Double_t MeanNormFactorTest = 0.0;
    
    Double_t TotalNumberOFFTrain = 0.0;
    Double_t TotalNumberOFFTest = 0.0;
    


    Double_t EventCounter = 0.0;// variable used to count events used to fill histograms
    // it is double (and not Int) because the contents of an histogram bin might not 
    // be an integer value... due to normalization factors !!
    Int_t tmpint = 0;
    Double_t tmpdouble = 0.0; // tmp variable double
    Double_t tmpdouble2 = 0.0; // tmp variable double
    Double_t tmperror = 0.0; // quantity used to store 
    // temporaly the bin error in OFF histogram
    Double_t tmpnormfactor = 0.0; // variable used to store normalization factor 
    // of the theta bin i (used for TRAIN and TEST sample)
             



    // Open root file

    

    
    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl
	  << "Opening root file " << fAlphaDistributionsRootFilename << endl;



    TFile rootfile (fAlphaDistributionsRootFilename, "READ");

    // Retrieving SuccessfulThetaBinsHisto from root file. 
    

    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" << endl
	  << "Retrieving SuccessfulThetaBinsHisto from root file and filling vector "
	  << " SuccessfulThetaBinsVector" << endl;


    cout << "Getting histogram " << SuccessfulThetaBinsHistoName << endl;

    SuccessfulThetaBinsHisto = (TH1F*) rootfile.Get(SuccessfulThetaBinsHistoName);
    


    // Check that bins in this histo correspond with the 
    // CosTheta intervals defined by vector fCosThetaRangeVector

	
    if (SuccessfulThetaBinsHisto -> GetNbinsX() != NThetaBins)
    {
	*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
	      << endl
	      << "Number of theta bins defined by vector fCosThetaRangeVector (" 
	      << NThetaBins << ") does not match with the bins in histogram "
	      << SuccessfulThetaBinsHistoName << "(" 
	      << SuccessfulThetaBinsHisto -> GetNbinsX() << ")" << endl
	      << "Function execution will be ABORTED." << endl;
	return kFALSE;
    }


    // Filling vector SuccessfulThetaBinsVector

    tmpint = 0;
    for (Int_t i = 0; i < NThetaBins; i++)
    {
	if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5)
	{
	    tmpint++;
	}
    }

    SuccessfulThetaBinsVector.Set(tmpint);
    
    tmpint = 0;
    for (Int_t i = 0; i < NThetaBins; i++)
    {
	if ((SuccessfulThetaBinsHisto -> GetBinContent(i+1)) > 0.5)
	{
	    if(tmpint < SuccessfulThetaBinsVector.GetSize())
	    {
		SuccessfulThetaBinsVector[tmpint] = i+1;
		tmpint++;
	    }
	    else
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		      << endl
		      << "Problem when filling vector 'SuccessfulThetaBinsVector'" 
		      << endl
		      << "Function execution will be ABORTED." << endl;
		return kFALSE;
	    }
	    

	}
    }


   


    // *********  HISTOGRAMS FOR TRAIN SAMPLE WILL BE FILLED *********///

    if (CombineTrainData)
    {

	// Retrieving NormFactorTrainHisto from root file. 
	
	NormFactorTrainHisto = (TH1F*) rootfile.Get(NormFactorTrainHistoName);
	// NormFactorTrainHisto-> DrawCopy();

	// Check that bins in this histo correspond with the 
	// CosTheta intervals defined by vector fCosThetaRangeVector

	
	if (NormFactorTrainHisto -> GetNbinsX() != NThetaBins)
	{
	    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		  << endl
		  << "Number of theta bins defined by vector fCosThetaRangeVector (" 
		  << NThetaBins << ") does not match with the bins in histogram "
		  << NormFactorTrainHistoName << "(" 
		  << NormFactorTrainHisto -> GetNbinsX() << ")" << endl
		  << "Function execution will be ABORTED." << endl;
	    return kFALSE;
	}

	// histo test 
	/*
	cout << "NORMFACTORTRAINHIST  TEST: BinCenter and Value" << endl;
	for (Int_t k = 0; k < NThetaBins; k++)
	{
	    
	    cout << NormFactorTrainHisto -> GetBinCenter(k+1) 
		 << ";  " << NormFactorTrainHisto -> GetBinContent(k+1)<< endl;
	}
	*/

	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	{
	    tmpdouble = NormFactorTrainHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]);
	    if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble || 
		fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble)
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		      << endl
		      << "Bins defined by vector fCosThetaRangeVector " 
		      << "do not match with the bins of histogram "
		      << NormFactorTrainHistoName << endl
		      << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1]
		      << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl
		      << "NormFactorHist Bin: " <<  tmpdouble << endl
		      << "Function execution will be ABORTED." << endl;

		return kFALSE;
	    }
	    
	}


	// Let's summ up all ON alpha distributions
	
	EventCounter = 0;

	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	{
	    if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		      << endl
		      << "Component " << SuccessfulThetaBinsVector[i]-1 
		      << " of fThetaRangeStringVector is EMPTY"
		      << endl
		      << "Function execution will be ABORTED." << endl;
		return kFALSE;

		
	    }

	

	    TString TmpHistoName = ("TrainAlphaAfterCuts");
	    TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    *fLog << "Getting Histogram with name " << endl
		  << TmpHistoName << endl;

	    if (i == 0)
	    {      
		CombinedAlphaTrainON = (TH1F*) rootfile.Get(TmpHistoName);
		NAlphaBins = CombinedAlphaTrainON -> GetNbinsX();

		// Name of histogram is set
		CombinedAlphaTrainON -> SetName(CombinedAlphaTrainONName);

		// update EventCounter
		EventCounter += CombinedAlphaTrainON -> GetEntries();

		// tmp
		cout << "NAlphaBins in histo CombinedAlphaTrainON: "
		     << NAlphaBins << endl;
		// endtmp
	    }
	    else
	    {
		TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName);
		
		// update EventCounter
		EventCounter += TmpTH1FHisto -> GetEntries();
		
		for (Int_t j = 1; j <= NAlphaBins; j++)
		{
		    // At some time alpha bins might not be the same for 
		    // the different alpha distributions. That's why
		    // I will check it before summing the alpha histo contents
		    tmpdouble =  CombinedAlphaTrainON->GetBinCenter(j) -
			TmpTH1FHisto->GetBinCenter(j);
		    tmpdouble = TMath::Abs(tmpdouble);
		    

		    if (tmpdouble > SmallQuantity)
		    {
			*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
			      << endl
			      << "Bins among the several alpha ON histograms "
			      << "for TRAIN sample do not match"
			      << endl
			      << "Function execution will be ABORTED." << endl;
			return kFALSE;
		    }
		    
		    tmpdouble = CombinedAlphaTrainON->GetBinContent(j) + 
			TmpTH1FHisto->GetBinContent(j);

		    
		    CombinedAlphaTrainON-> SetBinContent(j,tmpdouble);

		}
		
		// Dynamic memory allocated for TmpTH1FHisto is released


		//tmp
		cout << "Memory allocated by temporal histogram TmpTH1FHisto "
		     << "will be released" << endl;
		//endtmp
		delete TmpTH1FHisto;
	    
	    }


	}

	// Entries of histogram CombinedAlphaTrainON is set to 
	// the events counted in EventCounter

	tmpint = int (EventCounter);
	CombinedAlphaTrainON->SetEntries(tmpint);

	// Let's summ up all OFF alpha distributions
	
	EventCounter = 0;

	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	{
	    if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		      << endl
		      << "Component " << SuccessfulThetaBinsVector[i]-1 
		      << " of fThetaRangeStringVector is EMPTY"
		      << endl
		      << "Function execution will be ABORTED." << endl;
		return kFALSE;

		
	    }

	    TString TmpHistoName = ("TrainAlphaOFFAfterCuts");
	    TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    // Normalization constant for this theta bin is stored now in tmpnormfactor
	    tmpnormfactor = NormFactorTrainHisto -> GetBinContent(SuccessfulThetaBinsVector[i]);

	    // temp 
	    cout << "Normalization factor for bin "
		 << SuccessfulThetaBinsVector[i] << " is "
		 << tmpnormfactor << endl;

	    // endtemp

	    if (i == 0)
	    {
		CombinedAlphaTrainOFF = (TH1F*) rootfile.Get(TmpHistoName);
		NAlphaBins = CombinedAlphaTrainOFF -> GetNbinsX();
		

		// Name of histogram is set
		CombinedAlphaTrainOFF -> SetName(CombinedAlphaTrainOFFName);

		// Event counter is updated taking into account the 
		// normalization factor

		EventCounter += tmpnormfactor * CombinedAlphaTrainOFF -> GetEntries();

		// Contribution to total number of OFF events and mean normfactor
		// of this theta bin is computed

		TotalNumberOFFTrain = CombinedAlphaTrainOFF -> GetEntries();
		MeanNormFactorTrain = tmpnormfactor *
		    CombinedAlphaTrainOFF -> GetEntries();
		    
		
		
		// ****
		// Dimension of ErrorInCombinedAlphaTrainOFF is set
		ErrorInCombinedAlphaTrainOFF.Set(NAlphaBins);
		
		CombinedAlphaTrainOFFWithoutNormalization.Set(NAlphaBins);

		// Histogram now is normalized (to the respective ON data)
		// by using the normalization constants stored in 
		// histogram NormFactorTrainHisto

		// Error of the normalized bins is computed and stored in 
		// vector ErrorInCombinedAlphaTrainOFF
		
		for (Int_t j = 1; j <= NAlphaBins; j++)
		{

		    // Number of events (without normalization) are 
		    // stored in vector CombinedAlphaTrainOFFWithoutNormalization

			CombinedAlphaTrainOFFWithoutNormalization[j-1] = 
			    CombinedAlphaTrainOFF -> GetBinContent(j);
		   
		    // Bin content is set
		    tmpdouble2 = tmpnormfactor * 
			CombinedAlphaTrainOFF -> GetBinContent(j);

		    CombinedAlphaTrainOFF -> SetBinContent(j,tmpdouble2);

		    // (Bin Error)^2 is computed and stored in 
		    // ErrorInCombinedAlphaTrainOFF
		    // Bin error = Sqrt(BinContent) X Normalization Factor

		    tmperror = tmpdouble2*tmpnormfactor;
		    ErrorInCombinedAlphaTrainOFF[j-1] = tmperror;

		}

	    }
	    else
	    {
		TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName);


		// Event counter is updated taking into account the 
		// normalization factor

		EventCounter += tmpnormfactor * TmpTH1FHisto -> GetEntries();

		// Contribution to total number of OFF events and mean normfactor
		// of this theta bin is computed

		TotalNumberOFFTrain += TmpTH1FHisto -> GetEntries();
		MeanNormFactorTrain += tmpnormfactor *
		    TmpTH1FHisto -> GetEntries();

		

		for (Int_t j = 1; j <= NAlphaBins; j++)
		{
		    // At some time alpha bins might not be the same for 
		    // the different alpha distributions. That's why
		    // I will check it before summing the alpha histo contents
		    tmpdouble2 =  CombinedAlphaTrainOFF->GetBinCenter(j) -
			TmpTH1FHisto->GetBinCenter(j);
		    tmpdouble2 = TMath::Abs(tmpdouble2);

		    if (tmpdouble2 > SmallQuantity)
		    {
			*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
			      << endl
			      << "Bins among the several alpha OFF histograms "
			      << "for TRAIN sample do not match"
			      << endl
			      << "Function execution will be ABORTED." << endl;
			return kFALSE;
		    }
		    

		     // Number of events (without normalization) are 
		    // added in vector CombinedAlphaTrainOFFWithoutNormalization
		    
		    CombinedAlphaTrainOFFWithoutNormalization[j-1] += 
			TmpTH1FHisto -> GetBinContent(j);

		    // Histogram bin contents must be normalized (to the respective ON data)
		    // by using the normalization constants stored in 
		    // histogram NormFactorTrainHisto before being added 
		    // to CombinedAlphaTrainOFF

		   
		    tmpdouble2 = CombinedAlphaTrainOFF->GetBinContent(j) + 
			(TmpTH1FHisto->GetBinContent(j) * tmpnormfactor);

		    
		    CombinedAlphaTrainOFF-> SetBinContent(j,tmpdouble2);


		    // (Bin Error)^2 is computed and added to  
		    // ErrorInCombinedAlphaTrainOFF 
		    // Bin error = Sqrt(BinContent) X Normalization Factor

		    tmperror = TmpTH1FHisto->GetBinContent(j) * tmpnormfactor;
		    ErrorInCombinedAlphaTrainOFF[j-1] += tmperror;

		}
		
		// Dynamic memory allocated for TmpTH1FHisto is released

		delete TmpTH1FHisto;
	    
	    }

	    



	}

	


	// Mean Normalization factor is computed for the combined OFF histogram. 
	// This factor will be used when computing the significance via 
	// MHFindSignificance::FindSigmaONOFF().

	// The bin contents (and errors) of teh combined OFF TRAIN histo will be 
	// normalized (divided) by teh mean norm constant, so that 
	// the histogram for OFF data given as argument in FindSigmaONOFF
	// is equivalent to a non-normalized histogram. 
	
	MeanNormFactorTrain = MeanNormFactorTrain/TotalNumberOFFTrain;



	// Sqrt(Contents) is performed in ErrorInCombinedAlphaTrainOFF, 
	// in order to compute the real error in the bin content of histogram 
	// CombinedAlphaTrainOFF.

	// Then error and contents (normalized with mean normfctor) are
	// set to histogram bin 	

	for (Int_t j = 1; j <= NAlphaBins; j++)
	{	    
	    tmpdouble2 = 
		(CombinedAlphaTrainOFF -> GetBinContent(j))/MeanNormFactorTrain;

	    CombinedAlphaTrainOFF -> SetBinContent(j,tmpdouble2);

	    ErrorInCombinedAlphaTrainOFF[j-1] = 
		TMath::Sqrt(ErrorInCombinedAlphaTrainOFF[j-1])/MeanNormFactorTrain;
	    // Proper error treatment when adding quantities to histogram bin 

	    CombinedAlphaTrainOFF -> SetBinError(j,ErrorInCombinedAlphaTrainOFF[j-1]);

	    // ******  TMP **************
	    /*
	    // Error computation assuming a given (no normalized) bin content

	    tmpdouble = CombinedAlphaTrainOFF -> GetBinContent(j);

	    if (CombinedAlphaTrainOFFWithoutNormalization[j-1] < SmallQuantity)
	    {tmpnormfactor = 0;}
	    else
	    {tmpnormfactor = tmpdouble/CombinedAlphaTrainOFFWithoutNormalization[j-1];}

	    tmperror = TMath::Sqrt(CombinedAlphaTrainOFFWithoutNormalization[j-1]) *
		tmpnormfactor;

	    CombinedAlphaTrainOFF -> SetBinError(j,tmperror);
	    
	    */
	    // ******  END TMP **************			 

	    // tmp
	    /*
	    cout << CombinedAlphaTrainOFF -> GetBinContent(j) << " +/-  " 
		 << CombinedAlphaTrainOFF -> GetBinError(j)
		 << "  Number Events without normalization: " 
		 << CombinedAlphaTrainOFFWithoutNormalization[j-1] << endl;
	    */
	    // endtmp
	    
	}


	// Entries of histogram CombinedAlphaTrainON is set to 
	// the events counted during histogram filling

	EventCounter = EventCounter/MeanNormFactorTrain;
	tmpint = int (EventCounter);
	CombinedAlphaTrainOFF->SetEntries(tmpint);


	cout << "SILLY INFO FOR TRAIN SAMPLE: Total number of events Before nomralization and " 
	     << " re-re normalized" << endl
	     << TotalNumberOFFTrain << " - " << EventCounter << endl;

    }

    // *********  COMBINED  HISTOGRAMS FOR TRAIN SAMPLE ALREADY FILLED ********* ///




    // *********  HISTOGRAMS FOR TEST SAMPLE WILL BE FILLED *********  ///


    if (CombineTestData)
    {

	// Retrieving NormFactorTestHisto from root file. 
	
	NormFactorTestHisto = (TH1F*) rootfile.Get(NormFactorTestHistoName);
	// NormFactorTestHisto-> DrawCopy();

	// Check that bins in this histo correspond with the 
	// CosTheta intervals defined by vector fCosThetaRangeVector

	
	if (NormFactorTestHisto -> GetNbinsX() != NThetaBins)
	{
	    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		  << endl
		  << "Number of theta bins defined by vector fCosThetaRangeVector (" 
		  << NThetaBins << ") does not match with the bins in histogram "
		  << NormFactorTestHistoName << "(" 
		  << NormFactorTestHisto -> GetNbinsX() << ")" << endl
		  << "Function execution will be ABORTED." << endl;
	    return kFALSE;
	}

	// histo test 
	/*
	cout << "NORMFACTORTRAINHIST  TEST: BinCenter and Value" << endl;
	for (Int_t k = 0; k < NThetaBins; k++)
	{
	    
	    cout << NormFactorTestHisto -> GetBinCenter(k+1) 
		 << ";  " << NormFactorTestHisto -> GetBinContent(k+1)<< endl;
	}
	*/

	

	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	{
	    tmpdouble = NormFactorTestHisto -> GetBinCenter(SuccessfulThetaBinsVector[i]);
	    if (fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1] > tmpdouble || 
		fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] < tmpdouble)
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		      << endl
		      << "Bins defined by vector fCosThetaRangeVector " 
		      << "do not match with the bins of histogram "
		      << NormFactorTestHistoName << endl
		      << "CosTheta: " << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]-1]
		      << "-" << fCosThetaRangeVector[SuccessfulThetaBinsVector[i]] << endl
		      << "NormFactorHist Bin: " <<  tmpdouble << endl
		      << "Function execution will be ABORTED." << endl;

		return kFALSE;
	    }
	    
	}


	// Let's summ up all ON alpha distributions

	// Event counter (counts events used in the histogram filling) is 
	// initialized to zero

	EventCounter = 0;

	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	{
	    if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		      << endl
		      << "Component " << SuccessfulThetaBinsVector[i]-1 
		      << " of fThetaRangeStringVector is EMPTY"
		      << endl
		      << "Function execution will be ABORTED." << endl;
		return kFALSE;

		
	    }

	

	    TString TmpHistoName = ("TestAlphaAfterCuts");
	    TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    *fLog << "Getting Histogram with name " << endl
		  << TmpHistoName << endl;

	    if (i == 0)
	    {      
		CombinedAlphaTestON = (TH1F*) rootfile.Get(TmpHistoName);
		NAlphaBins = CombinedAlphaTestON -> GetNbinsX();

		// Name of histogram is set
		CombinedAlphaTestON -> SetName(CombinedAlphaTestONName);

		// EventCounter is updated
		
		EventCounter += CombinedAlphaTestON -> GetEntries();

		// tmp
		cout << "NAlphaBins in histo CombinedAlphaTestON: "
		     << NAlphaBins << endl;
		// endtmp
	    }
	    else
	    {
		TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName);

		// EventCounter is updated
		
		EventCounter += TmpTH1FHisto -> GetEntries();

		for (Int_t j = 1; j <= NAlphaBins; j++)
		{
		    // At some time alpha bins might not be the same for 
		    // the different alpha distributions. That's why
		    // I will check it before summing the alpha histo contents
		    tmpdouble =  CombinedAlphaTestON->GetBinCenter(j) -
			TmpTH1FHisto->GetBinCenter(j);
		    tmpdouble = TMath::Abs(tmpdouble);
		    

		    if (tmpdouble > SmallQuantity)
		    {
			*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
			      << endl
			      << "Bins among the several alpha ON histograms "
			      << "for TEST sample do not match"
			      << endl
			      << "Function execution will be ABORTED." << endl;
			return kFALSE;
		    }
		    
		    tmpdouble = CombinedAlphaTestON->GetBinContent(j) + 
			TmpTH1FHisto->GetBinContent(j);

		    
		    CombinedAlphaTestON-> SetBinContent(j,tmpdouble);

		}
		
		// Dynamic memory allocated for TmpTH1FHisto is released


		//tmp
		cout << "Memory allocated by temporal histogram TmpTH1FHisto "
		     << "will be released" << endl;
		//endtmp
		delete TmpTH1FHisto;
	    
	    }




	}

	// Entries of histogram CombinedAlphaTestON is set to 
	// the events counted during histogram filling
	tmpint = int (EventCounter);
	CombinedAlphaTestON->SetEntries(tmpint);


	// Let's summ up all OFF alpha distributions

	EventCounter = 0;

	for (Int_t i = 0; i < SuccessfulThetaBinsVector.GetSize(); i++)
	{
	    if (fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1].IsNull())
	    {
		*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
		      << endl
		      << "Component " << SuccessfulThetaBinsVector[i]-1 
		      << " of fThetaRangeStringVector is EMPTY"
		      << endl
		      << "Function execution will be ABORTED." << endl;
		return kFALSE;

		
	    }

	    TString TmpHistoName = ("TestAlphaOFFAfterCuts");
	    TmpHistoName += fThetaRangeStringVector[SuccessfulThetaBinsVector[i]-1];

	    // Normalization constant for this theta bin is stored now in tmpnormfactor
	    tmpnormfactor = NormFactorTestHisto -> GetBinContent(SuccessfulThetaBinsVector[i]);

	    // temp 
	    cout << "Normalization factor for bin "
		 << SuccessfulThetaBinsVector[i] << " is "
		 << tmpnormfactor << endl;

	    // endtemp

	    if (i == 0)
	    {
		CombinedAlphaTestOFF = (TH1F*) rootfile.Get(TmpHistoName);
		NAlphaBins = CombinedAlphaTestOFF -> GetNbinsX();
	       

		// Name of histogram is set
		CombinedAlphaTestOFF -> SetName(CombinedAlphaTestOFFName);

		// EventCounter is updated taking into account 
		// normalization factor
		
		EventCounter += tmpnormfactor*CombinedAlphaTestOFF -> GetEntries();

		
		// Contribution to total number of OFF events and mean normfactor
		// of this theta bin is computed

		TotalNumberOFFTest = CombinedAlphaTestOFF -> GetEntries();
		MeanNormFactorTest = tmpnormfactor *
		    CombinedAlphaTestOFF -> GetEntries();


		// Dimension of ErrorInCombinedAlphaTrainOFF is set
		ErrorInCombinedAlphaTestOFF.Set(NAlphaBins);

		CombinedAlphaTestOFFWithoutNormalization.Set(NAlphaBins);


		// Histogram now is normalized (to the respective ON data)
		// by using the normalization constants stored in 
		// histogram NormFactorTestHisto

		
		for (Int_t j = 1; j <= NAlphaBins; j++)
		{

		    // Vector containing number of events without 
		    // normalization is filled
		   
			CombinedAlphaTestOFFWithoutNormalization[j-1] = 
			    CombinedAlphaTestOFF -> GetBinContent(j);
		   


		    tmpdouble2 = tmpnormfactor * 
			CombinedAlphaTestOFF -> GetBinContent(j);

		    CombinedAlphaTestOFF -> SetBinContent(j,tmpdouble2);

		    // (Bin Error)^2 is computed and stored in 
		    // ErrorInCombinedAlphaTestOFF
		    // Bin error = Sqrt(BinContent) X Normalization Factor

		    tmperror = tmpdouble2*tmpnormfactor;
		    ErrorInCombinedAlphaTestOFF[j-1] = tmperror;


		}

	    }
	    else
	    {
		TmpTH1FHisto = (TH1F*) rootfile.Get(TmpHistoName);


		// EventCounter is updated taking into account 
		// normalization factor
		
		EventCounter += tmpnormfactor* TmpTH1FHisto-> GetEntries();

		// Contribution to total number of OFF events and mean normfactor
		// of this theta bin is computed


		TotalNumberOFFTest += TmpTH1FHisto -> GetEntries();
		MeanNormFactorTest += tmpnormfactor *
		    TmpTH1FHisto -> GetEntries();



		for (Int_t j = 1; j <= NAlphaBins; j++)
		{
		    // At some time alpha bins might not be the same for 
		    // the different alpha distributions. That's why
		    // I will check it before summing the alpha histo contents
		    tmpdouble2 =  CombinedAlphaTestOFF->GetBinCenter(j) -
			TmpTH1FHisto->GetBinCenter(j);
		    tmpdouble2 = TMath::Abs(tmpdouble2);

		    if (tmpdouble2 > SmallQuantity)
		    {
			*fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance()" 
			      << endl
			      << "Bins among the several alpha OFF histograms "
			      << "for TRAIN sample do not match"
			      << endl
			      << "Function execution will be ABORTED." << endl;
			return kFALSE;
		    }
		    

		    // Vector containing number of events without 
		    // normalization is updated

		    // Vector containing number of events without 
		    // normalization is filled
		   
			CombinedAlphaTestOFFWithoutNormalization[j-1] += 
			    TmpTH1FHisto -> GetBinContent(j);
		   



		    
		    
		    // Histogram bin contents must be normalized (to the respective ON data)
		    // by using the normalization constants stored in 
		    // histogram NormFactorTestHisto before being added 
		    // to CombinedAlphaTestOFF

		   
		    tmpdouble2 = CombinedAlphaTestOFF->GetBinContent(j) + 
			(TmpTH1FHisto->GetBinContent(j) * tmpnormfactor);

		    
		    CombinedAlphaTestOFF-> SetBinContent(j,tmpdouble2);


		     // (Bin Error)^2 is computed and added to  
		    // ErrorInCombinedAlphaTestOFF 
		    // Bin error = Sqrt(BinContent) X Normalization Factor

		    tmperror = TmpTH1FHisto->GetBinContent(j) * tmpnormfactor;
		    ErrorInCombinedAlphaTestOFF[j-1] += tmperror;


		}
		
		// Dynamic memory allocated for TmpTH1FHisto is released

		delete TmpTH1FHisto;
	    
	    }

	    



	}



	// Mean Normalization factor is computed for the combined OFF histogram. 
	// This factor will be used when computing the significance via 
	// MHFindSignificance::FindSigmaONOFF().

	// The bin contents (and errors) of teh combined OFF TEST histo will be 
	// normalized (divided) by teh mean norm constant, so that 
	// the histogram for OFF data given as argument in FindSigmaONOFF
	// is equivalent to a non-normalized histogram. 
	
	MeanNormFactorTest = MeanNormFactorTest/TotalNumberOFFTest;



	// Sqrt(Contents) is performed in ErrorInCombinedAlphaTestOFF, 
	// in order to compute the real error in the bin content of histogram 
	// CombinedAlphaTestOFF.

	// Then error and contents (normalized with mean normfctor) are
	// set to histogram bin 	




	
	for (Int_t j = 1; j <= NAlphaBins; j++)
	{

	    tmpdouble2 = 
		(CombinedAlphaTestOFF -> GetBinContent(j))/MeanNormFactorTest;

	    CombinedAlphaTestOFF -> SetBinContent(j,tmpdouble2);

	    ErrorInCombinedAlphaTestOFF[j-1] = 
		TMath::Sqrt(ErrorInCombinedAlphaTestOFF[j-1])/MeanNormFactorTest;

	    // Proper error treatment for the bin contents 
	    CombinedAlphaTestOFF -> SetBinError(j,ErrorInCombinedAlphaTestOFF[j-1]);
	    
	    
	    // ******  TMP **************
	    /*
	    // Error computation assuming a given (no normalized) bin content
	    
	    tmpdouble = CombinedAlphaTestOFF -> GetBinContent(j);

	    if (CombinedAlphaTestOFFWithoutNormalization[j-1] < SmallQuantity)
	    {tmpnormfactor = 0;}
	    else
	    {tmpnormfactor = tmpdouble/CombinedAlphaTestOFFWithoutNormalization[j-1];}
	   
	    tmperror = (TMath::Sqrt(CombinedAlphaTestOFFWithoutNormalization[j-1])) * tmpnormfactor;
	    
	    CombinedAlphaTestOFF -> SetBinError(j,tmperror);
	    
	    */
	    
	    // ******  ENDTMP **************
	    
	    // tmp
	    // cout << CombinedAlphaTestOFF -> GetBinContent(j) << " +/-  " 
	    //	 << CombinedAlphaTestOFF -> GetBinError(j) << endl;
	    
	    // endtmp

	}


	
	// Entries of histogram CombinedAlphaTestOFF is set to 
	// the events counted during histogram filling
	
	EventCounter = EventCounter/MeanNormFactorTest;
	tmpint = int(EventCounter);
	CombinedAlphaTestOFF->SetEntries(tmpint);


	cout << "SILLY INFO FOR TEST SAMPLE: Total number of events Before nomralization and " 
	     << " re-re normalized" << endl
	     << TotalNumberOFFTest << " - " << EventCounter << endl;

	/*
	// Sumw2() test

	cout << "********   TEST   ********* " << endl;
	cout << "Sumw2() is executed: " << endl;
	CombinedAlphaTestOFF -> Sumw2();
	for (Int_t j = 1; j <= NAlphaBins; j++)
	{
	     // tmp
	    cout << CombinedAlphaTestOFF -> GetBinContent(j) << " +/-  " 
		 << CombinedAlphaTestOFF -> GetBinError(j) << endl;
	    
	    // endtmp

	}
	  
	*/
	
	

    }



     // *********  COMBINED  HISTOGRAMS FOR TEST SAMPLE ALREADY FILLED ********* ///


    /*
    CombinedAlphaTrainON-> DrawCopy();
    gPad -> SaveAs("tmpONTrain.ps");

    CombinedAlphaTrainOFF-> DrawCopy();
    gPad -> SaveAs("tmpOFFTrain.ps");


    CombinedAlphaTestON-> DrawCopy();
    gPad -> SaveAs("tmpONTest.ps");

    CombinedAlphaTestOFF-> DrawCopy();
    gPad -> SaveAs("tmpOFFTest.ps");

    */
    

    // ****  COMPUTATION OF SIGNIFICANCES FOR COMBINED ALPHA HISTOGRAMS  ****** ///
    
    if (SuccessfulThetaBinsVector.GetSize() >= 1)
    {
	// There is, at least, ONE theta bin for which optimization of 
	// supercuts was SUCCESSFUL, and thus, it is worth to 
	// compute significance of signal

		
	// Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF
	

	
	// Maximum value of alpha below which signal is expected
	const Double_t alphasig = fAlphaSig; 

	// Minimum value of alpha for bkg region in ON data
	const Double_t alphabkgmin = fAlphaBkgMin;
	
	// Maximum value of alpha for bkg region in ON data
	const Double_t alphabkgmax = fAlphaBkgMax;
	
	// degree of polynomial function used to fit ON data in background region
	const Int_t    degree   =    2;

	// degree of polynomial function used to fit OFF data in all alpha region (0-90)
	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 final alpha plot for ON and OFF 
	// (with Nex and significance) in Psfile
   



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

	

	if (CombineTrainData)
	{

	    // Name of psfile where Alpha plot will be stored
	    TString psfilename = ("TRAINSampleCombined"); 
	    
	    findsig.FindSigmaONOFF(CombinedAlphaTrainON, CombinedAlphaTrainOFF, 
				   MeanNormFactorTrain, 
				   alphabkgmin, alphabkgmax, 
				   degree, degreeOFF,
				   alphasig, 
				   drawpoly, fitgauss, print, 
				   saveplots, psfilename);


	  fOverallSigmaLiMaTrain = double(findsig.GetSignificance());
	  fOverallNexTrain = double(findsig.GetNexONOFF());

	  *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl
		<< "After performing the combined analysis for TRAIN sample it was found " 
		<< "the following Nex and Significance: " <<  fOverallNexTrain 
		<< ", " << fOverallSigmaLiMaTrain << endl;
	      

	}


	if (CombineTestData)
	{

	    // Name of psfile where Alpha plot will be stored
	    TString psfilename = ("TESTSampleCombined"); 
	    
	    findsig.FindSigmaONOFF(CombinedAlphaTestON, CombinedAlphaTestOFF, 
				   MeanNormFactorTest, 
				   alphabkgmin, alphabkgmax, 
				   degree, degreeOFF,
				   alphasig, 
				   drawpoly, fitgauss, print, 
				   saveplots, psfilename);

	    fOverallSigmaLiMaTest = double(findsig.GetSignificance());
	    fOverallNexTest = double(findsig.GetNexONOFF());
	    
	    *fLog << "MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl
		  << "After performing the combined analysis for TEST sample it was found " 
		  << "the following Nex and Significance: " <<  fOverallNexTest 
		  << ", " << fOverallSigmaLiMaTest << endl;
	    
	}
	
	

	
    }


    return kTRUE;
}

    
    





