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

/////////////////////////////////////////////////////////////////////////////
//
// MHFindSignificanceONOFF
//
// determines the significance of a gamma signal in an |alpha| plot
// it uses real OFF events in the computation of excess events
//
// Input :  2 TH1 histogram (ON and OFF) of |alpha| : with 0 < |alpha| < 90 degrees
//  
//************TEMP ************************************************


//         alphamin, alphamax :     defining the background region
//         alphasig           :     defining the signal region for which
//                                  the significance is calculated
//         degree : the degree of the polynomial to be fitted to the background
//                  ( a0 + a1*x + a2*x**2 + a3*x**3 + ...) 
//
// Output : 
//
//
//   - the number of events in the signal region (Non)
//     the number of background events in the signal region (Nbg)
//     (counted number of events 'NoffSig' and fitted number of events 'NoffSigFitted 
//   - the number of excess events in the signal region (Nex = Non - Nbg)
//     (again, counted 'NexONOFF' and fitted 'NexONOFFFitted'
//   - thew effective number of background events (Noff), and gamma :
//     Nbg = gamma * Noff
//   - 
//   - the significance of the gamma signal according to Li & Ma               
//     3 significances are computed 
//        a) LiMa formula (17) using fitted quantities; fSigLiMa
//        b) LiMa formula (17) using counted quantities; fSigLiMa2
//        c) LiMa formula (5) using counted quantities.
//
// call member function 'FindSigmaONOFF' 
//      to fit the background and to determine the significance  
//
//
//
/////////////////////////////////////////////////////////////////////////////

////////// ***********  ENDTEMPORAL  *************************


#include "MHFindSignificanceONOFF.h"

#include <fstream>
#include <math.h>

#include <TArrayD.h>
#include <TArrayI.h>
#include <TH1.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TFitter.h>
#include <TMinuit.h>
#include <TPaveText.h>
#include <TStyle.h>
#include <TPostScript.h>

#include "MLog.h"
#include "MLogManip.h"
#include "MMinuitInterface.h"


ClassImp(MHFindSignificanceONOFF);

using namespace std;

const TString MHFindSignificanceONOFF::gsDefName  = "MHFindSignificanceONOFF";
const TString MHFindSignificanceONOFF::gsDefTitle = "Find Significance in alpha plot";



// --------------------------------------------------------------------------
//
// fcnpoly
//
// calculates the chi2 for the fit of the polynomial function 'poly' 
// to the histogram 'fhist'
//
// it is called by CallMinuit() (which is called in FitPolynomial()) 
//
// bins of fhist with huge errors are ignored in the calculation of the chi2
// (the huge errors were set in 'FitPolynomial()')
//

static void fcnpoly(Int_t &npar, Double_t *gin, Double_t &f, 
                    Double_t *par, Int_t iflag)
{
    TH1 *fhist = (TH1*)gMinuit->GetObjectFit();
    TF1 *fpoly = fhist->GetFunction("Poly");    


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

    Double_t chi2 = 0.0;

    Int_t nbins = fhist->GetNbinsX();
    Int_t mbins = 0;
    for (Int_t i=1; i<=nbins; i++)
    {
      Double_t content = fhist->GetBinContent(i);
      Double_t error   = fhist->GetBinError(i);
      Double_t center  = fhist->GetBinCenter(i);

      //-----------------------------
      // ignore unwanted points      if (error > 1.e19)
        continue;

      if (content <= 0.0)
      {
        gLog << "fcnpoly : bin with zero content; i, content, error = "
             << i << ",  " << content << ",  " << error << endl;
        continue;
      }        

      if (error <= 0.0)
      {
        gLog << "fcnpoly : bin with zero error; i, content, error = "
             << i << ",  " << content << ",  " << error << endl;
        continue;
      }        

      //-----------------------------
      mbins++;

      Double_t fu;
      fu = fpoly->EvalPar(&center, par);

      // the fitted function must not be negative
      if (fu <= 0.0)
      {
        chi2 = 1.e10;
        break;
      }

      Double_t temp = (content - fu) / error;
      chi2 += temp*temp;
    }

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

    f = chi2;

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

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


// --------------------------------------------------------------------------
//
// fcnpolyOFF
//
// calculates the chi2 for the fit of the polynomial function 'poly' 
// to the histogram 'fhist'
//
// it is called by CallMinuit() (which is called in FitPolynomial()) 
//
// bins of fhist with huge errors are ignored in the calculation of the chi2
// (the huge errors were set in 'FitPolynomial()')
//

static void fcnpolyOFF(Int_t &npar, Double_t *gin, Double_t &f, 
		       Double_t *par, Int_t iflag)
{
    TH1 *fhist = (TH1*)gMinuit->GetObjectFit();
    TF1 *fpoly = fhist->GetFunction("PolyOFF");    


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

    Double_t chi2 = 0.0;

    Int_t nbins = fhist->GetNbinsX();
    Int_t mbins = 0;
    for (Int_t i=1; i<=nbins; i++)
    {
      Double_t content = fhist->GetBinContent(i);
      Double_t error   = fhist->GetBinError(i);
      Double_t center  = fhist->GetBinCenter(i);

      //-----------------------------
      // ignore unwanted points
      if (error > 1.e19)
        continue;

      if (content <= 0.0)
      {
        gLog << "fcnpoly : bin with zero content; i, content, error = "
             << i << ",  " << content << ",  " << error << endl;
        continue;
      }        

      if (error <= 0.0)
      {
        gLog << "fcnpoly : bin with zero error; i, content, error = "
             << i << ",  " << content << ",  " << error << endl;
        continue;
      }        

      //-----------------------------
      mbins++;

      Double_t fu;
      fu = fpoly->EvalPar(&center, par);

      // the fitted function must not be negative
      if (fu <= 0.0)
      {
        chi2 = 1.e10;
        break;
      }

      Double_t temp = (content - fu) / error;
      chi2 += temp*temp;
    }

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

    f = chi2;

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

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




// --------------------------------------------------------------------------
//
// fcnpolygauss
//
// calculates the chi2 for the fit of the (polynomial+Gauss) function 
// 'PolyGauss' to the histogram 'fhist'
//
// it is called by CallMinuit() (which is called in FitGaussPoly()) 
//
// bins of fhist with huge errors are ignored in the calculation of the chi2
// (the huge errors were set in 'FitGaussPoly()')
//

static void fcnpolygauss(Int_t &npar, Double_t *gin, Double_t &f, 
                         Double_t *par, Int_t iflag)
{
    TH1 *fhist = (TH1*)gMinuit->GetObjectFit();
    TF1 *fpolygauss = fhist->GetFunction("PolyGauss");    


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

    Double_t chi2 = 0.0;

    Int_t nbins = fhist->GetNbinsX();
    Int_t mbins = 0;
    for (Int_t i=1; i<=nbins; i++)
    {
      Double_t content = fhist->GetBinContent(i);
      Double_t error   = fhist->GetBinError(i);
      Double_t center  = fhist->GetBinCenter(i);

      //-----------------------------
      // ignore unwanted points
      if (error > 1.e19)
        continue;

      if (content <= 0.0)
      {
        gLog << "fcnpolygauss : bin with zero content; i, content, error = "
             << i << ",  " << content << ",  " << error << endl;
        continue;
      }        

      if (error <= 0.0)
      {
        gLog << "fcnpolygauss : bin with zero error; i, content, error = "
             << i << ",  " << content << ",  " << error << endl;
        continue;
      }        

      //-----------------------------
      mbins++;

      Double_t fu;
      fu = fpolygauss->EvalPar(&center, par);

      // the fitted function must not be negative
      if (fu <= 0.0)
      {
        chi2 = 1.e10;
        break;
      }

      Double_t temp = (content - fu) / error;
      chi2 += temp*temp;
    }

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

    f = chi2;

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

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



// --------------------------------------------------------------------------
//
//  Constructor
//
MHFindSignificanceONOFF::MHFindSignificanceONOFF(const char *name, const char *title)
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();

//    fSigVsAlpha = NULL;

    fPoly   = NULL;
    fPolyOFF   = NULL;
    fPolyOFFNormalized   = NULL;
    fGPoly  = NULL;
    fGBackg = NULL;

    fHist     = NULL;
    fHistOrig = NULL;

    fHistOFF     = NULL;
    fHistOrigOFF = NULL;
    fHistOFFNormalized = NULL;

    // allow rebinning of the alpha plot
//    fRebin = kTRUE;

    fRebin = kFALSE;

    // allow reducing the degree of the polynomial
    fReduceDegree = kTRUE;

    // Low and Upper limits for the OFF alpha distribution fit
    // are set to 0 and 90 degrees respectively

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

    // Bool variable used to decide wether to print or not the results 
    // of the fit, significance, Nex... onto the final alpha plot. 
    // for the time being, this variable is set in the constructor. 
    // At some point, I might make it such it can be set externally...

    fPrintResultsOntoAlphaPlot = kTRUE;
    //fPrintResultsOntoAlphaPlot = kFALSE;



    fCanvas = NULL;
    
    fSavePlots = kFALSE; // By default plots are not saved in Postscriptfiles

    fPsFilename = NULL;

   

}

// --------------------------------------------------------------------------
//
//  Destructor. 
//
// =====>  it is not clear why one obtains sometimes a segmentation violation
//         when the destructor is active     <=======================
//
// therefore the 'return'statement
//

MHFindSignificanceONOFF::~MHFindSignificanceONOFF()
{


    
     return;
    
    *fLog << "destructor of MHFindSignificanceONOFF is called" << endl;
    

    fPsFilename = NULL;
  
    cout << "PS set to null... " << endl;

    delete fHist;
    delete fHistOFF;

    delete fHistOFFNormalized;

    cout << "Histos removed..." << endl;

    delete fPoly;
    delete fPolyOFF;
    delete fPolyOFFNormalized;
    delete fGPoly;
    delete fGBackg;
  
    cout << "Functions are also removed..." << endl;
    
    // delete fCanvas; if I removed fCanvas pointed memory address the 
    // program crashes  ???

    *fLog << "destructor of MHFindSignificanceONOFF finished successfully" << endl;

   
}

// --------------------------------------------------------------------------
//
//  Set flag fRebin 
//
//  if flag is kTRUE rebinning of the alpha plot is allowed
//
//
void MHFindSignificanceONOFF::SetRebin(Bool_t b)
{
  fRebin = b;

  *fLog << "MHFindSignificanceONOFF::SetRebin; flag fRebin set to " 
        << (b? "kTRUE" : "kFALSE") << endl;
}

// --------------------------------------------------------------------------
//
//  Set flag fReduceDegree 
//
//  if flag is kTRUE reducing of the degree of the polynomial is allowed
//
//
void MHFindSignificanceONOFF::SetReduceDegree(Bool_t b)
{
  fReduceDegree = b;

  *fLog << "MHFindSignificanceONOFF::SetReduceDegree; flag fReduceDegree set to " 
        << (b? "kTRUE" : "kFALSE") << endl;
}

// Function that returns one of the 3 LiMa sigmas. 
// The returned value is the one used in the optimization 
// and final alpha plots. 

// For the time being, if fUseFittedQuantities = kTRUE (default)
// fSigLiMa is used, otherwise, fSigLiMa2 is used. 

Double_t MHFindSignificanceONOFF::GetSignificance()
{

    if(fUseFittedQuantities)
    {   
	return fSigLiMa;
    }

    return fSigLiMa2;


    
}


// --------------------------------------------------------------------------
//
//  FindSigmaONOFF
//
//  calls FitPolynomialOFF  it gets bkg events in "signal" region 
//                          from histogram histOFF which is the alpha 
//                          distribution of OFF data NON normalized. 
//                          Normalization factor is also one of the 
//                          arguments.
//  calls DetExcess         to determine the number of excess events
//                          using the previously computed bkg events

//  calls SigmaLiMa         to determine the significance of the gamma signal
//                          in the range |alpha| < alphasig
//  calls FitGaussPoly      to fit a (polynomial+Gauss) function in the
//                          whole |alpha| region of ON - OFF diagram 
//
//

Bool_t MHFindSignificanceONOFF::FindSigmaONOFF(TH1 *fhistON,  
					       TH1 *fhistOFF,
					       Double_t NormFactor,
					       Double_t alphamin,
					       Double_t alphamax,
					       Int_t degreeON, 
					       Int_t degreeOFF,
					       Double_t alphasig, 
					       Bool_t drawpoly,   
					       Bool_t fitgauss, 
					       Bool_t print, 
					       Bool_t saveplots, 
					       //TPostScript* PsFile
					       const TString psfilename)
{
    
    //*fLog << "MHFindSignificanceONOFF::FindSigma;" << endl;
 
    // Pointer to object TPostScript where plots will be stored 
    // is copied into member varialbe 
    // NOT WORKING !!!
    // fPsFilename = PsFile;

    // Temporally (while TPostSctipt option is not working) psfiles 
    // will be produced by the standard way (Canvas.SaveAs())

    fPsFilenameString = psfilename;



    // "3 Li and Ma significances" are initialized to 0.0

    fSigLiMa = 0.0;
    fSigLiMa2 = 0.0;
    fSigLiMa3 = 0.0;
   
    
    // fNormFactor is set 
    
    fNormFactor = NormFactor;

    // Report when this histograms given in the 
    // arguments are empty
    
    Double_t tmpdouble= -1.0;
    
    tmpdouble = double(fhistON->GetEntries());
    if (tmpdouble < 0.5)
    {
	cout << "MHFindSignificanceONOFF::FindSigmaONOFF; ERROR " << endl
	     << "fhistON has ZERO entries" << endl;
   } 

    tmpdouble = double(fhistOFF->GetEntries());
    if (tmpdouble < 0.5)
    {
	cout << "MHFindSignificanceONOFF::FindSigmaONOFF; ERROR " << endl
	     << "fhistOFF has ZERO entries" << endl;
    } 



   
    // Variables set for alpha from ON events
    fHistOrig = fhistON;
    fHist = (TH1*)fHistOrig->Clone();
    fHist->SetName(fhistON->GetName());
    fDegree   = degreeON;


     // Variables set for alpha from OFF events

    fHistOrigOFF = fhistOFF;    
    fHistOFF = (TH1*)fHistOrigOFF->Clone();
    fHistOFF->SetName(fhistOFF->GetName());
    fDegreeOFF   = degreeOFF;

    if ( !fHist )
    {
	*fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; Clone of histogram could not be generated" 
	      << endl;
	return kFALSE;
    }


    
    if ( !fHistOFF )
    {
	*fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; Clone of OFF histogram could not be generated" 
	      << endl;
	return kFALSE;
    }



  
  fHist->Sumw2();
  fHist->SetXTitle("|alpha|  [\\circ]");
  fHist->SetYTitle("Counts");
  fHist->UseCurrentStyle();
  

  fHistOFF->Sumw2(); // if error has been set (via function SetBinError(j))
   // the errors set remain, i.e. are not overwritten with the sum of the square of weights.
  // Which means that this function will not have any effect.

  fHistOFF->SetXTitle("|alpha|  [\\circ]");
  fHistOFF->SetYTitle("Counts");
  fHistOFF->UseCurrentStyle();
  

  


/////////////////////////////////////

  fAlphamin = alphamin;
  fAlphamax = alphamax;
  fAlphammm = (alphamin+alphamax)/2.0;
  fAlphasig = alphasig;




  // UYpper limits for fit to OFF data set are taken also from alphamax
  // fAlphaminOFF is set in constructor. Inn principle, it WILL ALWAYS BE ZERO.

  fAlphamaxOFF = alphamax;

  
  fDraw     = drawpoly;
  fSavePlots = saveplots;
  fFitGauss = fitgauss;
  

  //--------------------------------------------
  // fit a polynomial in the backgr  region
  
  //*fLog << "MHFindSignificanceONOFF::FindSigma;  calling FitPolynomial()" << endl;
  if ( !FitPolynomialOFF())
  {
      *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; PolynomialOFF failed"
	    << endl;  
      return kFALSE;
  }
      
      
  //--------------------------------------------
  // calculate the number of excess events in the signal region
  
  //*fLog << "MHFindSignificanceONOFF::FindSigma;  calling DetExcess()" << endl;
  if ( !DetExcessONOFF())
  {
      *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; DetExcessONOFF failed"
	    << endl;  
      return kFALSE;
  }

  



  //--------------------------------------------
  // calculate the significance of the excess

  //*fLog << "MHFindSignificanceONOFF::FindSigma;  calling SigmaLiMa()" << endl;

  // For testing purposes "3 Li&Ma significances" will be computed
  // At some point, only one will remain


  Double_t siglima = 0.0;
  

  // Significance computed using effective number of OFF events  in signal 
  // region (fNoff) and gamma factor (fGama).
  // This is Wolfgang approach to the calulation of significance 
  // using Li&Ma formula and estimated OFF events from polynomial fit.
  
  if(fUseFittedQuantities) 
  {
      if ( !SigmaLiMa(fNon, fNoff, fGamma, &siglima) )
      {
	  *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; SigmaLiMa failed"
		<< endl;  
	  return kFALSE;
      }

      fSigLiMa = siglima;
  }
  else
  {
      fSigLiMa = 0.0;
  }



  // Significance computed using counted number of OFF events in signal 
  // region (fNoffSig) and normalization factor (fNormFactor).
  // This is the strictly speaking the significance in Li&Ma paper...
  
  
  if ( !SigmaLiMa(fNon, fNoffSig, fNormFactor, &siglima) )
  {
      *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; SigmaLiMa2 failed"
	    << endl;  
      return kFALSE;
  }
  fSigLiMa2 = siglima;
  
  
  // Significance computed using counted number of OFF events in signal 
  // region (fNoffSig) and normalization factor (fNormFactor). 
  // significance of gamma signal according to Li & Ma using  
  // formula (5)



  if ( !SigmaLiMaForm5(fNon, fNoffSig, fNormFactor, &siglima) )
      {
	  *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; SigmaLiMa failed"
		<< endl;  
	  return kFALSE;
      }

  fSigLiMa3 = siglima;
  
    
  
 
 
//--------------------------------------------
// calculate the error of the number of excess events
// using fitted quantities and counted quantities
 
 
// from fit to OFF histogram

 if (fSigLiMa > 0.0) 
 {fdNexONOFFFitted = fNexONOFFFitted / fSigLiMa;}
 else
 {fdNexONOFFFitted = 0.0;}


 // from counted OFF events
 if (fSigLiMa2 > 0.0) 
 {
     fdNexONOFF = fNexONOFF / fSigLiMa2;
 }
 else
 {
     fdNexONOFF = 0.0;
 }

  if (fDraw || fFitGauss)
  {
      
      // ------------------------------------------------
      // Polynomial fit to bkg region from ON data is performed, 
      // Because some of the quantities will be used as 
      // initial values in the PolyGAuss fit.
      
      // Besides, this function will modify the binning of fHist 
      // so that the number of events in each bin is big enough
      
      // This might change in future...
      
      if ( !FitPolynomial())
      {
	  *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; Polynomial fit failed"
		<< endl;  
	  return kFALSE;
      }
  }


  if (fDraw)
  {

      // Compute fHistOFFNormalized

      if (!ComputeHistOFFNormalized())
      {
	  *fLog << "MHFindSignificanceONOFF::ComputeHistOFFNormalized; Normalization of fHistOFF was not possible"
		<< endl;  
	  return kFALSE;
	  
      }
  }


 


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

  //*fLog << "MHFindSignificanceONOFF::FindSigma;  calling PrintPoly()" << endl;
  if (print)
    PrintPolyOFF();



  
  




  //--------------------------------------------
  // fit a (polynomial + Gauss) function to the ON-OFF alpha distribution 

  if (fFitGauss)
  {

    //--------------------------------------------------
    // delete objects from this fit
    // in order to have independent starting conditions for the next fit

      delete gMinuit;
      gMinuit = NULL;
    //--------------------------------------------------

       //*fLog << "MHFindSignificanceONOFF::FindSigma;  calling FitGaussPoly()" 
      //      << endl;
    if ( !FitGaussPoly() )
    {
      *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; FitGaussPoly failed"  
      	    << endl;  
      return kFALSE;
    }

    if (print)
    {
      //*fLog << "MHFindSignificanceONOFF::FindSigma;  calling PrintPolyGauss()" 
      //      << endl;
      PrintPolyGauss();
    }
  }

  //--------------------------------------------------
  // draw the histogram if requested

  if (fDraw)
  {
       // TEMPORALLY I will plot fHistOFF and fHistOFFNormalized (with the fits)

      
      if (!DrawHistOFF())
      {
	  *fLog << "MHFindSignificanceONOFF::DrawHistOFF; Drawing  of fHistOFF was not possible"
		<< endl;  
	  // return kFALSE;
	  
      }
      
     
      if (!DrawHistOFFNormalized())
      {
	  *fLog << "MHFindSignificanceONOFF::DrawHistOFFNormalized; Drawing  of fHistOFFNormalized was not possible"
		<< endl;  
	  // return kFALSE;
	  
      }




    //*fLog << "MHFindSignificanceONOFF::FindSigma;  calling DrawFit()" << endl;
      if ( !DrawFit() )
      {
	  *fLog << "MHFindSignificanceONOFF::FindSigmaONOFF; DrawFit failed"  
		<< endl;  
	  return kFALSE;
      }
  }
  

  //--------------------------------------------------
  // delete objects from this fit
  // in order to have independent starting conditions for the next fit

  delete gMinuit;
  gMinuit = NULL;
  //--------------------------------------------------
  
  return kTRUE;

}

 
//// ********************* end sigmaonoff ************************







// UNDER CONSTRUCTION

Bool_t MHFindSignificanceONOFF::SigmaVsAlphaONOFF(TH1 *fhistON, TH1 *fhistOFF,  
					   Double_t alphamin, Double_t alphamax, 
					   Int_t degree, Bool_t print)
{
    *fLog << "   MHFindSignificanceONOFF::SigmaVsAlphaONOFF still under construction !!!" << endl;

    return kFALSE;

}




// --------------------------------------------------------------------------
//
//  FitPolynomialOFF
//  - create a clone 'fHistOFF' of the |alpha| distribution 'fHistOrigOFF'
//  - fit a polynomial of degree 'fDegreeOFF' to the alpha distribution 
//    'fHistOFF' in the region alphaminOFF < |alpha| < alphamaxOFF



Bool_t MHFindSignificanceONOFF::FitPolynomialOFF()

{
  //--------------------------------------------------
  // check the histogram :
  //       - calculate initial values of the parameters
  //       - check for bins with zero entries
  //       - set minimum errors
  //       - save the original errors
  //       - set errors huge outside the fit range
  //         (in 'fcnpolyOFF' points with huge errors will be ignored)


  Double_t dummy = 1.e20;

  Double_t mean;
  Double_t rms;
  Double_t nclose;
  Double_t nfar;
  Double_t a2init = 0.0;
  TArrayD  saveError;

  Int_t nbins;
  Int_t nrebin = 1;


 //----------------   start while loop for rebinning   -----------------
  while(1)
  {

  fNzeroOFF   = 0;
  fMbinsOFF   = 0;
  fMlowOFF   = 0;
  fNoffTot = 0.0;

  // same variables (as in fitpoly to ON data) are used
  // here for naming the actual values for low/up limit for fit
  fAlphami =  10000.0;
  fAlphamm =  10000.0;
  fAlphama = -10000.0;

  mean   = 0.0;
  rms    = 0.0;
  nclose = 0.0;
  nfar   = 0.0;

  nbins = fHistOFF->GetNbinsX();
  saveError.Set(nbins);

  for (Int_t i=1; i<=nbins; i++)
  {
    saveError[i-1] = fHistOFF->GetBinError(i);

    // bin should be completely contained in the fit range
    // (fAlphaminOFF, fAlphamaxOFF)
    Double_t  xlo = fHistOFF->GetBinLowEdge(i);
    Double_t  xup = fHistOFF->GetBinLowEdge(i+1);

    if ( xlo >= fAlphaminOFF-fEps  &&  xlo <= fAlphamaxOFF+fEps  &&
	 xup >= fAlphaminOFF-fEps  &&  xup <= fAlphamaxOFF+fEps     )
    {
      fMbinsOFF++;

      if ( xlo < fAlphami )
        fAlphami = xlo;

      if ( xup > fAlphama )
        fAlphama = xup;

      Double_t content = fHistOFF->GetBinContent(i);
      // fNoffTot += content;

      mean += content;
      rms  += content*content;

// count events in low-alpha and high-alpha region
      if ( xlo >= fAlphammm-fEps  &&  xup >= fAlphammm-fEps)
      {
        nfar   += content;
        if ( xlo < fAlphamm )
          fAlphamm = xlo;
        if ( xup < fAlphamm )
          fAlphamm = xup; 
      }
      else
      {
        nclose += content;
        if ( xlo > fAlphamm )
          fAlphamm = xlo;
        if ( xup > fAlphamm )
          fAlphamm = xup; 
      }

      // count bins with zero entry
      if (content < 0.0001)
        fNzeroOFF++;
     
      // set minimum error
      // modified by K. Mase
      if (content < 0.0)
      {
        fMlowOFF += 1;
        fHistOFF->SetBinError(i, 3.0);
      }

      //*fLog << "Take : i, content, error = " << i << ",  " 
      //      << fHist->GetBinContent(i) << ",  "
      //      << fHist->GetBinError(i)   << endl;

      continue;
    }    
    // bin is not completely contained in the fit range : set error huge

    fHistOFF->SetBinError(i, dummy);

    //*fLog << "Omit : i, content, error = " << i << ",  " 
    //      << fHist->GetBinContent(i) << ",  " << fHist->GetBinError(i) 
    //      << endl;

  }

  // mean of entries/bin in the fit range
  if (fMbinsOFF > 0)
  {
    mean /= ((Double_t) fMbinsOFF);
    rms  /= ((Double_t) fMbinsOFF);
    rms = sqrt( rms - mean*mean );
  }

  

  // if there are no events in the background region
  //    there is no reason for rebinning
  //    and this is the condition for assuming a constant background (= 0)
  if (mean <= 0.0)
    break;

  Double_t helpmi = fAlphami*fAlphami*fAlphami;
  Double_t helpmm = fAlphamm*fAlphamm*fAlphamm;
  Double_t helpma = fAlphama*fAlphama*fAlphama;
  Double_t help   =   (helpma-helpmm) * (fAlphamm-fAlphami)    
	            - (helpmm-helpmi) * (fAlphama-fAlphamm);
  if (help != 0.0)
    a2init =  ( (fAlphamm-fAlphami)*nfar - (fAlphama-fAlphamm)*nclose )
                * 1.5 * fHistOFF->GetBinWidth(1) / help;
  else
    a2init = 0.0;


  //--------------------------------------------
  // rebin the histogram
  //   - if a bin has no entries 
  //   - or if there are too many bins with too few entries
  //   - or if the new bin width would exceed half the size of the 
  //     signal region

  if ( !fRebin  ||
       ( fNzeroOFF <= 0 && (Double_t)fMlowOFF<0.05*(Double_t)fMbinsOFF )  || 
       (Double_t)(nrebin+1)/(Double_t)nrebin * fHistOFF->GetBinWidth(1) 
                                                           > fAlphasig/2.0 )
  {
    //*fLog << "before break" << endl;
    break;
  }

  nrebin += 1;
  TString histname = fHistOFF->GetName();
  delete fHistOFF;
  fHistOFF = NULL;

  *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; rebin the |alpha|OFF plot, grouping "
        << nrebin << " bins together" << endl;

  // TH1::Rebin doesn't work properly
  //fHist = fHistOrig->Rebin(nrebin, "Rebinned");
  // use private routine RebinHistogram()
  fHistOFF = new TH1F;
  fHistOFF->Sumw2();
  fHistOFF->SetNameTitle(histname, histname);
  fHistOFF->UseCurrentStyle();

  // do rebinning such that x0 remains a lower bin edge
  Double_t x0 = 0.0;
  if ( !RebinHistogramOFF(x0, nrebin) )
  {
    *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; RebinHistgramOFF() failed" 
          << endl;
    return kFALSE;
  }

  fHistOFF->SetXTitle("|alpha|  [\\circ]");
  fHistOFF->SetYTitle("Counts");
  
  }
  //----------------   end of while loop for rebinning   -----------------


  // if there are still too many bins with too few entries don't fit
  // and assume a constant background

  fConstantBackg = kFALSE;
  if ( fNzeroOFF > 0  ||  (Double_t)fMlowOFF>0.05*(Double_t)fMbinsOFF ) 
  {
    *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; polynomial fit not possible,  fNzeroOFF, fMlowOFF, fMbinsOFF = "
          << fNzeroOFF << ",  " << fMlowOFF << ",  " << fMbinsOFF << endl;
    *fLog << "                    assume a constant background" << endl;

    fConstantBackg = kTRUE;
    fDegreeOFF        = 0;

    TString funcname = "PolyOFF";
    Double_t xmin =   0.0;
    Double_t xmax =  90.0;

    TString formula = "[0]";

    fPolyOFF = new TF1(funcname, formula, xmin, xmax);
    TList *funclist = fHistOFF->GetListOfFunctions();
    funclist->Add(fPolyOFF);

    //--------------------
    Int_t nparfree = 1;
    fChisqOFF         = 0.0; 
    fNdfOFF           = fMbinsOFF - nparfree;
    fProbOFF          = 0.0;
    fIstatOFF         = 0;

    fValuesOFF.Set(1);
    fErrorsOFF.Set(1);

    Double_t val, err;
    val = mean;
    err = rms;
    

    fPolyOFF->SetParameter(0, val);
    fPolyOFF->SetParError (0, err);

    fValuesOFF[0] = val;
    fErrorsOFF[0] = err; 

    fEmaOFF[0][0]  = err*err;
    fCorrOFF[0][0] = 1.0;
    //--------------------
    //--------------------------------------------------
    // reset the errors of the points in the histogram
    for (Int_t i=1; i<=nbins; i++)
    {
      fHistOFF->SetBinError(i, saveError[i-1]);
    }


    return kTRUE;
  }


  //===========   start loop for reducing the degree   ==================
  //              of the polynomial
  while (1)
  {
      //--------------------------------------------------
      // prepare fit of a polynomial :   (a0 + a1*x + a2*x**2 + a3*x**3 + ...)


      TString funcname = "PolyOFF";
      Double_t xmin =   0.0;
      Double_t xmax =  90.0;

      TString formula = "[0]";
      TString bra1     = "+[";
      TString bra2     =    "]";
      TString xpower   = "*x";
      TString newpower = "*x";
      for (Int_t i=1; i<=fDegreeOFF; i++)
      {
          formula += bra1;
          formula += i;
          formula += bra2;
          formula += xpower;

          xpower += newpower;
      }

      //*fLog << "FitPolynomial : formula = " << formula << endl;

      fPolyOFF = new TF1(funcname, formula, xmin, xmax);
      TList *funclist = fHistOFF->GetListOfFunctions();
      funclist->Add(fPolyOFF);

      //------------------------
      // attention : the dimensions must agree with those in CallMinuit()
      const UInt_t npar = fDegreeOFF+1;

      TString parname[npar];
      TArrayD vinit(npar);
      TArrayD  step(npar);
      TArrayD limlo(npar);
      TArrayD limup(npar);
      TArrayI   fix(npar);

      vinit[0] =   mean;
      vinit[2] = a2init;

      for (UInt_t j=0; j<npar; j++)
      {
          parname[j]  = "p";
          parname[j] += j+1;

          step[j] = vinit[j] != 0.0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
      }

      // limit the first coefficient of the polynomial to positive values
      // because the background must not be negative
      limup[0] = fHistOFF->GetEntries();

      // use the subsequernt loop if you want to apply the
      // constraint : uneven derivatives (at alpha=0) = zero
      for (UInt_t j=1; j<npar; j+=2)
      {
          vinit[j] = 0;
          step[j]  = 0;
          fix[j]   = 1;
      }

      //*fLog << "FitPolynomial : before CallMinuit()" << endl;

      MMinuitInterface inter;
      const Bool_t rc = inter.CallMinuit(fcnpolyOFF, parname, vinit, step,
                                         limlo, limup, fix, fHistOFF, "Migrad",
                                         kFALSE);

      //*fLog << "FitPolynomial : after CallMinuit()" << endl;

      if (rc != 0)
      {
          //  *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit failed"
          //        << endl;
          //  return kFALSE;
      }

      //-------------------
      // get status of minimization
      Double_t fmin   = 0;
      Double_t fedm   = 0;
      Double_t errdef = 0;
      Int_t    npari  = 0;
      Int_t    nparx  = 0;

      if (gMinuit)
          gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fIstatOFF);

      *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; fmin, fedm, errdef, npari, nparx, fIstat = "
            << fmin << ",  " << fedm << ",  " << errdef << ",  " << npari
            << ",  " << nparx << ",  " << fIstatOFF << endl;


      //-------------------
      // store the results

      Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
      fChisqOFF         = fmin;
      fNdfOFF          = fMbinsOFF - nparfree;
      fProbOFF          = TMath::Prob(fChisqOFF, fNdfOFF);


      // get fitted parameter values and errors
      fValuesOFF.Set(npar);
      fErrorsOFF.Set(npar);

      for (Int_t j=0; j<=fDegreeOFF; j++)
      {
          Double_t val, err;
          if (gMinuit)
              gMinuit->GetParameter(j, val, err);

          fPolyOFF->SetParameter(j, val);
          fPolyOFF->SetParError(j, err);

          fValuesOFF[j] = val;
          fErrorsOFF[j] = err;
      }

      // if the highest coefficient (j0) of the polynomial
      // is consistent with zero reduce the degree of the polynomial

      Int_t j0 = 0;
      for (Int_t j=fDegreeOFF; j>1; j--)
      {
          // ignore fixed parameters
          if (fErrorsOFF[j] == 0)
              continue;

          // this is the highest coefficient
          j0 = j;
          break;
      }

      if (!fReduceDegree || j0==0 || TMath::Abs(fValuesOFF[j0]) > fErrorsOFF[j0])
          break;

      // reduce the degree of the polynomial
      *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; reduce the degree of the polynomial from "
          << fDegreeOFF << " to " << (j0-2) << endl;
      fDegreeOFF = j0 - 2;

      funclist->Remove(fPolyOFF);
      //if (fPoly)
      delete fPolyOFF;
      fPolyOFF = NULL;

      // delete the Minuit object in order to have independent starting
      // conditions for the next minimization
      //if (gMinuit)
      delete gMinuit;
      gMinuit = NULL;
  }
  //===========   end of loop for reducing the degree   ==================
  //              of the polynomial


  //--------------------------------------------------
  // get the error matrix of the fitted parameters

  if (fIstatOFF >= 1)
  {
      // error matrix was calculated
      if (gMinuit)
          gMinuit->mnemat(&fEmatOFF[0][0], fNdimOFF);

      // copy covariance matrix into a matrix which includes also the fixed
      // parameters
      TString  name;
      Double_t bnd1, bnd2, val, err;
      Int_t    jvarbl;
      Int_t    kvarbl;
      for (Int_t j=0; j<=fDegreeOFF; j++)
      {
          if (gMinuit)
              gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);

          for (Int_t k=0; k<=fDegreeOFF; k++)
          {
              if (gMinuit)
                  gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);

              fEmaOFF[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fEmatOFF[jvarbl-1][kvarbl-1];
          }
      }
  }
  else
  {
      // error matrix was not calculated, construct it
      *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; error matrix not defined"
          << endl;
      for (Int_t j=0; j<=fDegreeOFF; j++)
      {
          for (Int_t k=0; k<=fDegreeOFF; k++)
              fEmaOFF[j][k] = 0;

          fEmaOFF[j][j] = fErrorsOFF[j]*fErrorsOFF[j];
      }
  }



  
  //--------------------------------------------------
  // calculate correlation matrix
  for (Int_t j=0; j<=fDegreeOFF; j++)
  {
      for (Int_t k=0; k<=fDegreeOFF; k++)
      {
          const Double_t sq = fEmaOFF[j][j]*fEmaOFF[k][k];
          fCorrOFF[j][k] = 
	      sq == 0 ? 0 : (fEmaOFF[j][k] / TMath::Sqrt(fEmaOFF[j][j]*fEmaOFF[k][k]));
      }
  }
  
  //--------------------------------------------------
  // reset the errors of the points in the histogram
  for (Int_t i=1; i<=nbins; i++)
  {
      fHistOFF->SetBinError(i, saveError[i-1]);
  }
  
  return kTRUE;


}








// --------------------------------------------------------------------------
//
//  FitPolynomial
//
//  - create a clone 'fHist' of the |alpha| distribution 'fHistOrig'
//  - fit a polynomial of degree 'fDegree' to the alpha distribution 
//    'fHist' in the region alphamin < |alpha| < alphamax
//
//  in pathological cases the histogram is rebinned before fitting
//     (this is done only if fRebin is kTRUE)
//
//  if the highest coefficient of the polynomial is compatible with zero
//     the fit is repeated with a polynomial of lower degree
//     (this is done only if fReduceDegree is kTRUE)
//
//

Bool_t MHFindSignificanceONOFF::FitPolynomial()
{
  //--------------------------------------------------
  // check the histogram :
  //       - calculate initial values of the parameters
  //       - check for bins with zero entries
  //       - set minimum errors
  //       - save the original errors
  //       - set errors huge outside the fit range
  //         (in 'fcnpoly' points with huge errors will be ignored)


  Double_t dummy = 1.e20;

  Double_t mean;
  Double_t rms;
  Double_t nclose;
  Double_t nfar;
  Double_t a2init = 0.0;
  TArrayD  saveError;

  Int_t nbins;
  Int_t nrebin = 1;

  //----------------   start while loop for rebinning   -----------------
  while(1)
  {

  fNzero   = 0;
  fMbins   = 0;
  fMlow    = 0;
  fNbgtot  = 0.0;

  fAlphami =  10000.0;
  fAlphamm =  10000.0;
  fAlphama = -10000.0;

  mean   = 0.0;
  rms    = 0.0;
  nclose = 0.0;
  nfar   = 0.0;

  nbins = fHist->GetNbinsX();
  saveError.Set(nbins);

  for (Int_t i=1; i<=nbins; i++)
  {
    saveError[i-1] = fHist->GetBinError(i);

    // bin should be completely contained in the fit range
    // (fAlphamin, fAlphamax)
    Double_t  xlo = fHist->GetBinLowEdge(i);
    Double_t  xup = fHist->GetBinLowEdge(i+1);

    if ( xlo >= fAlphamin-fEps  &&  xlo <= fAlphamax+fEps  &&
	 xup >= fAlphamin-fEps  &&  xup <= fAlphamax+fEps     )
    {
      fMbins++;

      if ( xlo < fAlphami )
        fAlphami = xlo;

      if ( xup > fAlphama )
        fAlphama = xup;

      Double_t content = fHist->GetBinContent(i);
      fNbgtot += content;

      mean += content;
      rms  += content*content;

      // count events in low-alpha and high-alpha region
      if ( xlo >= fAlphammm-fEps  &&  xup >= fAlphammm-fEps)
      {
        nfar   += content;
        if ( xlo < fAlphamm )
          fAlphamm = xlo;
        if ( xup < fAlphamm )
          fAlphamm = xup; 
      }
      else
      {
        nclose += content;
        if ( xlo > fAlphamm )
          fAlphamm = xlo;
        if ( xup > fAlphamm )
          fAlphamm = xup; 
      }

      // count bins with zero entry
      /*if (content <= 0.0)
        fNzero++;*/
     
      // set minimum error
      /*if (content < 9.0)
      {
        fMlow += 1;
        fHist->SetBinError(i, 3.0);
	}*/

      //*fLog << "Take : i, content, error = " << i << ",  " 
      //      << fHist->GetBinContent(i) << ",  "
      //      << fHist->GetBinError(i)   << endl;

      continue;
    }    
    // bin is not completely contained in the fit range : set error huge

    fHist->SetBinError(i, dummy);

    //*fLog << "Omit : i, content, error = " << i << ",  " 
    //      << fHist->GetBinContent(i) << ",  " << fHist->GetBinError(i) 
    //      << endl;

  }

  // mean of entries/bin in the fit range
  if (fMbins > 0)
  {
    mean /= ((Double_t) fMbins);
    rms  /= ((Double_t) fMbins);
  }

  rms = sqrt( rms - mean*mean );

  // if there are no events in the background region
  //    there is no reason for rebinning
  //    and this is the condition for assuming a constant background (= 0)
  if (mean <= 0.0)
    break;

  Double_t helpmi = fAlphami*fAlphami*fAlphami;
  Double_t helpmm = fAlphamm*fAlphamm*fAlphamm;
  Double_t helpma = fAlphama*fAlphama*fAlphama;
  Double_t help   =   (helpma-helpmm) * (fAlphamm-fAlphami)    
	            - (helpmm-helpmi) * (fAlphama-fAlphamm);
  if (help != 0.0)
    a2init =  ( (fAlphamm-fAlphami)*nfar - (fAlphama-fAlphamm)*nclose )
                * 1.5 * fHist->GetBinWidth(1) / help;
  else
    a2init = 0.0;


  //--------------------------------------------
  // rebin the histogram
  //   - if a bin has no entries 
  //   - or if there are too many bins with too few entries
  //   - or if the new bin width would exceed half the size of the 
  //     signal region

  if ( !fRebin  ||
       ( fNzero <= 0 && (Double_t)fMlow<0.05*(Double_t)fMbins )  || 
       (Double_t)(nrebin+1)/(Double_t)nrebin * fHist->GetBinWidth(1) 
                                                           > fAlphasig/2.0 )
  {
    //*fLog << "before break" << endl;
    break;
  }

  nrebin += 1;
  TString histname = fHist->GetName();
  delete fHist;
  fHist = NULL;

  *fLog << "MHFindSignificanceONOFF::FitPolynomial; rebin the |alpha| plot, grouping "
        << nrebin << " bins together" << endl;

  // TH1::Rebin doesn't work properly
  //fHist = fHistOrig->Rebin(nrebin, "Rebinned");
  // use private routine RebinHistogram()
  fHist = new TH1F;
  fHist->Sumw2();
  fHist->SetNameTitle(histname, histname);
  fHist->UseCurrentStyle();

  // do rebinning such that x0 remains a lower bin edge
  Double_t x0 = 0.0;
  if ( !RebinHistogram(x0, nrebin) )
  {
    *fLog << "MHFindSignificanceONOFF::FitPolynomial; RebinHistgram() failed" 
          << endl;
    return kFALSE;
  }

  fHist->SetXTitle("|alpha|  [\\circ]");
  fHist->SetYTitle("Counts");

  }
  //----------------   end of while loop for rebinning   -----------------


  // if there are still too many bins with too few entries don't fit
  // and assume a constant background

  fConstantBackg = kFALSE;
  if ( fNzero > 0  ||  (Double_t)fMlow>0.05*(Double_t)fMbins ) 
  {
    *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit not possible,  fNzero, fMlow, fMbins = "
          << fNzero << ",  " << fMlow << ",  " << fMbins << endl;
    *fLog << "                    assume a constant background" << endl;

    fConstantBackg = kTRUE;
    fDegree        = 0;

    TString funcname = "Poly";
    Double_t xmin =   0.0;
    Double_t xmax =  90.0;

    TString formula = "[0]";

    fPoly = new TF1(funcname, formula, xmin, xmax);
    TList *funclist = fHist->GetListOfFunctions();
    funclist->Add(fPoly);

    //--------------------
    Int_t nparfree = 1;
    fChisq         = 0.0; 
    fNdf           = fMbins - nparfree;
    fProb          = 0.0;
    fIstat         = 0;

    fValues.Set(1);
    fErrors.Set(1);

    Double_t val, err;
    val = mean;
    err = sqrt( mean / (Double_t)fMbins );

    fPoly->SetParameter(0, val);
    fPoly->SetParError (0, err);

    fValues[0] = val;
    fErrors[0] = err; 

    fEma[0][0]  = err*err;
    fCorr[0][0] = 1.0;
    //--------------------

    //--------------------------------------------------
    // reset the errors of the points in the histogram
    for (Int_t i=1; i<=nbins; i++)
    {
      fHist->SetBinError(i, saveError[i-1]);
    }


    return kTRUE;
  }


  //===========   start loop for reducing the degree   ==================
  //              of the polynomial
  while (1)
  {
      //--------------------------------------------------
      // prepare fit of a polynomial :   (a0 + a1*x + a2*x**2 + a3*x**3 + ...)
 
      TString funcname = "Poly";
      Double_t xmin =   0.0;
      Double_t xmax =  90.0;

      TString formula = "[0]";
      TString bra1     = "+[";
      TString bra2     =    "]";
      TString xpower   = "*x";
      TString newpower = "*x";
      for (Int_t i=1; i<=fDegree; i++)
      {
          formula += bra1;
          formula += i;
          formula += bra2;
          formula += xpower;

          xpower += newpower;
      }

      //*fLog << "FitPolynomial : formula = " << formula << endl;

      fPoly = new TF1(funcname, formula, xmin, xmax);
      TList *funclist = fHist->GetListOfFunctions();
      funclist->Add(fPoly);

      //------------------------
      // attention : the dimensions must agree with those in CallMinuit()
      const UInt_t npar = fDegree+1;

      TString parname[npar];
      TArrayD vinit(npar);
      TArrayD  step(npar);
      TArrayD limlo(npar);
      TArrayD limup(npar);
      TArrayI   fix(npar);

      vinit[0] =   mean;
      vinit[2] = a2init;

      for (UInt_t j=0; j<npar; j++)
      {
          parname[j]  = "p";
          parname[j] += j+1;

          step[j] = vinit[j] != 0.0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
      }

      // limit the first coefficient of the polynomial to positive values
      // because the background must not be negative
      limup[0] = fHist->GetEntries();

      // use the subsequernt loop if you want to apply the
      // constraint : uneven derivatives (at alpha=0) = zero
      for (UInt_t j=1; j<npar; j+=2)
      {
          vinit[j] = 0;
          step[j]  = 0;
          fix[j]   = 1;
      }

      //*fLog << "FitPolynomial : before CallMinuit()" << endl;

      MMinuitInterface inter;
      const Bool_t rc = inter.CallMinuit(fcnpoly, parname, vinit, step,
                                         limlo, limup, fix, fHist, "Migrad",
                                         kFALSE);

      //*fLog << "FitPolynomial : after CallMinuit()" << endl;

      if (rc != 0)
      {
          //  *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit failed"
          //        << endl;
          //  return kFALSE;
      }


      //-------------------
      // get status of minimization
      Double_t fmin   = 0;
      Double_t fedm   = 0;
      Double_t errdef = 0;
      Int_t    npari  = 0;
      Int_t    nparx  = 0;

      if (gMinuit)
          gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fIstat);

      *fLog << "MHFindSignificanceONOFF::FitPolynomial; fmin, fedm, errdef, npari, nparx, fIstat = "
            << fmin << ",  " << fedm << ",  " << errdef << ",  " << npari
            << ",  " << nparx << ",  " << fIstat << endl;


      //-------------------
      // store the results

      Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
      fChisq         = fmin;
      fNdf           = fMbins - nparfree;
      fProb          = TMath::Prob(fChisq, fNdf);


      // get fitted parameter values and errors
      fValues.Set(npar);
      fErrors.Set(npar);

      for (Int_t j=0; j<=fDegree; j++)
      {
          Double_t val, err;
          if (gMinuit)
              gMinuit->GetParameter(j, val, err);

          fPoly->SetParameter(j, val);
          fPoly->SetParError(j, err);

          fValues[j] = val;
          fErrors[j] = err;
      }


      //--------------------------------------------------
      // if the highest coefficient (j0) of the polynomial
      // is consistent with zero reduce the degree of the polynomial

      Int_t j0 = 0;
      for (Int_t j=fDegree; j>1; j--)
      {
          // ignore fixed parameters
          if (fErrors[j] == 0)
              continue;

          // this is the highest coefficient
          j0 = j;
          break;
      }

      if (!fReduceDegree || j0==0 || TMath::Abs(fValues[j0]) > fErrors[j0])
          break;

      // reduce the degree of the polynomial
      *fLog << "MHFindSignificanceONOFF::FitPolynomial; reduce the degree of the polynomial from "
          << fDegree << " to " << (j0-2) << endl;
      fDegree = j0 - 2;

      funclist->Remove(fPoly);
      //if (fPoly)
      delete fPoly;
      fPoly = NULL;

      // delete the Minuit object in order to have independent starting
      // conditions for the next minimization
      //if (gMinuit)
      delete gMinuit;
      gMinuit = NULL;
  }
  //===========   end of loop for reducing the degree   ==================
  //              of the polynomial


  //--------------------------------------------------
  // get the error matrix of the fitted parameters


  if (fIstat >= 1)
  {
      // error matrix was calculated
      if (gMinuit)
          gMinuit->mnemat(&fEmat[0][0], fNdim);

      // copy covariance matrix into a matrix which includes also the fixed
      // parameters
      TString  name;
      Double_t bnd1, bnd2, val, err;
      Int_t    jvarbl;
      Int_t    kvarbl;
      for (Int_t j=0; j<=fDegree; j++)
      {
          if (gMinuit)
              gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);

          for (Int_t k=0; k<=fDegree; k++)
          {
              if (gMinuit)
                  gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);

              fEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fEmat[jvarbl-1][kvarbl-1];
          }
      }
  }
  else
  {
      // error matrix was not calculated, construct it
      *fLog << "MHFindSignificanceONOFF::FitPolynomial; error matrix not defined"
          << endl;
      for (Int_t j=0; j<=fDegree; j++)
      {
          for (Int_t k=0; k<=fDegree; k++)
              fEma[j][k] = 0;

          fEma[j][j] = fErrors[j]*fErrors[j];
      }
  }


  //--------------------------------------------------
  // calculate correlation matrix
  for (Int_t j=0; j<=fDegree; j++)
      for (Int_t k=0; k<=fDegree; k++)
      {
          const Double_t sq = fEma[j][j]*fEma[k][k];
          fCorr[j][k] = sq==0 ? 0 : fEma[j][k] / TMath::Sqrt(fEma[j][j]*fEma[k][k]);
      }


  //--------------------------------------------------
  // reset the errors of the points in the histogram
  for (Int_t i=1; i<=nbins; i++)
      fHist->SetBinError(i, saveError[i-1]);


  return kTRUE;
}

// --------------------------------------------------------------------------
//
// ReBinHistogramOFF
//
// rebin the histogram 'fHistOrig' by grouping 'nrebin' bins together 
// put the result into the histogram 'fHist'
// the rebinning is made such that 'x0' remains a lower bound of a bin
//

Bool_t MHFindSignificanceONOFF::RebinHistogramOFF(Double_t x0, Int_t nrebin)
{
  //-----------------------------------------
  // search bin i0 which has x0 as lower edge

  Int_t i0 = -1;
  Int_t nbold = fHistOrigOFF->GetNbinsX();
  for (Int_t i=1; i<=nbold; i++)
  {
      if (TMath::Abs(fHistOrigOFF->GetBinLowEdge(i) - x0) < 1.e-4 )
      {
          i0 = i;
          break;
      }
  }

  if (i0 == -1)
  {
    i0 = 1;
    *fLog << "MHFindsignificanceOFF::RebinOFF; no bin found with " << x0
          << " as lower edge,  start rebinning with bin 1" << endl;
  }

  Int_t istart = i0 - nrebin * ( (i0-1)/nrebin );

  //-----------------------------------------
  // get new bin edges

  const Int_t    nbnew = (nbold-istart+1) / nrebin;
  const Double_t xmin  = fHistOrigOFF->GetBinLowEdge(istart);
  const Double_t xmax  = xmin + (Double_t)nbnew * nrebin * fHistOrigOFF->GetBinWidth(1);
  fHistOFF->SetBins(nbnew, xmin, xmax);

  *fLog << "MHFindSignificanceONOFF::ReBinOFF; x0, i0, nbold, nbnew, xmin, xmax = "
        << x0 << ",  " << i0 << ",  " << nbold << ",  " << nbnew << ",  "
        << xmin << ",  " << xmax << endl;

  //-----------------------------------------
  // get new bin entries

  for (Int_t i=1; i<=nbnew; i++)
  {
      Int_t j = nrebin*(i-1) + istart;

      Double_t content = 0;
      Double_t error2  = 0;
      for (Int_t k=0; k<nrebin; k++)
      {
          content += fHistOrigOFF->GetBinContent(j+k);
          error2  += fHistOrigOFF->GetBinError(j+k) * fHistOrigOFF->GetBinError(j+k);
      }
      fHistOFF->SetBinContent(i, content);
      fHistOFF->SetBinError  (i, sqrt(error2));
  }
  fHistOFF->SetEntries( fHistOrigOFF->GetEntries() );

  return kTRUE;
}




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



//
// ReBinHistogram
//
// rebin the histogram 'fHistOrig' by grouping 'nrebin' bins together 
// put the result into the histogram 'fHist'
// the rebinning is made such that 'x0' remains a lower bound of a bin
//

Bool_t MHFindSignificanceONOFF::RebinHistogram(Double_t x0, Int_t nrebin)
{
  //-----------------------------------------
  // search bin i0 which has x0 as lower edge

  Int_t i0 = -1;
  Int_t nbold = fHistOrig->GetNbinsX();
  for (Int_t i=1; i<=nbold; i++)
  {
      if (TMath::Abs(fHistOrig->GetBinLowEdge(i) - x0) < 1.e-4 )
      {
          i0 = i;
          break;
      }
  }

  if (i0 == -1)
  {
    i0 = 1;
    *fLog << "MHFindsignificance::Rebin; no bin found with " << x0
          << " as lower edge,  start rebinning with bin 1" << endl;
  }

  Int_t istart = i0 - nrebin * ( (i0-1)/nrebin );

  //-----------------------------------------
  // get new bin edges

  const Int_t    nbnew = (nbold-istart+1) / nrebin;
  const Double_t xmin  = fHistOrig->GetBinLowEdge(istart);
  const Double_t xmax  = xmin + (Double_t)nbnew * nrebin * fHistOrig->GetBinWidth(1);
  fHist->SetBins(nbnew, xmin, xmax);

  *fLog << "MHFindSignificanceONOFF::ReBin; x0, i0, nbold, nbnew, xmin, xmax = "
        << x0 << ",  " << i0 << ",  " << nbold << ",  " << nbnew << ",  "
        << xmin << ",  " << xmax << endl;

  //-----------------------------------------
  // get new bin entries

  for (Int_t i=1; i<=nbnew; i++)
  {
      Int_t j = nrebin*(i-1) + istart;

      Double_t content = 0;
      Double_t error2  = 0;
      for (Int_t k=0; k<nrebin; k++)
      {
          content += fHistOrig->GetBinContent(j+k);
          error2  += fHistOrig->GetBinError(j+k) * fHistOrig->GetBinError(j+k);
      }
      fHist->SetBinContent(i, content);
      fHist->SetBinError  (i, sqrt(error2));
  }
  fHist->SetEntries( fHistOrig->GetEntries() );

  return kTRUE;
}

// --------------------------------------------------------------------------
//
//  FitGaussPoly
//
//  fits a (Gauss + polynomial function) to the alpha distribution 'fhist' 
//
//
Bool_t MHFindSignificanceONOFF::FitGaussPoly()
{
  *fLog << "Entry FitGaussPoly" << endl;

  //--------------------------------------------------
  // check the histogram :
  //       - calculate initial values of the parameters
  //       - check for bins with zero entries
  //       - set minimum errors
  //       - save the original errors
  //       - set errors huge outside the fit range
  //         (in 'fcnpoly' points with huge errors will be ignored)


  Double_t dummy = 1.e20;

  fGNzero   = 0;
  fGMbins   = 0;

  //------------------------------------------
  // if a constant background has been assumed (due to low statistics)
  // fit only in the signal region
  if ( !fConstantBackg )
  {
    fAlphalow = 0.0;
    fAlphahig = fAlphamax;
  }
  else
  {
    fAlphalow = 0.0;
    fAlphahig = 2.0*fAlphasig>25.0 ? 25.0 : 2.0*fAlphasig;
  }
  //------------------------------------------


  fAlphalo =  10000.0;
  fAlphahi = -10000.0;


  Int_t nbins = fHist->GetNbinsX();
  TArrayD saveError(nbins);

  for (Int_t i=1; i<=nbins; i++)
  {
    saveError[i-1] = fHist->GetBinError(i);

    // bin should be completely contained in the fit range
    // (fAlphalow, fAlphahig)
    Double_t  xlo = fHist->GetBinLowEdge(i);
    Double_t  xup = fHist->GetBinLowEdge(i+1);

    if ( xlo >= fAlphalow-fEps  &&  xlo <= fAlphahig+fEps  &&
	 xup >= fAlphalow-fEps  &&  xup <= fAlphahig+fEps     )
    {
      fGMbins++;

      if ( xlo < fAlphalo )
        fAlphalo = xlo;

      if ( xup > fAlphahi )
        fAlphahi = xup;

      Double_t content = fHist->GetBinContent(i);


      // count bins with zero entry
      if (content <= 0.0)
        fGNzero++;
     
      // set minimum error
      if (content < 9.0)
        fHist->SetBinError(i, 3.0);

      //*fLog << "Take : i, content, error = " << i << ",  " 
      //      << fHist->GetBinContent(i) << ",  "
      //      << fHist->GetBinError(i)   << endl;

      continue;
    }    
    // bin is not completely contained in the fit range : set error huge

    fHist->SetBinError(i, dummy);

    //*fLog << "Omit : i, content, error = " << i << ",  " 
    //      << fHist->GetBinContent(i) << ",  " << fHist->GetBinError(i) 
    //      << endl;

  }


  // if a bin has no entries don't fit
  if (fGNzero > 0)
  {
    *fLog << "MHFindSignificanceONOFF::FitGaussPoly; out of " << fGMbins 
          << " bins there are " << fGNzero
          << " bins with zero entry" << endl;

    fGPoly = NULL;
    return kFALSE;
  }


  //--------------------------------------------------
  // prepare fit of a (polynomial+Gauss) :   
  // (a0 + a1*x + a2*x**2 + a3*x**3 + ...) + A*exp( -0.5*((x-x0)/sigma)**2 )

  TString funcname = "PolyGauss";
  Double_t xmin =   0.0;
  Double_t xmax =  90.0;

  TString xpower   = "*x";
  TString newpower = "*x";

  TString formulaBackg = "[0]";
  for (Int_t i=1; i<=fDegree; i++)
      formulaBackg += Form("+[%d]*x^%d", i, i);

  const TString formulaGauss = 
        Form("[%d]/[%d]*exp(-0.5*((x-[%d])/[%d])^2)",
             fDegree+1, fDegree+3, fDegree+2, fDegree+3);

  TString formula = formulaBackg;
  formula += "+";
  formula += formulaGauss;

  *fLog << "FitGaussPoly : formulaBackg = " << formulaBackg << endl;
  *fLog << "FitGaussPoly : formulaGauss = " << formulaGauss << endl;
  *fLog << "FitGaussPoly : formula = " << formula << endl;

  fGPoly = new TF1(funcname, formula, xmin, xmax);
  TList *funclist = fHist->GetListOfFunctions();
  funclist->Add(fGPoly);

  fGBackg = new TF1("Backg", formulaBackg, xmin, xmax);
  //funclist->Add(fGBackg); 

  //------------------------
  // attention : the dimensions must agree with those in CallMinuit()
  Int_t npar = fDegree+1 + 3;

  TString parname[npar];
  TArrayD vinit(npar);
  TArrayD  step(npar); 
  TArrayD limlo(npar); 
  TArrayD limup(npar); 
  TArrayI   fix(npar);


  // take as initial values for the polynomial 
  // the result from the polynomial fit
  for (Int_t j=0; j<=fDegree; j++)
    vinit[j] = fPoly->GetParameter(j);

  Double_t sigma = 8;
  vinit[fDegree+1] = 2.0 * fNexONOFF * fHist->GetBinWidth(1) / TMath::Sqrt(TMath::Pi()*2);
  vinit[fDegree+2] = 0;
  vinit[fDegree+3] = sigma;

  *fLog << "FitGaussPoly : starting value for Gauss-amplitude = " 
        << vinit[fDegree+1] << endl;

  for (Int_t j=0; j<npar; j++)
  {
      parname[j]  = "p";
      parname[j] += j+1;

      step[j] = vinit[j]!=0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
  }

  // limit the first coefficient of the polynomial to positive values
  // because the background must not be negative
  limup[0] = fHist->GetEntries()*10;

  // limit the sigma of the Gauss function
  limup[fDegree+3] = 20;


  // use the subsequernt loop if you want to apply the
  // constraint : uneven derivatives (at alpha=0) = zero
  for (Int_t j=1; j<=fDegree; j+=2)
  {
      vinit[j] = 0;
      step[j]  = 0;
      fix[j]   = 1;
  }

  // fix position of Gauss function
  vinit[fDegree+2] = 0;
  step[fDegree+2]  = 0;
  fix[fDegree+2]   = 1;
   
  // if a constant background has been assumed (due to low statistics)
  // fix the background
  if (fConstantBackg)
  {
      step[0] = 0;
      fix[0]  = 1;
  }

  MMinuitInterface inter;
  const Bool_t rc = inter.CallMinuit(fcnpolygauss, parname, vinit, step,
                                     limlo, limup, fix, fHist, "Migrad",
                                     kFALSE);

  if (rc != 0)
  {
  //  *fLog << "MHFindSignificanceONOFF::FitGaussPoly; (polynomial+Gauss) fit failed"
  //        << endl;
  //  return kFALSE;
  }


  //-------------------
  // get status of the minimization
  Double_t fmin;
  Double_t fedm;
  Double_t errdef;
  Int_t    npari;
  Int_t    nparx;

  if (gMinuit)
    gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, fGIstat);

  *fLog << "MHFindSignificanceONOFF::FitGaussPoly; fmin, fedm, errdef, npari, nparx, fGIstat = "
        << fmin << ",  " << fedm << ",  " << errdef << ",  " << npari
        << ",  " << nparx << ",  " << fGIstat << endl;


  //-------------------
  // store the results

  Int_t nparfree  = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0;
  fGChisq         = fmin; 
  fGNdf           = fGMbins - nparfree;
  fGProb          = TMath::Prob(fGChisq, fGNdf);


  // get fitted parameter values and errors
  fGValues.Set(npar);
  fGErrors.Set(npar);

  for (Int_t j=0; j<npar; j++)
  {
    Double_t val, err;
    if (gMinuit)
      gMinuit->GetParameter(j, val, err);

    fGPoly->SetParameter(j, val);
    fGPoly->SetParError(j, err);

    fGValues[j] = val;
    fGErrors[j] = err; 

    if (j <=fDegree)
    {
      fGBackg->SetParameter(j, val);
      fGBackg->SetParError(j, err);
    }
  }

  fSigmaGauss  = fGValues[fDegree+3];
  fdSigmaGauss = fGErrors[fDegree+3];
  // fitted total number of excess events 
  fNexGauss = fGValues[fDegree+1] * TMath::Sqrt(TMath::Pi()*2) /
                                         (fHist->GetBinWidth(1)*2 );
  fdNexGauss = fNexGauss * fGErrors[fDegree+1]/fGValues[fDegree+1];

  //--------------------------------------------------
  // get the error matrix of the fitted parameters


  if (fGIstat >= 1)
  {
    // error matrix was calculated
    if (gMinuit)
      gMinuit->mnemat(&fGEmat[0][0], fGNdim);

    // copy covariance matrix into a matrix which includes also the fixed
    // parameters
    TString  name;
    Double_t bnd1, bnd2, val, err;
    Int_t    jvarbl;
    Int_t    kvarbl;
    for (Int_t j=0; j<npar; j++)
    {
        if (gMinuit)
            gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);

        for (Int_t k=0; k<npar; k++)
        {
            if (gMinuit)
                gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);

            fGEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fGEmat[jvarbl-1][kvarbl-1];
        }
    }
  }
  else
  {
    // error matrix was not calculated, construct it
    *fLog << "MHFindSignificanceONOFF::FitPolynomial; error matrix not defined" 
          << endl;
    for (Int_t j=0; j<npar; j++)
    {
        for (Int_t k=0; k<npar; k++)
            fGEma[j][k] = 0;

        fGEma[j][j] = fGErrors[j]*fGErrors[j];
    }
  }


  //--------------------------------------------------
  // calculate correlation matrix
  for (Int_t j=0; j<npar; j++)
  {
    for (Int_t k=0; k<npar; k++)
    {
        const Double_t sq = fGEma[j][j]*fGEma[k][k];
        fGCorr[j][k] = sq==0 ? 0 : fGEma[j][k] / sqrt( fGEma[j][j]*fGEma[k][k] );
    }
  }


  //--------------------------------------------------
  // reset the errors of the points in the histogram
  for (Int_t i=1; i<=nbins; i++)
    fHist->SetBinError(i, saveError[i-1]);

  return kTRUE;

}



// --------------------------------------------------------------------------
//
//  DetExcessONOFF
//
//  using the result of the polynomial fit (fValuesOFF), DetExcessONOFF determines
//
//  - the total number of events in the signal region (fNon)
//  - the number of backgound events (fitted) in the signal region (fNoffSigFitted)
//  - the number of OFF (normalized evetns) in the alpha OFF distribution (fNoffSig)
//  - the number of excess events (fNexONOFF)
//  - the number of excess events using the fitted OFF (fNexONOFFFitted)
//  - the effective number of background events (fNoff), and fGamma :
//    fNoffSigFitted = fGamma * fNoff;  fdNoffSigFitted = fGamma * sqrt(fNoff);
//
//  It assumed that the polynomial is defined as
//               a0 + a1*x + a2*x**2 + a3*x**3 + ..
//
//  and that the alpha distribution has the range 0 < alpha < 90 degrees
//

Bool_t MHFindSignificanceONOFF::DetExcessONOFF()
{
  //*fLog << "MHFindSignificanceONOFF::DetExcessONOFF;" << endl;

//--------------------------------------------
  // calculate the total number of events (fNon) in the signal region

  fNon  = 0.0;
  fdNon = 0.0;

  Double_t alphaup = -1000.0; 
  Double_t binwidth = fHist->GetBinWidth(1);

  Int_t nbins = fHist->GetNbinsX();
  for (Int_t i=1; i<=nbins; i++)
  {
    Double_t  xlo = fHist->GetBinLowEdge(i);
    Double_t  xup = fHist->GetBinLowEdge(i+1);
 
    // bin must be completely contained in the signal region
    if ( xlo <= (fAlphasig+fEps)  &&  xup <= (fAlphasig+fEps)    )
    {
      Double_t width = fabs(xup-xlo);
      if (fabs(width-binwidth) > fEps)
      {
        *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alpha plot has variable binning, which is not allowed" 
          << endl;
        return kFALSE;
      }

      if (xup > alphaup)
        alphaup = xup;

      fNon  += fHist->GetBinContent(i);
      fdNon += fHist->GetBinError(i) * fHist->GetBinError(i);
    }
  }
  fdNon = sqrt(fdNon);

  *fLog << "MHFindSignificamceONOFF::DetExcessONOFF:" << endl
	<< "fNon = " << fNon << " ; fdNon = " << fdNon << endl;


  // tmp
  //if (fNon == 0)
  //{
  //    cout << "ERROR... WITH COUNTED OF ON EVETNS: fAlphasig, fEps = "
//	   << fAlphasig << ", " << fEps << endl;
  //    fHist-> DrawCopy();
  //   gPad -> SaveAs("HistON.ps");
  //}
  // endtmp


  // the actual signal range is :
  if (alphaup == -1000.0)
    return kFALSE;

  fAlphasi = alphaup;


  //*fLog << "fAlphasi, fNon, fdNon, binwidth, fDegree = " << fAlphasi << ",  "
  //      << fNon << ",  " << fdNon << ",  " << binwidth << ",  "
  //      << fDegree << endl;


  // Calculate the number of OFF events in the signal region
  // ___________________________________________________
  

  fNoffSig = 0.0;
  fdNoffSig = 0.0;

  Double_t alphaupOFF = -1000.0; 
  Double_t binwidthOFF = fHistOFF->GetBinWidth(1);

  Int_t nbinsOFF = fHistOFF->GetNbinsX();

  for (Int_t i=1; i<=nbinsOFF; i++)
  {
    Double_t xlo = fHistOFF->GetBinLowEdge(i);
    Double_t xup = fHistOFF->GetBinLowEdge(i+1);
 
    // bin must be completely contained in the signal region
    if ( xlo <= (fAlphasig+fEps)  &&  xup <= (fAlphasig+fEps)    )
    {
      Double_t width = fabs(xup-xlo);
      if (fabs(width-binwidthOFF) > fEps)
      {
        *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alphaOFF plot has variable binning, which is not allowed" 
          << endl;
        return kFALSE;
      }

      if (xup > alphaupOFF)
        alphaup = xup;

      fNoffSig  += fHistOFF->GetBinContent(i);
      fdNoffSig += fHistOFF->GetBinError(i) * fHistOFF->GetBinError(i);
    }
  }
  fdNoffSig = sqrt(fdNoffSig);

  // tmp
  //if (fNoffSig == 0)
  //{
  //   cout << "ERROR... WITH COUNTED OF OFF EVETNS: fAlphasig, fEps = "
//	   << fAlphasig << ", " << fEps << endl;
  //     fHistOFF-> DrawCopy();
  //   gPad -> SaveAs("HistOFF.ps");
  //}
  //endtmp

  // the actual signal range is :
  if (alphaup == -1000.0)
    return kFALSE;

  fAlphasiOFF = alphaup;

  if (fabs(fAlphasiOFF - fAlphasi) > fEps)
  {
      *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; fAlphasiOFF ("
	    << fAlphasiOFF << ") is not equal to fAlphasi (" 
	    << fAlphasi << "), this is something that should not happen" 
	    << endl;

    //return kFALSE; It might happen in pathological cases (very few OFF)
      // and I want to see the results, anyhow
  }



   // Calculate the number of OFF events in the total OFF region 
  // defined by fAlphaminOFF and fAlphamaxOFF 
  // ___________________________________________________
  
  fNoffTot = 0.0;
  fdNoffTot = 0.0;

  
 
  for (Int_t i=1; i<=nbinsOFF; i++)
  {
    Double_t xlo = fHistOFF->GetBinLowEdge(i);
    Double_t xup = fHistOFF->GetBinLowEdge(i+1);
 
    // bin must be completely contained in the signal region
    if ( xlo >= (fAlphaminOFF-fEps)  &&  xup <= (fAlphamaxOFF+fEps)    )
    {
      Double_t width = fabs(xup-xlo);
      if (fabs(width-binwidthOFF) > fEps)
      {
        *fLog << "MHFindSignificanceONOFF::DetExcessONOFF; alphaOFF plot has variable binning, which is not allowed" 
          << endl;
        return kFALSE;
      }

      fNoffTot  += fHistOFF->GetBinContent(i);
      fdNoffTot += fHistOFF->GetBinError(i) * fHistOFF->GetBinError(i);
    }
  }
  fdNoffTot = sqrt(fdNoffTot);









  //--------------------------------------------
  // calculate the number of OFF fitted events (fNoffSigFitted) in the signal region
  // and its error (fdNoffSigFitted) 

 

 if (fUseFittedQuantities) 
 {
     //--------------------------------------------
     // calculate the number of OFF fitted events (fNoffSigFitted) in the signal region
     // and its error (fdNoffSigFitted) 
     
     Double_t fac = 1.0/binwidthOFF;
     
     fNoffSigFitted        = 0.0;
     Double_t altothejplus1 = fAlphasi; // Limit for signal found for ON data is used. 
     for (Int_t j=0; j<=fDegreeOFF; j++)
     {
	 fNoffSigFitted += fValuesOFF[j] * altothejplus1 / ((Double_t)(j+1));
	 altothejplus1 *= fAlphasi;
     }
     fNoffSigFitted *= fac;
     
     // derivative of fNoffSigFitted 
     Double_t facj;
     Double_t fack;
     
     Double_t sum = 0.0;
     altothejplus1 = fAlphasi;
     for (Int_t j=0; j<=fDegreeOFF; j++)
     {
	 facj = altothejplus1 / ((Double_t)(j+1));
	 
	 Double_t altothekplus1 = fAlphasi;    
	 for (Int_t k=0; k<=fDegreeOFF; k++)
	 {
	     fack = altothekplus1 / ((Double_t)(k+1));
	     sum   += facj * fack * fEmaOFF[j][k];
	     altothekplus1 *= fAlphasi;
	 }
	 altothejplus1 *= fAlphasi;
     }
     sum  *= fac*fac;
     
     if (sum < 0.0)
     {
	 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; error squared is negative" 
	       << endl;
	 return kFALSE;
     }
     
     fdNoffSigFitted = sqrt(sum);
     
     
     // We can now compare fNoffSig with fNoffSigFitted (and their errors)
     // NUmbers should agree within 10 % (errors within 20%)
     
     if (fabs(fNoffSig - fNoffSigFitted) > 0.1 * fNoffSigFitted)
     {
	 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; number of OFF events and Fitted number of OFF events in signal region do not agree (within 10 %)" << endl;
	 
	 *fLog << "fNoffSig = " << fNoffSig << " ; fNoffSigFitted = " << fNoffSigFitted << endl;
	 
	 
	 // return kFALSE; NOt yet...
     }
     
     /*
       if (fabs(fdNoffSig - fdNoffSigFitted) > 0.2 * fdNoffSigFitted)
       {
       *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; error in number of OFF events and error in Fitted number of OFF events in signal region do not agree (within 20 %)" 
       << endl;
       
       
       *fLog << "fdNoffSig = " << fdNoffSig << " ; fdNoffSigFitted = " << fdNoffSigFitted << endl;
       
       
       //return kFALSE; NOt yet...
       }
       
     */
     
    
     
     // Calculate the number of OFF events in the whole fit region (fAlphaminOFF-fAlphamaxOFF)
     // ___________________________________________________
     
     
     
     fNoffTotFitted = 0.0;
     
     altothejplus1 = fAlphamaxOFF; // Limit for OFF data fit 
     for (Int_t j=0; j<=fDegreeOFF; j++)
     {
	 fNoffTotFitted += fValuesOFF[j] * altothejplus1 / ((Double_t)(j+1));
	 altothejplus1 *= fAlphamaxOFF;
     }
     fNoffTotFitted *= fac;
     
     // derivative of fdNoffTotFitted 
     
     
     sum = 0.0;
     altothejplus1 = fAlphamaxOFF;
     for (Int_t j=0; j<=fDegreeOFF; j++)
     {
	 facj = altothejplus1 / ((Double_t)(j+1));
	 
	 Double_t altothekplus1 = fAlphamaxOFF;    
	 for (Int_t k=0; k<=fDegreeOFF; k++)
	 {
	     fack = altothekplus1 / ((Double_t)(k+1));
	     
	     sum   += facj * fack * fEmaOFF[j][k];
	     altothekplus1 *= fAlphamaxOFF;
	 }
	 altothejplus1 *= fAlphamaxOFF;
     }
     sum  *= fac*fac;

     if (sum < 0.0)
     {
	 *fLog << "MHFindsignificanceONOFF::DetExcessONOFF; error squared is negative" 
	       << endl;
	 return kFALSE;
     }
     
     fdNoffTotFitted = sqrt(sum);
     


 }

else
 {
     fNoffSigFitted = 0.0;
     fdNoffSigFitted = 0.0;
     fNoffTotFitted = 0.0;
     fdNoffTotFitted = 0.0;
 }



 *fLog << "MHFindSignificamceONOFF::DetExcessONOFF; INFO ABOOUT COMPUTED OFF EVENTS..." << endl
       << "fNoffSig = " << fNoffSig << "; fdNoffSig =  " << fdNoffSig << endl
       << "fNoffSigFitted = " << fNoffSigFitted 
       << "; fdNoffSigFitted =  " << fdNoffSigFitted << endl;
 




  //--------------------------------------------
  // calculate the number of excess events in the signal region

  fNexONOFF = fNon - fNoffSig*fNormFactor;
  fNexONOFFFitted = fNon - fNoffSigFitted*fNormFactor;
  
  
  *fLog << "MHFindSignificamceONOFF::DetExcessONOFF;" << endl
	<< "fNexONOFF (= fNon - fNoffSig*fNormFactor)  = " << fNexONOFF << endl
	<< "fNexONOFFFitted (= fNon - fNoffSigFitted*fNormFactor)  = " << fNexONOFFFitted << endl;


  //--------------------------------------------
  // calculate the effective number of background events (fNoff) , and fGamma :
  // fNbg = fGamma * fNoff = fNormFactor* fNoffSigFitted;   
  // dfNbg = fGamma * sqrt(fNoff) = fNormFactor * fdNoffSigFitted;
  
  if (fNoffSigFitted < 0.0)
  {
      *fLog << "MHFindSignificamceONOFF::DetExcessONOFF; number of fitted OFF events in signal region is negative,  fNoffSigFitted, fdNoffSigFitted = "
	    << fNoffSigFitted   << ",  " << fdNoffSigFitted << endl;
      
      fGamma = 1.0;
      fNoff  = 0.0;
      return kFALSE;
  }
  
  if (fNoffSigFitted > 0.0)
  {
      fGamma = fNormFactor * fdNoffSigFitted*fdNoffSigFitted/fNoffSigFitted;
      fNoff  = fNormFactor * fNoffSigFitted/fGamma; 
  }
  else
  {
      fGamma = 1.0;
      fNoff  = 0.0;
  }
  
   *fLog << "MHFindSignificamceONOFF::DetExcessONOFF: " << endl
	 << "fGamma = " << fGamma << " ; fNoff = " << fNoff << endl;

  
  return kTRUE;
}





// --------------------------------------------------------------------------
//
//  SigmaLiMa
//
//  calculates the significance according to Li & Ma
//  ApJ 272 (1983) 317; formula 17
//
Bool_t MHFindSignificanceONOFF::SigmaLiMa(Double_t non,   Double_t noff, 
                                     Double_t gamma, Double_t *siglima)
{
  if (gamma <= 0.0  ||  non <= 0.0  ||  noff <= 0.0)
  {
    *siglima = 0.0;
    return kFALSE;
  }

  Double_t help1 = non  * log( (1.0+gamma)*non  / (gamma*(non+noff)) );
  Double_t help2 = noff * log( (1.0+gamma)*noff / (       non+noff ) );
  *siglima = sqrt( 2.0 * (help1+help2) );

  Double_t nex = non - gamma*noff;
  if (nex < 0.0)
    *siglima = - *siglima;

  //*fLog << "MHFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = "
  //      << non << ",  " << noff << ",  " << gamma << ",  " << *siglima << endl;

  return kTRUE;
}



//  calculates the significance according to Li & Ma
//  ApJ 272 (1983) 317; formula 5
//
Bool_t MHFindSignificanceONOFF::SigmaLiMaForm5(Double_t non,   Double_t noff, 
					       Double_t gamma, Double_t *siglima)
{
  if (gamma <= 0.0  ||  non <= 0.0  ||  noff <= 0.0)
  {
    *siglima = 0.0;
    return kFALSE;
  }

  Double_t nex = non - (gamma*noff);
  Double_t tmp = non + (gamma*gamma)*noff;
  tmp = TMath::Sqrt(tmp);

  *siglima = nex/tmp;

  if (nex < 0.0)
    *siglima = - *siglima;

  //*fLog << "MHFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = "
  //      << non << ",  " << noff << ",  " << gamma << ",  " << *siglima << endl;

  return kTRUE;
}




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

// Following function computes a clone of fHistOFF and normalizes 
// contents, errors and fPolyOFF (if exists) with the fNormFactor. 
// This normalized OFF hist will be used when plotting OFF data 
// together with ON data.

Bool_t MHFindSignificanceONOFF::ComputeHistOFFNormalized()
{

    
    if (!fHist)
    {
	*fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fHist does not exist, normalization of HistOFF can not be performed properly..."
	      << endl;
	return kFALSE;
	
    }

    

    if (!fHistOFF)
    {
	*fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fHistOFF does not exist, hence can not be normalized" 
          << endl;
    return kFALSE;
	
    }


    if (fNormFactor <= 0)
    {
	*fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fNormFactor is ZERO or NEGATIVE, it might have not been defined yet..."
          << endl;
    return kFALSE;

    }


    Double_t BinWidthAlphaON = fHist -> GetBinWidth(1);
    Double_t BinWidthAlphaOFF = fHistOFF -> GetBinWidth(1);
    Double_t BinWidthRatioONOFF = BinWidthAlphaON/BinWidthAlphaOFF;


    *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; INFO about alpha ON, OFF histo bins"
	  << endl
	//  << "fHist bin width = " << BinWidthAlphaON 
	  << " fHistOFF bin width = " << BinWidthAlphaOFF << endl;



    TString fHistOFFNormalizedName = fHistOFF -> GetName();
    fHistOFFNormalizedName += (" (Normalized)");
    // fHistOFFNormalized = (TH1*) fHistOFF -> Clone();
    // fHistOFFNormalized -> SetNameTitle(fHistOFFNormalizedName, fHistOFFNormalizedName);


    Int_t nbinsOFFNormalized = 0;
    Int_t nbinsOFF = 0;
    Double_t xlow = 0.0;
    Double_t xup = 0.0;
    Double_t content = 0.0;
    Double_t error = 0.0;
    Double_t BinCenter = 0.0;
    

    // Bins for normalized OFF histo will be the ones of ON histo

    
    nbinsOFF = fHistOFF -> GetNbinsX();
    nbinsOFFNormalized = nbinsOFF;
    xlow = fHistOFF -> GetBinLowEdge(1);
    xup = fHistOFF -> GetBinLowEdge(nbinsOFFNormalized);
    xup = xup + BinWidthAlphaON;

    *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; Limits for fHistOFFNormalized: " 
	  << "nbins, xlow, xup: " << nbinsOFFNormalized << ", " 
	  <<  xlow << ", " << xup << endl;
    
    fHistOFFNormalized = new TH1F (fHistOFFNormalizedName, fHistOFFNormalizedName, 
				  nbinsOFFNormalized, xlow, xup);
    

    // fHistOFFNormalized is filled with data from fHistOFF, 
    // taken into account the possible bin difference

    
    for (Int_t i = 1; i <= nbinsOFF; i++)
    {
	BinCenter = fHistOFF -> GetBinCenter(i);
	fHistOFFNormalized -> Fill (BinCenter, fHistOFF -> GetBinContent(i));
	fHistOFFNormalized -> SetBinError(i, fHistOFF -> GetBinError(i));
    }
	    

    

    for (Int_t i = 1; i <= nbinsOFFNormalized; i++)
    {
	content = fNormFactor * fHistOFFNormalized -> GetBinContent(i);
	error = fNormFactor * fHistOFFNormalized -> GetBinError(i);
	
	fHistOFFNormalized -> SetBinContent (i, content);
	fHistOFFNormalized -> SetBinError (i, error);
    }
	

    // Number of entries is obtained from histOFF.
    // and set to histOFFNoramlized; otherwise, the number 
    // of entries in  histOFFNoramlized would be "nbins"

    Double_t entries = fNormFactor * (fHistOFF -> GetEntries());
    fHistOFFNormalized -> SetEntries(entries);




    // If polynomial fit has been performed for fHistOFF, 
    // it is defined a new polyfunction for fHistOFFNormalized, 
    // which will be the polyfunction of fHistOFF normalized.
    // Function will be added to the function list of fHistOFFNormalized


    
    if (fPolyOFF == NULL)
    {
	*fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fPolyOFF does not exist..."
	      << endl;
    }

    if (fPolyOFF)
    {
	*fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; fPolyOFF exists... will be also normalized and included in the list of functions of fHistOFFNormalized"
	      << endl;
	
	

	
	// Normalization of the function using fNormFactor and 
	// BinWidthON/BinWidthOFF relation of alpha ON/OFF histograms
	// This makes possible to plot it together with ON alpha histo

	TString FunctionName("PolyOFFNormalized");

	Double_t xmin = fAlphaminOFF;
	Double_t xmax = fAlphamaxOFF;

	TString formula = "[0]";
	TString bra1     = "+[";
	TString bra2     =    "]";
	TString xpower   = "*x";
	TString newpower = "*x";
	for (Int_t i=1; i<=fDegreeOFF; i++)
	{
	    formula += bra1;
	    formula += i;
	    formula += bra2;
	    formula += xpower;
	    
	    xpower += newpower;
	}

      *fLog << "MHFindsignificanceONOFF::ComputeHistOFFNormalized; formula = " << formula << endl;


	
	fPolyOFFNormalized = new TF1 (FunctionName, formula, xmin, xmax);
	
	
	Double_t Parameter = 0.0;
	Double_t ParameterError = 0.0;

	*fLog << " MHFindsignificanceONOFF::ComputeHistOFFNormalized; Fit parameters info: " << endl;
	for (Int_t i = 0; i <= fDegreeOFF; i++)
	{
	    Parameter = fNormFactor * BinWidthRatioONOFF * fValuesOFF[i];
	    ParameterError = fNormFactor * BinWidthRatioONOFF * fErrorsOFF[i];

	    fPolyOFFNormalized -> SetParameter(i, Parameter);
	    fPolyOFFNormalized -> SetParError(i,ParameterError);

	   

	    // Parameters are shown :
	    
	    *fLog  << " fValuesOFF[" << i<< "] = " << fValuesOFF[i] 
		   << " ; Parameter for  fPolyOFFNormalized = " << Parameter  << endl;

	}
	    

	TList *funclist = fHistOFFNormalized->GetListOfFunctions();
	
	// temporal...
	//*fLog << "INFO concerning list of functions of fHistOFFNormalized :" << endl
	//       << "List before adding OFF Normal., after adding it and after removing fPolyOFF..." 
	//       << endl;
	
	//funclist-> Print();
	funclist-> Add(fPolyOFFNormalized);
	
	//funclist-> Print();
	
		

    }

    return kTRUE;

}

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

Bool_t MHFindSignificanceONOFF::DrawHistOFF()
{
    if (fHistOFF == NULL )
    {
	*fLog << "MHFindSignificanceONOFF::DrawHistOFF; fHistOFF = NULL" << endl;
	return kFALSE;
    }

//    fPsFilename -> NewPage();

    // PLOT DISABLE
    /*
    TCanvas* CanvasHistOFF = new TCanvas(fHistOFF->GetName(), fHistOFF->GetName(), 600, 600);

    //gStyle->SetOptFit(1011);

    gROOT->SetSelectedPad(NULL);    
    gStyle->SetPadLeftMargin(0.1);
    gStyle -> SetOptStat(1);

    CanvasHistOFF->cd();


    if (fHistOFF)
    {
      fHistOFF->DrawCopy();
    }

    // TF1 *fpolyOFF = fHistOFF->GetFunction("PolyOFF");    
    if (fPolyOFF == NULL)
	*fLog << "MHFindSignificanceONOFF::DrawHistOFF; fpolyOFF = NULL" << endl;

    if (fPolyOFF)
    {
      // 2, 1 is red and solid
      fPolyOFF->SetLineColor(2);
      fPolyOFF->SetLineStyle(1);
      fPolyOFF->SetLineWidth(2);
      fPolyOFF->DrawCopy("same");
    }

    CanvasHistOFF -> Update();

    
    */

    return kTRUE;
}


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

Bool_t MHFindSignificanceONOFF::DrawHistOFFNormalized()
{
    if (fHistOFFNormalized == NULL )
    {
	*fLog << "MHFindSignificanceONOFF::DrawHistOFF; fHistOFFNormalized = NULL" << endl;
	return kFALSE;
    }

    //   fPsFilename -> NewPage();

    // PLOT DISABLE TO PERFORM GRID ANALYSIS
    /*
    TCanvas* CanvasHistOFFNormalized = new TCanvas(fHistOFFNormalized->GetName(), 
						   fHistOFFNormalized->GetName(), 600, 600);

    //gStyle->SetOptFit(1011);

    // gROOT->SetSelectedPad(NULL);    
    gStyle->SetPadLeftMargin(0.1);
    gStyle -> SetOptStat(1);

    CanvasHistOFFNormalized->cd();
    

    if (fHistOFFNormalized)
    {
      fHistOFFNormalized->DrawCopy();
    }       

    // TF1 *fpolyOFFNormalized = fHistOFFNormalized->GetFunction("PolyOFFNormalized");    
    if (fPolyOFFNormalized == NULL)
	*fLog << "MHFindSignificanceONOFF::DrawHistOFF; fPolyOFFNormalized = NULL" << endl;

    if (fPolyOFFNormalized)
    {
      // 2, 1 is red and solid
      fPolyOFFNormalized->SetLineColor(2);
      fPolyOFFNormalized->SetLineStyle(1);
      fPolyOFFNormalized->SetLineWidth(2);
      fPolyOFFNormalized->DrawCopy("same");
    }

    CanvasHistOFFNormalized -> Update();
*/

    return kTRUE;
    
}


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

Bool_t MHFindSignificanceONOFF::DrawFit(const Option_t *opt)
{
    if (fHistOFFNormalized == NULL ||  fHist == NULL)
      *fLog << "MHFindSignificanceONOFF::DrawFit; fHistOFFNormalized = NULL or fHist == NULL" << endl;

    //fPsFilename -> NewPage();


    // PLOT DISABLE TO PERFORM GRID ANALYSIS
    // I DO SAVE PS FILE, BUT CANVAS IS DELETED AFTERWARDS

    fCanvas = new TCanvas("Alpha", "Alpha plot", 600, 600);
    fCanvas -> SetFillColor(10);
    



    //gStyle->SetOptFit(1011);

    gROOT->SetSelectedPad(NULL);    
    gStyle -> SetFrameFillColor(10);
    gStyle->SetPadLeftMargin(0.15);
    gStyle -> SetOptStat(1);
    
    fCanvas->cd();
    
    if (fHist)
      {
	fHist -> SetTitle("Alpha Plot");
	fHist-> SetTitleOffset(1.5, "Y");
        fHist-> SetLineColor(kRed);
	fHist-> DrawCopy();
	
      }


    if (fHistOFFNormalized)
    {
	TF1 *fpoly = fHistOFFNormalized->GetFunction("PolyOFFNormalized");    
	if (fpoly == NULL)
	    *fLog << "MHFindSignificanceONOFF::DrawFit; fPolyOFFNormalized = NULL" << endl;

	if (fpoly)
	{
	    // 2, 1 is red and solid
	    fpoly->SetLineColor(2);
	    fpoly->SetLineStyle(1);
	    fpoly->SetLineWidth(2);
	    fpoly->DrawCopy("same");
	}


//	fHistOFFNormalized->SetLineColor(kRed);
//	fHistOFFNormalized->DrawCopy("same");

    }

    if (fFitGauss)
    {
      TF1 *fpolygauss = fHist->GetFunction("PolyGauss");    
      if (fpolygauss == NULL)
        *fLog << "MHFindSignificanceONOFF::DrawFit; fpolygauss = NULL" << endl;

      if (fpolygauss)
      {
        // 4, 1 is blue and solid
        fpolygauss->SetLineColor(4);
        fpolygauss->SetLineStyle(1);
        fpolygauss->SetLineWidth(4);
        fpolygauss->DrawCopy("same");
      }
 
      TF1 *fbackg = fHist->GetFunction("Backg");    
      if (fbackg == NULL)
        *fLog << "MHFindSignificanceONOFF::DrawFit; fbackg = NULL" << endl;

    if (fHistOFFNormalized)
      {
      fHistOFFNormalized->SetFillStyle(3002);
      fHistOFFNormalized->SetFillColor(9);
      fHistOFFNormalized->DrawCopy("Psamebar");}


      if (fbackg)
      {
        // 10, 1 is white and solid
        //fbackg->SetLineColor(kRed);
        fbackg->SetLineColor(196);
        fbackg->SetLineStyle(1);
        fbackg->SetLineWidth(4);
	fbackg->DrawCopy("same"); // I do not want to draw it... already too many things.
      }
    }

//    fCanvas->SaveAs("plot.root");
//    fCanvas->SaveAs("plot.ps");

    //-------------------------------
    // print results onto the figure
    
    
    
    TPaveText *pt = new TPaveText(0.30, 0.35, 0.70, 0.90, "NDC");
    char tx[100];
    
    sprintf(tx, "Results of polynomial fit to OFF (order %2d) :", fDegreeOFF);
    TText *t1 = pt->AddText(tx);
    t1->SetTextSize(0.03);
    t1->SetTextColor(2);
    
    sprintf(tx, "   (%6.2f< |alpha| <%6.2f [\\circ])", fAlphaminOFF, fAlphamaxOFF);
    pt->AddText(tx);
    
    sprintf(tx, "   chi2 = %8.2f,  Ndof = %4d,  Prob = %6.2f", 
	    fChisqOFF, fNdfOFF, fProbOFF);
    pt->AddText(tx);
    
    sprintf(tx, "   OFF events (fit)= %8.1f #pm %8.1f", 
	    fNoffTotFitted, fdNoffTotFitted);
    pt->AddText(tx);
    
    sprintf(tx, "   OFF events (meas) = %8.1f #pm %8.1f", fNoffTot, fdNoffTot);
    pt->AddText(tx);
    
    sprintf(tx, "   OFF Normalization Factor (= Non/Noff) = %4.4f", fNormFactor);
    pt->AddText(tx);
    
    
    
    
    //sprintf(tx, "     ");
    //pt->AddText(tx);
    
    //--------------
    sprintf(tx, "Results for |alpha|< %6.2f [\\circ] :", fAlphasi);
    TText *t6 = pt->AddText(tx);
    t6->SetTextSize(0.03);
    t6->SetTextColor(8);
    
    sprintf(tx, "   Non = %8.1f #pm %8.1f", fNon, fdNon);
    pt->AddText(tx);
    
    
    if(fUseFittedQuantities)
      {
	//// **************************************************
	///// PRINT INFORMATION ABOUT FITTED QUANTITIES  /////////
	
	
	
	Double_t NoffFitNormalized = fNoffSigFitted * fNormFactor;
	Double_t ErrorNoffFitNormalized = fdNoffSigFitted * fNormFactor;
	Double_t SignificanceUsed = GetSignificance();
	
	sprintf(tx, "   Noff Fitted (Normalized)  = %8.1f #pm %8.1f", 
		NoffFitNormalized, ErrorNoffFitNormalized);
	pt->AddText(tx);
	
	
	sprintf(tx, "   Nex (ON - OFF Fitted) = %8.1f #pm %8.1f", 
		fNexONOFFFitted, fdNexONOFFFitted);
	pt->AddText(tx);
	
	
	sprintf(tx, "   Gamma = %4.4f,   Effective Noff (i.e. fNoff) = %6.1f", 
		fGamma, fNoff);
	pt->AddText(tx);
	
	
	Double_t ratio = fNoffSigFitted>0.0 ? fNexONOFFFitted/(fNoffSigFitted*fNormFactor) : 0.0;
	sprintf(tx, "   Significance = %6.2f,    Nex/(Nbg*NormFactor) = %6.2f", 
		SignificanceUsed, ratio);
	pt->AddText(tx);
	
	
      }
    
    else
      {
	//// **************************************************
	///// PRINT INFORMATION ABOUT MEASURED QUANTITIES  /////////
	
	
	Double_t NoffNormalized = fNoffSig * fNormFactor;
	Double_t ErrorNoffNormalized = fdNoffSig * fNormFactor;
	Double_t SignificanceUsed = GetSignificance();
	
	sprintf(tx, "   Noff measured (Normalized)  = %8.1f #pm %8.1f", 
		NoffNormalized, ErrorNoffNormalized);
	pt->AddText(tx);
	
	sprintf(tx, "   Nex (ON - OFF measured) = %8.1f #pm %8.1f", 
		fNexONOFF, fdNexONOFF);
	pt->AddText(tx);
	
	Double_t ratio = fNoffSig>0.0 ? fNexONOFF/(fNoffSig*fNormFactor) : 0.0;
	sprintf(tx, "   Significance = %6.2f,    Nex/(Nbg*NormFactor) = %6.2f", 
		SignificanceUsed, ratio);
	pt->AddText(tx);
	
      }
    
    /*
      // Temporally I will also show ALL SIGMALIMA COMPUTED.
      
      sprintf(tx, 
      "  fSigLiMa1 =  %6.2f, fSigLiMa2 =  %6.2f, fSigLiMa3 =  %6.2f", 
      fSigLiMa,fSigLiMa2, fSigLiMa3);
      pt->AddText(tx);
    */
    
    
    //--------------
    if (fFitGauss)
      {
	sprintf(tx, "Results of (polynomial+Gauss) fit  :");
	TText *t7 = pt->AddText(tx);
	t7->SetTextSize(0.03);
	t7->SetTextColor(4);
	
	sprintf(tx, "   chi2 = %8.2f,  Ndof = %4d,  Prob = %6.2f", 
		fGChisq, fGNdf, fGProb);
	pt->AddText(tx);
	
	sprintf(tx, "   Sigma of Gauss = %8.1f #pm %8.1f  [\\circ]", 
		fSigmaGauss, fdSigmaGauss);
	pt->AddText(tx);
	
	sprintf(tx, "   total no.of excess events = %8.1f #pm %8.1f", 
		fNexGauss, fdNexGauss);
	pt->AddText(tx);
      }
    //--------------
    
    pt->SetFillStyle(0);
    pt->SetBorderSize(0);
    pt->SetTextAlign(12);
    
    
    if(fPrintResultsOntoAlphaPlot)
      {
	pt->Draw();
      }
    fCanvas->Modified();
    fCanvas->Update();
    
    //    fPsFilename -> NewPage();



    if (fSavePlots)
    {
	// ********************************************
	// TMP solution while the TPostScript thing is not working. 
	// PsFileName for storing these histograms is derived 
	// from fPsFilenameString. 
	
	
	cout << "Alpha plot with ON-OFF data will be saved in PostScript file  " ;
	
	
	if (!fPsFilenameString.IsNull())
	{
	    TString filename  = (fPsFilenameString);
	    TString filename2 = (fPsFilenameString);
	    // Train or Test Sample is specified outside 
	    // class MHFindSignificanceONOFF, and included in 
	    // fPsFilenameString
	    
	    filename  += ("AlphaPlot.ps");
	    filename2 += ("AlphaPlot.C");
	    cout << filename << endl;	
	    fCanvas -> SaveAs(filename);
	    fCanvas -> SaveAs(filename2);
	}
	
	// END OF TEMPORAL SOLUTION
	// ********************************************
	
    }
    
    // Canvvas deleted to allow for GRID analysis

    delete fCanvas;
	

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Print the results of the polynomial fit to the alpha OFF distribution
// 
//
void MHFindSignificanceONOFF::PrintPolyOFF(Option_t *o) 
{
  *fLog << "---------------------------" << endl;
  *fLog << "MHFindSignificanceONOFF::PrintPolyOFF :" << endl; 

  *fLog << "fAlphaminOFF, fAlphamaxOFF, fDegreeOFF "
        << fAlphaminOFF << ",  " << fAlphamaxOFF << ",  " << fDegreeOFF <<  endl;

  *fLog << "fMbinsOFF, fNzeroOFF, fIstatOFF = " << fMbinsOFF << ",  "
        << fNzeroOFF << ",  " << fIstatOFF << endl;

  *fLog << "fChisqOFF, fNdfOFF, fProbOFF = " << fChisqOFF << ",  " 
        << fNdfOFF << ",  " << fProbOFF << endl; 

  *fLog << "fNon; fNoffSigFitted, fdNoffSigFitted; fNoffSig, fdNoffSig = "
        << fNon << ";  " << fNoffSigFitted << ",  " << fdNoffSigFitted 
	<< ";  " << fNoffSig << ",  " << fdNoffSig <<  endl;

  Double_t sigtoback = fNoffSigFitted >0.0 ? fNexONOFFFitted/(fNoffSigFitted*fNormFactor) : 0.0;

  
  *fLog << "fNexONOFFFitted, fdNexONOFFFitted, fGamma, fNoff, fSigLiMa, sigtoback = "
        << fNexONOFFFitted << ",  " << fdNexONOFFFitted << ",  " 
	<< fGamma << ",  " << fNoff 
        << ",  " << fSigLiMa << ",  "
	<< sigtoback << endl;


  Double_t sigtoback2 = fNoffSig >0.0 ? fNexONOFF/(fNoffSig*fNormFactor) : 0.0;
  
  *fLog << "fNexONOFF, fdNexONOFF, fNormFactor, fNoffSig, fSigLiMa2, sigtoback2 = "
        << fNexONOFF << ",  " << fdNexONOFF << ",  " 
	<< fNormFactor << ",  " << fNoffSig 
        << ",  " << fSigLiMa2 << ",  "
	<< sigtoback2 << endl;

 *fLog << "---------------------------" << endl;
}

// --------------------------------------------------------------------------
//
// Print the results of the polynomial fit to the alpha distribution
// 
//
void MHFindSignificanceONOFF::PrintPoly(Option_t *o) 
{
  *fLog << "---------------------------" << endl;
  *fLog << "MHFindSignificanceONOFF::PrintPoly :" << endl; 

  *fLog << "fAlphami, fAlphama, fDegree, fAlphasi = "
        << fAlphami << ",  " << fAlphama << ",  " << fDegree << ",  " 
        << fAlphasi << endl;

  *fLog << "fMbins, fNzero, fIstat = " << fMbins << ",  "
        << fNzero << ",  " << fIstat << endl;

  *fLog << "fChisq, fNdf, fProb = " << fChisq << ",  " 
        << fNdf << ",  " << fProb << endl; 

  *fLog << "fNon, fNbg, fdNbg, fNbgtot, fNbgtotFitted, fdNbgtotFitted = "
        << fNon << ",  " << fNbg << ",  " << fdNbg << ",  " << fNbgtot 
        << ",  " << fNbgtotFitted << ",  " << fdNbgtotFitted << endl;

  Double_t sigtoback = fNbg>0.0 ? fNex/fNbg : 0.0;
  *fLog << "fNex, fdNex, fGamma, fNoff, fSigLiMa, sigtoback = "
        << fNex << ",  " << fdNex << ",  " << fGamma << ",  " << fNoff 
        << ",  " << fSigLiMa << ",  " << sigtoback << endl;

  //------------------------------------
  // get errors

  /*
  Double_t eplus; 
  Double_t eminus; 
  Double_t eparab; 
  Double_t gcc;
  Double_t errdiag;


  if ( !fConstantBackg )
  {
    *fLog << "parameter value     error     eplus     eminus    eparab   errdiag   gcc"
          << endl; 

    for (Int_t j=0; j<=fDegree; j++)
    {
      if (gMinuit)
        gMinuit->mnerrs(j, eplus, eminus, eparab, gcc);
      errdiag = sqrt(fEma[j][j]);
      *fLog << j << "  " << fValues[j] << "  "   << fErrors[j] << "  "
            << eplus     << "  "       << eminus << "  " << eparab     << "  " 
            <<  errdiag  << "  "       << gcc    << endl;
    }
  }  
  else
  {
    *fLog << "parameter value     error     errdiag "
          << endl; 

    for (Int_t j=0; j<=fDegree; j++)
    {
      errdiag = sqrt(fEma[j][j]);
      *fLog << j << "  " << fValues[j] << "  "   << fErrors[j] << "  "
            <<  errdiag  << endl;
    }
  }  
  */

  //----------------------------------------
  /*
  *fLog << "Covariance matrix :" << endl;
  for (Int_t j=0; j<=fDegree; j++)
  {
    *fLog << "j = " << j << " :   ";
    for (Int_t k=0; k<=fDegree; k++)
    {
      *fLog << fEma[j][k] << "   ";
    }
    *fLog << endl;
  }

  *fLog << "Correlation matrix :" << endl;
  for (Int_t j=0; j<=fDegree; j++)
  {
    *fLog << "j = " << j << " :   ";
    for (Int_t k=0; k<=fDegree; k++)
    {
      *fLog << fCorr[j][k] << "   ";
    }
    *fLog << endl;
  }
  */

  *fLog << "---------------------------" << endl;
}

// --------------------------------------------------------------------------
//
// Print the results of the (polynomial+Gauss) fit to the alpha distribution
// 
//
void MHFindSignificanceONOFF::PrintPolyGauss(Option_t *o) 
{
  *fLog << "---------------------------" << endl;
  *fLog << "MHFindSignificanceONOFF::PrintPolyGauss :" << endl; 

  *fLog << "fAlphalo, fAlphahi = "
        << fAlphalo << ",  " << fAlphahi << endl;

  *fLog << "fGMbins, fGNzero, fGIstat = " << fGMbins << ",  "
        << fGNzero << ",  " << fGIstat << endl;

  *fLog << "fGChisq, fGNdf, fGProb = " << fGChisq << ",  " 
        << fGNdf << ",  " << fGProb << endl; 


  //------------------------------------
  // get errors

  Double_t eplus; 
  Double_t eminus; 
  Double_t eparab; 
  Double_t gcc;
  Double_t errdiag;

  *fLog << "parameter value     error     eplus     eminus    eparab   errdiag   gcc"
        << endl; 
  for (Int_t j=0; j<=(fDegree+3); j++)
  {
    if (gMinuit)
      gMinuit->mnerrs(j, eplus, eminus, eparab, gcc);
    errdiag = sqrt(fGEma[j][j]);
    *fLog << j << "  " << fGValues[j] << "  "   << fGErrors[j] << "  "
          << eplus     << "  "       << eminus << "  " << eparab     << "  " 
          <<  errdiag  << "  "       << gcc    << endl;
  }

  
  *fLog << "Covariance matrix :" << endl;
  for (Int_t j=0; j<=(fDegree+3); j++)
  {
    *fLog << "j = " << j << " :   ";
    for (Int_t k=0; k<=(fDegree+3); k++)
    {
      *fLog << fGEma[j][k] << "   ";
    }
    *fLog << endl;
  }

  *fLog << "Correlation matrix :" << endl;
  for (Int_t j=0; j<=(fDegree+3); j++)
  {
    *fLog << "j = " << j << " :   ";
    for (Int_t k=0; k<=(fDegree+3); k++)
    {
      *fLog << fGCorr[j][k] << "   ";
    }
    *fLog << endl;
  }

  *fLog << "---------------------------" << endl;
}

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