/* ======================================================================== *\
!
! *
! * 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): Marcols Lopez,  August 2004      <mailto:marcos@gae.ucm.es>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHMyFindSignificance
//
//
//  Non: # events in the signal region (0,fAlphasig) of the ON-Data hist.
//  Nbg: # events in the signal region (0,fAlphasig) of the OFF-Data hist.
//  Nex = Non - Nbg
//
/////////////////////////////////////////////////////////////////////////////
#include "MHMyFindSignificanceONOFF.h"

#include <TH1.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TPaveText.h>
#include <TGaxis.h>
#include <TStyle.h>
#include <TLine.h>
#include <TString.h>
#include <TGraph.h>
#include <TPaveStats.h>

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


ClassImp(MHMyFindSignificanceONOFF);

using namespace std;


// --------------------------------------------------------------------------
//
//  Constructor
//
MHMyFindSignificanceONOFF::MHMyFindSignificanceONOFF(const char *name, const char *title) 
    : fNon(0.0), fdNon(0.0), fNbg(0.0), fdNbg(0.0), fNex(0.0), fdNex(0.0), 
      fSigLiMa(0.0), 
      fNbgFit(0.0), fdNbgFit(0.0), fNexFit(0.0), fdNexFit(0.0), fGamma(0.0), 
      fNoff(0.0), fSigLiMaFit(0.0),
      fBestSigma(0.0), fBestAlphaCut(0.0), fBestExcess(0.0)
{
    fName  = name  ? name  : "MHMyFindSignificanceONOFF";
    fTitle = title ? title : "Find Significance in alpha plot";
   
    fHistON            = NULL;
    fHistOrigON        = NULL;
    fHistOFF           = NULL;
    fHistOrigOFF       = NULL;
    fHistOFFNormalized = NULL;

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

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




// ---------------------------------------------------------------------------
//
// Calculate  
//
//   - fNon:       number of events in the signal region from ON data 
//
//   and 3 different estimation of the background:
//
//   - fNbgOFF:    number of events in the signal region from  OFF data
//   - fNbgFitOFF: estimation of # events in the signal region from a fit of
//                 the OFF data in the Bkg region.
//   - fNbgFitON:  estimation of # events in the signal region from a fit of
//                 the ON data in the Bkg region.
//   - fNex
//
Bool_t MHMyFindSignificanceONOFF::FindSigmaONOFF(TH1 *histON, TH1 *histOFF, 
         Double_t alphasig, Double_t alphamin, Double_t alphamax, Int_t degree)
{

    fHistOrigON = histON;
    fHistOrigOFF = histOFF;

    fAlphasig = alphasig;  
    fAlphaminON = alphamin;
    fAlphamaxON = alphamax;
    fAlphaminOFF = 0;  // <--------------------
    fAlphamaxOFF = alphamax;
    fDegree   = degree;


 

    // -----------------------------------------------------------------
    // 
    // Calculate the normalization factor in range (alphanin,alphamax)
    //
    // FIXME: This can be done directly with: 
    //        fFindsigON.GetNbtot()/fFindsigOFF.GetNbtot()
    //        but for that the polynomial for OFF-Data has to be fitted
    //        (alphamin, alphamax) instead of (0,alphamax)
    //
    CalcNormalizationFactor(fHistOrigON, fHistOrigOFF, fAlphaminON, fAlphamaxON);


    // ----------------------------------------------------------------------- 
    //
    // Use MHFindSignificance to find the # evts in the signal and Bkg region
    // and to fit the histograms. 
    // The ON hist is fitted to a polynomial in reg. (fAlphaminON,fAlphamaxON)
    // and to a gaussian+polynomial.
    // The OFF hist is only fitted to a polynomial in(fAlphaminOFF,fAlphamaxOFF
    //
    fFindsigON.SetReduceDegree(fReduceDegree);
    fFindsigOFF.SetReduceDegree(fReduceDegree);
    fFindsigON.SetRebin(fRebin);
    fFindsigOFF.SetRebin(fRebin);
 
    const Bool_t rc1 = fFindsigON.FindSigma(fHistOrigON, fAlphaminON, fAlphamaxON, fDegree, fAlphasig, 0,1,0);   
   
    const Bool_t rc2 = fFindsigOFF.FindSigma(fHistOrigOFF, fAlphaminOFF,fAlphamaxOFF, fDegree, fAlphasig, 0,0,0);
  
    if (rc1==kFALSE || rc2==kFALSE)
    {
	*fLog << err << dbginf << "MHFindSignificance retuns kFALSE." << endl;
	return kFALSE;
    }


    // Check consistency
    if (fabs(fFindsigON.GetAlphasi() - fFindsigOFF.GetAlphasi()) > fEps)
    {
	*fLog << "fAlphasiOFF (" <<  fFindsigOFF.GetAlphasi()<< ") is not equal to fAlphasi ("
	    << fFindsigON.GetAlphasi() << "), this is something that should not happen" 
	    << endl;
	//return kFALSE; It might happen in pathological cases (very few OFF)
    }


    //
    // Get results of the fits
    //
    fHistON = fFindsigON.GetHist();  // don't use the original hist to draw
    fHistOFF = fFindsigOFF.GetHist();
    // From ON-Data
    fNon  = fFindsigON.GetNon();
    fdNon = fFindsigON.GetdNon();
    // From OFF-Data
    fNbg     = fFindsigOFF.GetNon();   // counts in signal reg. OFF hist NO 
    fdNbg    = fFindsigOFF.GetdNon();  //normalized !!!
    // From OFF-Data polynomial fit
    fNbgFit  = fFindsigOFF.GetNbg(); 
    fdNbgFit = fFindsigOFF.GetdNbg();
    fGamma   =  fFindsigOFF.GetGamma() * fNormFactor; //<-----------
    fNoff    =  fFindsigOFF.GetNoff();


    // 
    // Compare the counter number of bg events, NbgMea, and the fitted, NbgFit
    //
    if (fabs(fNbg - fNbgFit) > 0.1 * fNbgFit)
    {
	*fLog << err << dbginf << "number of OFF events and Fitted number of OFF events in signal region do not agree (within 10 %)" << endl;
	
	*fLog << "fNbg = " << fNbg << " ; fNbgFit = " << fNbgFit << endl;
    }



    // ---------------------------------------------------
    //
    // calculate the number of excess events in the signal region
    //
    fNex = fNon - fNbg * fNormFactor;
    fNexFit = fNon - fNbgFit * fNormFactor;
 

    //
    // Calculate Significance  
    //
    // Significance computed using counted number of OFF-data events in signal 
    // region (fNbg) and normalization factor (fNormFactor).
    // This is the strictly speaking the significance in Li&Ma paper...
    //
    if ( !SigmaLiMa(fNon, fNbg, fNormFactor, &fSigLiMa) )
    {
	*fLog << err << dbginf << "SigmaLiMa failed" << endl;  
	return kFALSE;
    }
    //
    // 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 ( !SigmaLiMa(fNon, fNoff, fGamma , &fSigLiMaFit) )
    {
	*fLog << err << dbginf << "SigmaLiMaFit failed"  << endl;  
	return kFALSE;
    }
    

    //
    // calculate the error of the number of excess events
    // using fitted quantities and counted quantities
    //
    fdNex = fNex / fSigLiMa;
    fdNexFit = fNexFit / fSigLiMaFit;



    //
    // Create histogram of OFF data scaled by the normalization factor
    //
    NormalizedHistOFF();


    return kTRUE;
}



// --------------------------------------------------------------------------
//
//  Create a new histogram to store the hist of OFF data scale with by the
//  normalization factor
//
void MHMyFindSignificanceONOFF::NormalizedHistOFF()
{

    //
    // We need the following factor in case that one of the histogram 
    // (ON or OFF) is rebin when calculatin the sifnigicance
    //
    const Double_t BinWidthRatioONOFF = fHistON->GetBinWidth(1)/fHistOFF->GetBinWidth(1) ;


    //
    // Clone OFF histogram and scale it
    //
    fHistOFFNormalized = (TH1*)fHistOFF->Clone();
    fHistOFFNormalized->SetName(Form("%s (Normalized)",fHistOFF->GetName()));

    fHistOFFNormalized->Scale(fNormFactor*BinWidthRatioONOFF); 

    
    //
    // Scale the fit function of OFF data
    // 
    TF1* fPolyOFFNormalized = (TF1*)fFindsigOFF.GetPoly()->Clone(); 

    fPolyOFFNormalized->SetLineColor(2);
    fPolyOFFNormalized->SetLineWidth(2);
   

    const Int_t npar = fPolyOFFNormalized->GetNpar();
    Double_t Parameter = 0.0;
    Double_t ParameterError = 0.0;  
    for(Int_t i=0; i<npar; i++)
    {
	Parameter = fNormFactor*BinWidthRatioONOFF *fPolyOFFNormalized->GetParameter(i);
	ParameterError = fNormFactor*BinWidthRatioONOFF *fPolyOFFNormalized->GetParError(i);
	
	fPolyOFFNormalized -> SetParameter(i, Parameter);
	fPolyOFFNormalized -> SetParError(i,ParameterError);
    }

    TList *funclistOFF = fHistOFFNormalized->GetListOfFunctions();
    funclistOFF->Clear();  // delete the functions cloned from fHistOFF
    funclistOFF->Add(fPolyOFFNormalized);

   
}



// --------------------------------------------------------------------------
//
//  Draw hte Alpha plot for the ON data, and on top of it, the normalized Alpha
//  plot for the OFF data.
//
void MHMyFindSignificanceONOFF::Draw(Option_t* option)
{
    
    if(!fHistON || !fHistOFF)
	return;
    
    

    //
    // If a canvas exists use it, otherwise Creates a new canvas
    //
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);
    
    //  pad->Divide(2,1);

    AppendPad("");

    pad->cd(1);
    
    //
    // Draw histograms and fit functions
    //
    fHistON->SetTitle("");
    fHistON->SetMarkerStyle(21);
    fHistON->Draw("p"); 

    fHistOFFNormalized->SetLineColor(2);
    fHistOFFNormalized->SetMarkerStyle(25);
    fHistOFFNormalized->SetMarkerColor(2);
    fHistOFFNormalized->Draw("Psame");

    TList *funclistON = fHistON->GetListOfFunctions();
    TF1* polyON = (TF1*)funclistON->FindObject("Poly");
    TF1* gaussON = (TF1*)funclistON->FindObject("PolyGauss");
    TF1* backgON = (TF1*)funclistON->FindObject("Backg");

    if(gaussON)
    {
	gaussON->SetLineColor(4);
	gaussON->SetLineWidth(4);
    }
    if(backgON)
	funclistON->Remove(backgON);
    if(polyON)
    {
	polyON->SetLineStyle(2);
	polyON->SetLineColor(2);
	polyON->SetLineWidth(2);
    } 

   
    //
    // Draw results onto the figure
    //
    DrawResults(option);


    gPad->Modified();
    gPad->Update();




 //     //
//      // Draw the difference of ON-OFF
//      //
//      pad->cd(2);
//      TH1D* h = (TH1D*)fHistON->Clone();
//      h->GetListOfFunctions()->Clear();
//      h->SetTitle("Alpha Plot ON-OFF");
//      h->Add(fHistOFFNormalized,-1.);
//      h->SetFillColor(17);
//      h->Draw("hist");

//      TLine* line = new TLine(0, 0, 90, 0);
//      line->Draw();

//      DrawResults("short");

}



// --------------------------------------------------------------------------
//
//  Draw the results on top the histograms pad
//
void MHMyFindSignificanceONOFF::DrawResults(Option_t* option)
{

    TString opt(option);

    TPaveText *pt;
    char tx[100];


    if(opt.Contains("short",TString::kIgnoreCase))
    {
	pt= new TPaveText(0.24, 0.59, 0.86, 0.87, "NDC");
	

	//  sprintf(tx, "   Nbg measured (Normalized)  = %8.1f #pm %8.1f", 
//   		 fNbg * fNormFactor,  fdNbg * fNormFactor);
//   	pt->AddText(tx);
	
//   	sprintf(tx, "   Nex (ON - OFF measured) = %8.1f #pm %8.1f",fNex,fdNex);
//   	pt->AddText(tx);

//  	Double_t ratio = fNbg>0.0 ? fNex/(fNbg*fNormFactor) : 0.0;
//   	sprintf(tx, "   Significance = %6.2f,    Nex/(Nbg*NormFactor) = %6.2f",  		fSigLiMa, ratio);
//   	pt->AddText(tx);
	//  sprintf(tx, "  Measured OFF background:");
//  	TText *t7 = pt->AddText(tx);
//  	t7->SetTextSize(0.02);
//  	t7->SetTextColor(8);
	
//  	sprintf(tx, "Results for |alpha|< %6.2f [\\circ] using fit to OFF:",
//  		fFindsigON.GetAlphasi());
//  	TText *t6 = pt->AddText(tx);
//  	t6->SetTextSize(0.03);
//  	t6->SetTextColor(8);

	//  sprintf(tx, "Noff = %6.1f", fNoff);
//  	pt->AddText(tx);

//  	sprintf(tx, "Nbg = %8.1f", 
//  		fNbgFit*fNormFactor);
//  	pt->AddText(tx);
	
//  	sprintf(tx, "Nex = %8.1f", fNexFit);
//  	pt->AddText(tx);

	
	sprintf(tx, "T_{obs} = 64 min.");
	pt->AddText(tx);

	sprintf(tx, "Zenith Angle = 7-17 #circ");
	pt->AddText(tx);
	
	sprintf(tx, "Significance = %6.2f #sigma", fSigLiMaFit);
	TText *t8 =pt->AddText(tx);
	//t8->SetTextSize(0.1);
	t8->SetTextColor(8);
    }

    else if(opt.Contains("all",TString::kIgnoreCase))
    {
	pt = new TPaveText(0.20, 0.33, 0.88, 0.90, "NDC");
	
	//
	// Fit to OFF data
	//
	sprintf(tx, "Results of polynomial fit to OFF (order %2d) :", fDegree);
	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", 
     	  fFindsigOFF.GetChisq(), fFindsigOFF.GetNdf(), fFindsigOFF.GetProb());
	pt->AddText(tx);
	
	sprintf(tx, "   Nbgtot(fit)= %8.1f #pm %8.1f", 
	        fFindsigOFF.GetNbgtotFitted(), fFindsigOFF.GetdNbgtotFitted());
	pt->AddText(tx);
	
	sprintf(tx, "   Nbgtot(meas) = %8.1f", fFindsigOFF.GetNbgtot());
	pt->AddText(tx);
	
	
	//
	// Results 
	//
	sprintf(tx, "Results for |alpha|< %6.2f [\\circ] :",
		fFindsigON.GetAlphasi());
	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);
	
	
	//  PRINT INFORMATION ABOUT MEASURED QUANTITIES  
	sprintf(tx, "  Measured OFF background:");
	TText *t7 = pt->AddText(tx);
	t7->SetTextSize(0.02);
	t7->SetTextColor(8);
	
	sprintf(tx, "   Nbg(meas)*NormFactor  = %8.1f #pm %8.1f", 
 		 fNbg*fNormFactor, fdNbg*fNormFactor);
 	pt->AddText(tx);
	
 	sprintf(tx, "   Nex (meas) = %8.1f #pm %8.1f", fNex, fdNex);
 	pt->AddText(tx);

	sprintf(tx,"   Normalization Factor (= Non/Noff) = %4.4f",fNormFactor);
	pt->AddText(tx);

	Double_t ratio = fNbg>0.0 ? fNex/(fNbg*fNormFactor) : 0.0;
 	sprintf(tx, "   Significance(mea) = %6.2f,    Nex/(Nbg*NormFactor) = %6.2f", fSigLiMa, ratio);
 	pt->AddText(tx);
	

	// PRINT INFORMATION ABOUT FITTED QUANTITIES 
	sprintf(tx, "  Fitted OFF background:");
	TText *t8 = pt->AddText(tx);
	t8->SetTextSize(0.02);
	t8->SetTextColor(8);
	
	sprintf(tx, "   Nbg(fit)*NormFactor  = %8.1f #pm %8.1f", 
		fNbgFit*fNormFactor, fdNbgFit*fNormFactor);
	pt->AddText(tx);
	
	sprintf(tx, "   Nex(fit) = %8.1f #pm %8.1f", fNexFit, fdNexFit);
	pt->AddText(tx);
	
	sprintf(tx, "   Effective Noff (i.e. fNoff) = %6.1f,   Gamma = %4.4f", 
		 fNoff, fGamma);
	pt->AddText(tx);
	
	ratio = fNbgFit>0.0 ? fNexFit/(fNbgFit*fNormFactor):0.0;
	sprintf(tx, "   Significance(fit) = %6.2f,    Nex/(Nbg*NormFactor) = %6.2f", fSigLiMaFit, ratio);
	pt->AddText(tx);
	

	//
	// Fit to ON data (polynomial+Gaus)
	//
	sprintf(tx, "Results of (polynomial+Gauss) fit to ON:");
	TText *t9 = pt->AddText(tx);
 	t9->SetTextSize(0.03);
 	t9->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]", 
 		 fFindsigON.GetSigmaGauss(), fFindsigON.GetdSigmaGauss());
 	pt->AddText(tx);
	
       	sprintf(tx, "   total no.of excess events = %8.1f #pm %8.1f", 
	      	fFindsigON.GetNexGauss(), fFindsigON.GetdNexGauss());
 	pt->AddText(tx);
    }    
    else 
	return;

	
    pt->SetFillStyle(0);
    pt->SetBorderSize(0);
    pt->SetTextAlign(12);
    
    pt->Draw();
}


// --------------------------------------------------------------------------
//
//  Set flag fRebin 
//
void MHMyFindSignificanceONOFF::SetRebin(Bool_t b)
{
  fRebin = b;

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

// --------------------------------------------------------------------------
//
//  Set flag fReduceDegree 
//
void MHMyFindSignificanceONOFF::SetReduceDegree(Bool_t b)
{
  fReduceDegree = b;

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

// ---------------------------------------------------------------------------
//
//  Given an histogram, count the number of events within a given range. 
//  Is equivaltent to TH1::Integral(binlo,binhi) but making sure that the bins
//  are fully contain in the interval (within a given tolerance).
//
Int_t MHMyFindSignificanceONOFF::CountEventsInRange(TH1* hist, Double_t alphalo, Double_t alphaup, Double_t* numevents, Double_t* numeventserror)
{    
  
    const Int_t nbins =  hist->GetNbinsX();
    const Double_t binwidth = hist->GetBinWidth(1);

    Double_t nevt = 0;
    Double_t nevtErr = 0; 

    //
    // Loop over all the bins
    // 
    for(Int_t i=1; i<=nbins; i++)
    {
	Double_t xlo = hist->GetBinLowEdge(i);
	Double_t xup = hist->GetBinLowEdge(i+1);

	// bin must be completely contained in the signal region
	if( xlo >= (alphalo-fEps) && xup <= (alphaup+fEps))
	{
	    // check all bins all of same width
	    const Double_t width = fabs(xup-xlo);
	    if (fabs(width-binwidth) > fEps)
	    {
		*fLog << err << GetName() << " alpha plot has variable binning, which is not allowed...exit."  << endl;
		return kFALSE;
	    }


	    nevt += hist->GetBinContent(i);
	    nevtErr += hist->GetBinError(i) * hist->GetBinError(i);
	}
    }

    nevtErr = TMath::Sqrt(nevtErr);
    

    // Return results
    *numevents = nevt;
    *numeventserror = nevtErr;
}


// ---------------------------------------------------------------------------
//
//  Counts the number of events in the Bkg region in the ON_Data and OFF-Data 
//  histograms. Then compute the normalization factor as:
//    NormFactor = NoffON / NoffNOFF
//
Int_t MHMyFindSignificanceONOFF::CalcNormalizationFactor(TH1* hon, TH1* hoff, Double_t alphalo, Double_t alphaup)
{
    Double_t nevton = 0;
    Double_t nevtoff = 0;
    Double_t dnevton = 0;
    Double_t dnevtoff = 0;
  
    // counts events in the region (alphalo,alphaup) for the ON-Data hist
    CountEventsInRange(hon, alphalo, alphaup,&nevton,&dnevton);

    // counts events in the region (alphalo,alphaup) for the OFF-Data hist
    CountEventsInRange(hoff, alphalo, alphaup,&nevtoff,&dnevtoff);


    // Normalization factor and its error
    fNormFactor = nevton/nevtoff;
    
    fNormFactorError = 1./nevton + 1./nevtoff;
    fNormFactorError = TMath::Sqrt(fNormFactorError);
    fNormFactorError *= fNormFactor; // <--------------


    *fLog << "non = " << nevton << " noff = " << nevtoff << endl;
    *fLog << "Normalization factor = " << fNormFactor << " +- " 
	  << fNormFactorError << endl;

    return kTRUE;
}



// --------------------------------------------------------------------------
//
//  SigmaLiMa
//
//  calculates the significance according to Li & Ma
//  ApJ 272 (1983) 317
//
Bool_t MHMyFindSignificanceONOFF::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 << "MHMyFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = "
  //      << non << ",  " << noff << ",  " << gamma << ",  " << *siglima << endl;

  return kTRUE;
}


// --------------------------------------------------------------------------
//
//  Sigma Vs Alpha
//
//
Bool_t MHMyFindSignificanceONOFF::SigmaVsAlpha(TH1 *histON, TH1 *histOFF,  Double_t alphamin, Double_t alphamax, Int_t degree, Bool_t draw)
{


    //--------------------------------------------
    // Create histogram
    
    //Int_t    nsteps    =  15;
    Int_t    nsteps    = histON->FindBin(40); 
  
    TH1D* fSigVsAlpha = new TH1D("SigVsAlpha","Sigma vs Alpha", nsteps, 0.0, 30);
     MH::SetBinning(fSigVsAlpha, histON->GetXaxis());
     fSigVsAlpha->SetXTitle("upper edge of signal region in |alpha|  [\\circ]");
    fSigVsAlpha->SetYTitle("Significance of gamma signal");


    TGraph* grSigVsAlpha = new TGraph;
    TGraph* grNexVsAlpha = new TGraph;
 
  


    //--------------------------------------------
    // loop over different signal regions
    for (Int_t i=1; i<=nsteps; i++)
    {
	//Double_t alphasig = fSigVsAlpha->GetBinCenter(i);
	//Double_t alphasig = fSigVsAlpha->GetBinLowEdge(i+1);
	Double_t alphasig = histON->GetBinLowEdge(i+1);

	
	gLog.SetNullOutput(1);
	//FindSigmaONOFF(histON, histOFF, alphasig, alphamin, alphamax, degree);
	 MHMyFindSignificanceONOFF findsig;
	 findsig.SetRebin(fRebin);
	 findsig.SetReduceDegree(fReduceDegree);
  
	 findsig.FindSigmaONOFF( histON, histOFF, alphasig, alphamin, alphamax, degree);
	 const Double_t sig = findsig.GetSigLiMa();
	 const Double_t Nex = findsig.GetNex();

	 fSigVsAlpha->SetBinContent(i, sig);

	 if( sig>fBestSigma )
	 {
	     fBestSigma = sig;
	     fBestAlphaCut = alphasig;
	     fBestExcess = Nex;
	 }

	 grSigVsAlpha->SetPoint(i-1,alphasig,sig);
	 grNexVsAlpha->SetPoint(i-1,alphasig,Nex);

	 cout << "****************** " << i<< " " << alphasig << " " << sig<< " " << Nex << "*************"<<  endl;
	 gLog.SetNullOutput(0);
    }


// fSigVsAlpha->DrawCopy();


    if (!draw)
	return kTRUE;

    // -----------------------------------------------------------------
    //
    // Draw
    // 
    TCanvas *c = new TCanvas("SigVsAlpha", "Sigma vs Alpha", 600, 600);
    c->SetGridx();
    c->SetGridy();
    
    Double_t margin = 0.12; // give more space for drawing the axis titles
    c->SetLeftMargin(margin);
    c->SetRightMargin(margin);


    //
    // Draw Sig Vs Alpha
    //
    grSigVsAlpha->SetMarkerStyle(8);
    //grSigVsAlpha->SetMarkerColor(2);
    //grSigVsAlpha->SetLineColor(2);
    grSigVsAlpha->Draw("apl");
    grSigVsAlpha->GetHistogram()->SetXTitle("Alpha Cut [\\circ]");
    grSigVsAlpha->GetHistogram()->SetYTitle("Significance");
    grSigVsAlpha->GetHistogram()->GetYaxis()->SetTitleOffset(1.3);


    //
    // Draw second graph on a transparent pad
    //
    TPad *overlay = new TPad("overlay","",0,0,1,1);
    overlay->SetFillStyle(4000);
    overlay->SetFillColor(0);
    overlay->SetFrameFillStyle(4000);
    overlay->SetLeftMargin(margin);
    overlay->SetRightMargin(margin);  
    overlay->Draw();
    overlay->cd();
  
    grNexVsAlpha->SetMarkerStyle(8);
    grNexVsAlpha->SetMarkerColor(kRed);
    grNexVsAlpha->SetLineColor(kRed);
    grNexVsAlpha->Draw("apl"); 

    // Trick to not draw the graph Y-Axis
    grNexVsAlpha->GetHistogram()->GetYaxis()->SetNdivisions(0);
    
     
    //
    // Draw an axis on the right side
    //
    c->Update(); //<-- Imprescindible for updating the gPad coordinates

    TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
			    gPad->GetUxmax(), gPad->GetUymax(),  
			    gPad->GetUymin(),gPad->GetUymax() ,510,"+L");
    axis->SetTitle("Excess events");
    axis->SetTitleOffset(1.3);
    //axis->SetLineColor(kRed);
    axis->SetLabelColor(kRed);
    axis->SetTextColor(kRed);
    axis->Draw();
   

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

    return kTRUE;
}
