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

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHExcessEnergyThetaONOFF                                                      //
//                                                                          //
//  3D-histogram in alpha vs. E-est and Theta                                //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MHExcessEnergyThetaONOFF.h"

#include <TCanvas.h>
#include <THStack.h>
#include <TLegend.h>
#include <TStyle.h>

#include "MMcEvt.hxx"
#include "MHillasSrc.h"
#include "MEnergyEst.h"
#include "MPointingPos.h"
#include "MRawRunHeader.h"
#include "MHAlphaEnergyTheta.h"
#include "MHFindSignificance.h"

#include "MBinning.h"
#include "MParList.h"

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

#include <TGraphErrors.h>

#include "MHMyFindSignificanceONOFF.h"


ClassImp(MHExcessEnergyThetaONOFF);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. It sets name and title of the histogram. 
//
MHExcessEnergyThetaONOFF::MHExcessEnergyThetaONOFF(const char *name, const char *title)
  :  fHist("","",10,0,1, 10,0,1)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHExcessEnergyThetaONOFF";
    fTitle = title ? title : "# Excess events vs. E and Theta";


    fHist.SetDirectory(NULL);

    fHist.SetTitle("# Excess events vs. E and Theta");
    fHist.SetXTitle("E [GeV]");
    fHist.SetYTitle("\\Theta [\\circ]");
    fHist.SetZTitle("# Excess events");



}

// --------------------------------------------------------------------------
//
// Set binnings and prepare filling of the histogram
// 
Bool_t MHExcessEnergyThetaONOFF::SetupFill(const MParList *plist)
{

 //   fEnergy = (MEnergyEst*)plist->FindObject("MEnergyEst");
//    if (!fEnergy)
//    {
//        *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
//        return kFALSE;
//    }
   

//    //
//    // Binning
//    //
//    MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE");
//    MBinning* binsalphaflux  = (MBinning*)plist->FindObject("BinningAlphaFlux");
//    MBinning* binstheta  = (MBinning*)plist->FindObject("BinningTheta");
//    if (!binsenergy || !binsalphaflux || !binstheta)
//    {
//        *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl;
//        return kFALSE;      
//    }

// //   SetBinning(&fHist, binsalphaflux, binsenergy, binstheta);
//    SetBinning(&fHist, binsenergy, binstheta, binsalphaflux);

//    fHist.Sumw2(); 


   return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Look in the parlist for MMcEvt or MPointingPos depending on the run type. 
//
Bool_t MHExcessEnergyThetaONOFF::ReInit(MParList *pList)
{

//     // This must be done in ReInit because in PreProcess the
//     // runheaders are not available
//     const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
//     if (!runheader)
//     {
//         *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
// 	*fLog << warn << dbginf << "Warning - Assuming data type file...searching for MPointingPos..." << endl;
//     }
   

//     if (runheader && runheader->IsMonteCarloRun())
//     {	
// 	fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
// 	if (!fMcEvt)
// 	{
// 	    *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
// 	    return kFALSE;
// 	}
//     }
//     else
//     {
// 	fPointingPos = (MPointingPos*)pList->FindObject("MPointingPos");
// 	if (!fPointingPos)
// 	{
// 	    *fLog << err << dbginf << "PointingPos not found... aborting." << endl;
// 	    return kFALSE;
// 	}
//     }
  
    return kTRUE;
}



// --------------------------------------------------------------------------
//
// Fill the histogram
// 
Bool_t MHExcessEnergyThetaONOFF::Fill(const MParContainer *par, const Stat_t w)
{
 //    MHillasSrc &hil = *(MHillasSrc*)par;

//   //   fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(),
// //                fMcEvt->GetTelescopeTheta()*kRad2Deg, w);
 

//     const Double_t Zd = (fMcEvt) ? fMcEvt->GetTelescopeTheta()*kRad2Deg : fPointingPos->GetZd() ;

//     fHist.Fill(fEnergy->GetEnergy(), Zd, hil.GetAlpha(), w);
	
//     //  *fLog<< Zd << " " << fEnergy->GetEnergy() << " " <<  hil.GetAlpha() << endl;

    return kTRUE;
}



// --------------------------------------------------------------------------
//
//  Get the Alpha distribution for each bin in energy and theta,
//  fit it with MHFindSignificance to get the number of Excess events
//  and it error.
// 
void MHExcessEnergyThetaONOFF::Calc(MHAlphaEnergyTheta* hAlphaON, MHAlphaEnergyTheta* hAlphaOFF)
{  
    *fLog << endl;
    fLog->Separator("Excess Events Calculation");
    

    TH3D* aetON  = (TH3D*)hAlphaON->GetHist();
    TH3D* aetOFF = (TH3D*)hAlphaOFF->GetHist();
  
    const TAxis* axisEnergy = aetON->GetXaxis();
    const TAxis* axisTheta  = aetON->GetYaxis();
    const Int_t energyBins = aetON->GetXaxis()->GetNbins();
    const Int_t thetaBins =  aetON->GetYaxis()->GetNbins();;
  
    MH::SetBinning(&fHist, axisEnergy, axisTheta);
 


    //
    // Loop over Theta bins
    //
    for(Int_t iy=1; iy<=thetaBins; iy++)   
    {
	TCanvas* c= new TCanvas(Form("zbin %d",iy),Form("zbin %d",iy));
	c->Divide(4,3);

	//
	// Loop over Energy bins
	//
	for(Int_t ix=1; ix<=energyBins; ix++)
	{ 
	    const Double_t e = aetON->GetXaxis()->GetBinCenter(ix); 
	    const Double_t e1 = aetON->GetXaxis()->GetBinLowEdge(ix);     
	    const Double_t e2 =	aetON->GetXaxis()->GetBinLowEdge(ix+1);

	    //
	    // Get Alpha plot for that Theta & Energy bin
	    //
	    TH1* halphaON = aetON->ProjectionZ("AlphaTempON",ix,ix,iy,iy);
	    halphaON->SetTitle(Form("Energy Bin = (%.0f, %.0f) GeV",e1,e2));
	    halphaON->SetDirectory(NULL);

	    TH1* halphaOFF = aetOFF->ProjectionZ("AlphaTempOFF",ix,ix,iy,iy);
	    halphaOFF->SetTitle(Form("Energy Bin = (%.0f, %.0f) GeV",e1,e2));
	    halphaOFF->SetDirectory(NULL);

	    TCanvas *c = new TCanvas;
	    halphaOFF->Draw();


	    //
	    // Fit Alpha plot to get Nex
	    //
	    MHMyFindSignificanceONOFF findsig;
	    findsig.SetRebin(kFALSE); 
	    //findsig.SetRebin(kTRUE);
	    double alphasig=25;
	    
	//      if(ix<5) 
//  		alphasig = 20;
//  	    else
//  		alphasig = 20;
	     
//  	  switch (ix)
//  	    {
//  	    case 1:
//  		alphasig = 20;
//  		break;
//  	    case 2:
//  		alphasig = 20;
//  		break;
//  	    case 3:
//  		alphasig = 20;
//  		break;
//  	    case 4:
//  		alphasig = 20;
//  		break;
//  	    case 5:
//  		alphasig = 20;
//  		break;
//  	    case 6:
//  		alphasig = 10;
//  		break;
//  	    case 7:
//  		alphasig = 10;
//  		break;
//  	    case 8:
//  		alphasig = 10;
//  		break;
//  	    case 9:
//  		alphasig = 10;
//  		break;
//  	    case 10:
		
//  		alphasig = 10;
//  		break;
//  	    }

  const Double_t AlphaMinON = 30;
    const Double_t AlphaMaxON = 90;
    const Double_t AlphaMinOFF = 0;
    const Double_t AlphaMaxOFF = AlphaMaxON;
    const Double_t Degree = 4;

	  fLog->SetNullOutput(kTRUE);
	    

	    Bool_t rc =  findsig.FindSigmaONOFF(halphaON, halphaOFF, alphasig, AlphaMinON,AlphaMaxON,Degree);
 
	    //  Double_t SigmaGauss = findsig.GetSigmaGauss(); 

	  //     alphasig = SigmaGauss*3;	 
//    	    rc = findsig.FindSigma(halpha, 25,90, degree, alphasig, 0,1,0);

	    fLog->SetNullOutput(kFALSE);

	

	    if(ix<=20)
	    {
		c->cd(ix);
		gPad->SetGridx();
		gPad->SetGridy();
		
		findsig.Draw("all");
	    }

	    Double_t Nex = 0;  
	    Double_t dNex = 0;
	    Double_t Sig = 0;

	    if(!rc)
	    {
		cout << "ERROR: Fit not possible. " << endl;
	    }
	    // else
	    {
		//SigmaGauss = findsig.GetSigmaGauss(); 
		Nex = findsig.GetNex();  
		dNex = findsig.GetdNex();    
		Sig = findsig.GetSigLiMa();
	    }
	    // Avoid to have a negative number of excess events, which 
            // will give a negative flux which doesn't make sense
	    if(Nex<0)
	    {
		Nex = 0;
		dNex = 0;
		Sig = 0;
	    }

	    fHist.SetBinContent(ix,iy, Nex);
	    fHist.SetBinError(ix,iy, dNex);


	    // Print 
	    *fLog<< endl 
		 << " Theta Bin = "  << iy << endl
		 << "  Energy bin width E = "  << ix << endl
		 << "  Alphasig = " <<alphasig << endl
		 << "    N = " <<  Nex    <<" +- " <<  dNex    << endl
		 << "  Significane = " <<  Sig  << endl;
	}
    }


    *fLog << "Total Number of Excess evts = " << fHist.Integral() << endl;
    

}



// --------------------------------------------------------------------------
//
// Draw the histogram
// 
void MHExcessEnergyThetaONOFF::Draw(Option_t *opt)
{
    TCanvas *c1 = new TCanvas("# Excess events vs. E and Theta","# Excess events vs. E and Theta");
    c1->SetLogx();
    c1->SetLogz();

    fHist.SetStats(0);
    fHist.Draw("lego");

    TCanvas *c2 = new TCanvas;
    c2->SetLogx();

    TH1* h=   fHist.ProjectionX();
    h->SetStats(0);
    h->SetTitle("No. Excess events Vs. Energy");
    h->SetXTitle("E [GeV]");
    h->SetYTitle("# Excess events");
  
    h->Draw("e1");


    /////////////
    TCanvas *c3 = new TCanvas;
   

    Double_t res[10] ={0.224281,0.30083, 0.397355, 0.455888, 0.482985, 0.479416, 0.441429, 0.424654, 0.366658, 0.287384};

    TGraphErrors* graphIntegral = new TGraphErrors(h);
    
    for(int i=0; i<graphIntegral->GetN(); i++)
    {
	Double_t errorX = graphIntegral->GetErrorX(i);
	Double_t errorY = graphIntegral->GetErrorY(i);
	errorX *= res[i];

       	graphIntegral->SetPointError(i,errorX,errorY);
    }

  
    c3->SetLogx();
    c3->SetLogy();
    c3->SetGridx();
    c3->SetGridy();

   graphIntegral->SetMarkerStyle(8);
     graphIntegral->GetHistogram() ->SetXTitle("E [GeV]"); 
    graphIntegral->GetHistogram() ->SetYTitle("# Excess events");
      graphIntegral->Draw("Ap");

}