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

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

#include "MHExcessEnergyTheta.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"

ClassImp(MHExcessEnergyTheta);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. It sets name and title of the histogram. 
//
MHExcessEnergyTheta::MHExcessEnergyTheta(const char *name, const char *title)
  :  fHist("","",10,0,1, 10,0,1)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHExcessEnergyTheta";
    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 MHExcessEnergyTheta::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 MHExcessEnergyTheta::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 MHExcessEnergyTheta::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 MHExcessEnergyTheta::Calc(MHAlphaEnergyTheta* hAlpha)
{  
    TH3D* aet = (TH3D*)hAlpha->GetHist();
  
    const TAxis* axisEnergy = aet->GetXaxis();
    const TAxis* axisTheta  = aet->GetYaxis();
    const Int_t energyBins = aet->GetXaxis()->GetNbins();
    const Int_t thetaBins =  aet->GetYaxis()->GetNbins();;
  
    MH::SetBinning(&fHist, axisEnergy, axisTheta);
 
    for(Int_t iy=1; iy<=thetaBins; iy++)   
    {
	TCanvas* c= new TCanvas(Form("zbin %d",iy),Form("zbin %d",iy));
	c->Divide(4,3);

	for(Int_t ix=1; ix<=energyBins; ix++)
	{ const Double_t e = aet->GetXaxis()->GetBinCenter(ix); 
	    const Double_t e1 = aet->GetXaxis()->GetBinLowEdge(ix);     
	    const Double_t e2 =	aet->GetXaxis()->GetBinLowEdge(ix+1);

	    //TH1* halpha = aet->ProjectionZ("AlphaTemp",ix,ix,iy,iy);
	    TH1* halpha = aet->ProjectionZ("AlphaTemp",ix,ix);
	    halpha->SetTitle(Form("Energy Bin = (%.0f, %.0f) GeV",e1,e2));
	    halpha->SetDirectory(NULL);

	    MHFindSignificance findsig;
	    //findsig.SetRebin(kFALSE);
	    double alphasig=20;
	    
	    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 = 20;
		break;
	    case 7:
		alphasig = 15;
		break;
	    case 8:
		alphasig = 10;
		break;
	    case 9:
		alphasig = 5;
		break;
	    case 10:
		alphasig = 5;
		break;
	    }


	    fLog->SetNullOutput(kTRUE);
	    
	    Bool_t rc = findsig.FindSigma(halpha, 30,90, 2, alphasig, 0,1,0);
	    Double_t SigmaGauss = findsig.GetSigmaGauss(); 

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

	    fLog->SetNullOutput(kFALSE);

	    if(ix<=20)
	    {
		c->cd(ix);
		findsig.DrawFit("brief");
	    }

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



// --------------------------------------------------------------------------
//
// Draw the histogram
// 
void MHExcessEnergyTheta::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");
}
