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

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

#include "MHFlux.h"

#include <TCanvas.h>
#include <THStack.h>
#include <TLegend.h>
#include <TStyle.h>
#include <TAxis.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TPaveText.h>

#include "MHillasSrc.h"
#include "MEnergyEst.h"
#include "MPointingPos.h"
#include "MRawRunHeader.h"

#include "MHExcessEnergyTheta.h"
#include "MHMcCollectionArea.h"
#include "MHEffectiveOnTime.h"

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

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

ClassImp(MHFlux);

using namespace std;

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


    fHist.SetDirectory(NULL); 
    fHist.SetName("Flux vs. E and Theta");
    fHist.SetTitle("Flux vs. E and Theta");
    fHist.SetXTitle("E [GeV]");
    fHist.SetYTitle("\\Theta [\\circ]");
    fHist.SetZTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]");

    fAverageFlux.SetDirectory(NULL);
    fAverageFlux.SetName("Average Flux");
    fAverageFlux.SetTitle("Average Flux");
    fAverageFlux.SetXTitle("E [GeV]"); 
    fAverageFlux.SetYTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]");
}

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

 
   return kTRUE;
}

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



// --------------------------------------------------------------------------
//
// Fill the histogram
// 
Bool_t MHFlux::Fill(const MParContainer *par, const Stat_t w)
{

    return kTRUE;
}


// --------------------------------------------------------------------------
//
//  Calc
// 
void MHFlux::Calc(MHExcessEnergyTheta* hExcess, MHMcCollectionArea* hColArea, MHEffectiveOnTime* hEffTime)
{
    const TH2D* hex  = hExcess->GetHist();  
    const TH1D* hca  = hColArea->GetHist(); 
    const TH1D* heot = &hEffTime->GetHEffOnTheta();

    TAxis* axisEnergy = hex->GetXaxis();
    const TAxis* axisTheta  = hex->GetYaxis();
    const Int_t energyBins = hex->GetXaxis()->GetNbins();
    const Int_t thetaBins =  hex->GetYaxis()->GetNbins();;

    MH::SetBinning(&fHist,axisEnergy, axisTheta);
 
 
    //
    // Calculate flux for each energy and theta
    //
    for (Int_t iy=1; iy<=thetaBins; iy++) // loop on theta
    { 
        // Get Effective Time [sec] and its error
	Double_t t = heot->GetBinContent(iy);
	const Double_t dt = heot->GetBinError(iy); 

	if (t < 1e-3) 
	    t = 0.0;

	for (Int_t ix=1; ix<=energyBins; ix++) 
	{   
	    const Double_t  n = hex->GetBinContent(ix,iy);
	    const Double_t dn = hex->GetBinError(ix,iy);
	    
	    // Get AreaEff and its error
	    Double_t energy = axisEnergy->GetBinCenter(ix);	    
	    Int_t bin = hca->GetXaxis()->FindBin(energy);	

	    // Get NumberExcessEvents and its error
	    const Double_t  a = hca->GetBinContent(bin)*1e4; //cm^2
	    const Double_t da = hca->GetBinError(bin) *1e4; //cm^2
 
	    // energy bin width in TeV
	    const Double_t en = axisEnergy->GetBinWidth(ix)*1e-3; //TeV

	    // 
	    // Check that we can calculate the flux for the current bin
	    //
	    if (t==0) 
		cout << "No_Ton ";
	    if (a==0) 
		cout << "No_Aeff ";
	    if (n==0) 
		cout << "No_Events ";
	    if ((t == 0) || (a == 0) || (n == 0)) {
		cout << endl;
		continue;
	    }
 	    
	    //
	    // Flux calculation and its error
	    //
	    const Double_t flux = n/(en*t*a);
      
	    // error  propagation formula
	    const Double_t errN = dn/(en*a*t);
	    const Double_t errA = da * n/(en*t*a*a);
	    const Double_t errT = dt * n/(en*a*t*t);      
	    const Double_t error = sqrt(errN*errN + errA*errA + errT*errT); 

	    cout << dn << " " << en << " " << a << " " << t << endl;

	    fHist.SetBinContent(ix,iy,flux);
	    fHist.SetBinError(ix,iy,error);

	}  //energy    
    } //theta
    fHist.Print("all");
}



// --------------------------------------------------------------------------
//
// Draw the histogram
// 
void MHFlux::Draw(Option_t *opt)
{

    // --------------------------------------------------------------------
    //
    // Draw lego plot of Flux vs. Energy and Theta
    //
    TCanvas *c1 = new TCanvas("Flux vs. E and Theta","Flux vs. E and Theta");
    c1->SetLogx();
    c1->SetLogz(); 
    fHist.SetStats(0);
    fHist.Draw("lego");


    // --------------------------------------------------------------------
    //
    // Draw the Flux for each Theta bin
    //
    TCanvas *c2 = new TCanvas("Fluxes for each Theta bin","Fluxes for each Theta bin");
    c2->SetLogx();
    c2->SetLogy();
    c2->SetGridx();
    c2->SetGridy();


    THStack* hs = new THStack("Fluxes for each Theta bin","Fluxes for each Theta bin");

    TLegend * leg = new TLegend(0.73,0.65,0.89,0.89);
   
    TAxis* yaxis = fHist.GetYaxis();
    const Int_t nbiny = fHist.GetYaxis()->GetNbins();
    
    for(Int_t iy=1; iy<=nbiny; iy++)
    {

	TH1D* h1 = fHist.ProjectionX(Form("%d",iy),iy,iy,"e"); //<----- Option e is very important, otherwise the errors are not copied

	if(h1->GetEntries()==0)
	    continue;

	h1->SetLineColor(iy);
	hs->Add(h1,"e1");
	leg->AddEntry(h1,Form("\\theta = %.0f",yaxis->GetBinCenter(iy)),"l");

// 	TCanvas *c = new TCanvas();
// 	c->SetLogx();
// 	c->SetLogy();
// 	h1->DrawCopy("e"); 
// 	h1->Print("all");
    }



    // --------------------------------------------------------------------
    //
    // Calculate and Draw the Flux average on Theta bins
    //
    fAverageFlux.SetStats(0);

    MH::SetBinning(&fAverageFlux,fHist.GetXaxis());

    for(int ix=1; ix<=fHist.GetXaxis()->GetNbins(); ix++) // energy 
    {
	Double_t sumw = 0;
	Double_t sumcontents = 0;
	Double_t sumerrors = 0;

	for(int iy=1; iy<=fHist.GetYaxis()->GetNbins(); iy++) // theta
	{
	    Double_t weight = fHist.GetYaxis()->GetBinWidth(iy);
	    sumw += weight;
	    sumcontents += fHist.GetBinContent(ix,iy)*weight;
	    sumerrors += fHist.GetBinError(ix,iy)*weight;
	}

	fAverageFlux.SetBinContent(ix,sumcontents/sumw);
	fAverageFlux.SetBinError(ix,sumerrors/sumw);
    }

  //   for(int ix=1; ix<=fHist.GetXaxis()->GetNbins(); ix++) // energy 
//     {
// 	Double_t sumw = 0;
// 	Double_t sumcontents = 0;

// 	cout << " energy bin "<<endl; 
// 	for(int iy=1; iy<=fHist.GetYaxis()->GetNbins(); iy++) // theta
// 	{

// 	    Double_t bincontent = fHist.GetBinContent(ix,iy);
// 	    Double_t binerror = fHist.GetBinError(ix,iy);
	    
// 	    if( bincontent == 0 || binerror == 0 )
// 		continue;
// 	    cout << binerror << endl;

// 	    Double_t weight = 1/(binerror*binerror);
// 	    sumw += weight;
// 	    sumcontents += bincontent*weight;
	
// 	    cout << " theta bin "  << fHist.GetBinContent(ix,iy)<< "  " <<weight
// 		 << endl;
// 	}
// 	cout << "*****************" << sumcontents << "   "<< sumw << endl;

// 	if(sumcontents == 0 || sumw == 0 )
// 		continue;

// 	fAverageFlux.SetBinContent(ix,sumcontents/sumw);
// 	fAverageFlux.SetBinError(ix,TMath::Sqrt(1/sumw));
//     }


   

    fAverageFlux.SetMarkerStyle(8);
    fAverageFlux.SetLineColor(6);
    hs->Add(&fAverageFlux,"pe1");
    leg->AddEntry(&fAverageFlux,"Average on Theta","l");
  
    
    c2->cd();
    hs->Draw("nostack");
    leg->Draw();

 
    TCanvas *c3 = new TCanvas("Average Flux","Average Flux");
    c3->SetLogx();
    c3->SetLogy();
    c3->SetGridx();
    c3->SetGridy();
    fAverageFlux.Draw();
    fAverageFlux.Print("all");

    //
    // Fix the Average Flux to a power law
    //
    TF1* fluxfit = new TF1("f1","[0]*pow(x,-[1])",90,1500);
    fluxfit->SetParNames("f0","a");
    fluxfit->SetParameter(0,5.10986e-05);
    fluxfit->SetParameter(1,2.4);
    fluxfit->SetTitle("Flux fit");
    fluxfit->SetLineColor(27);
    fluxfit->SetLineWidth(3);

    fAverageFlux.Fit("f1","R");  
   
  
    //
    // Draw the Crab spectrum measured by HEGRA between 500 GeV and 80 TeV
    //
    TF1* CrabFlux = new TF1("CrabFlux","[0]*pow(x/1000.,-[1])",350,2000);
    CrabFlux->SetParameter(0,2.83e-11);
    CrabFlux->SetParameter(1,2.62);
    CrabFlux->SetLineStyle(2); 
    CrabFlux->SetLineColor(4);
    CrabFlux->Draw("same");
   
    //
    // Draw formula
    //
    TPaveText* func = new TPaveText(0.16, 0.22, 0.67, 0.28,"NDC");
    func->AddText(Form("#frac{dF}{dE} = %.2e * E^{-%.2f} [#frac{ph}{cm^{2} s TeV}]",fluxfit->GetParameter(0),fluxfit->GetParameter(1)));
    func->SetFillStyle(0);
    func->SetBorderSize(0);
    func->Draw();  
    

//     //
//     // Draw "Preliminary"
//     //
//     TPaveText* lab = new TPaveText(0.33, 0.83, 0.68, 0.89,"NDC");
//     lab->AddText("preliminary");
//     lab->SetTextColor(2);
//     lab->SetFillStyle(0);
//     lab->SetBorderSize(0);
//     lab->Draw();  
     


    // ---------------------------------------------------------------------
    //
    // Integral flux
    // 
    TH1D *hIntegral = (TH1D*)fAverageFlux.Clone();
    hIntegral->GetListOfFunctions()->Clear();

    Int_t nbinsx = fAverageFlux.GetNbinsX();
    
    for(int i=1; i<=nbinsx; i++)
    {
	cout <<"Integral Flux: Binwidth:" << hIntegral->GetBinWidth(i) << endl;
	
	hIntegral->SetBinContent(i,hIntegral->GetBinContent(i)*hIntegral->GetBinWidth(i)*1e-3);
	hIntegral->SetBinError(i,hIntegral->GetBinError(i)*hIntegral->GetBinWidth(i)*1e-3);
    }


    for(int i=nbinsx-1; i>=1; i--)
    {
	Double_t integralsofar = hIntegral->GetBinContent(i+1);
	Double_t current = hIntegral->GetBinContent(i);
	
	Double_t currentE = hIntegral->GetBinError(i);
	Double_t Esofar = hIntegral->GetBinError(i+1);

	hIntegral->SetBinContent(i,(current+integralsofar));
	hIntegral->SetBinError(i,TMath::Sqrt(currentE*currentE+Esofar*Esofar));
    }

    
    hIntegral->SetTitle("Integral Flux"); 
    hIntegral->SetXTitle("E [GeV]"); 
    hIntegral->SetYTitle("Integral Flux [s^{-1} cm^{-2}]");
        
    
    TCanvas *c20 = new TCanvas();
    c20->SetLogx();
    c20->SetLogy();
    c20->SetGridx();
    c20->SetGridy();

    hIntegral->Draw();


    

    // --------------------------------------------------------------------
    //
    //  E^2 * Flux
    //
    TH1D *hEscaledFlux = (TH1D*)fAverageFlux.Clone();
   hEscaledFlux->GetListOfFunctions()->Clear();
   
    nbinsx = hEscaledFlux->GetNbinsX();
    
    for(int i=1; i<=nbinsx; i++)
    {
	
	Double_t energy = hEscaledFlux->GetBinLowEdge(i)*1e-3; // TeV
	Double_t Flux = hEscaledFlux->GetBinContent(i);
	Double_t dFlux = hEscaledFlux->GetBinError(i);

	hEscaledFlux->SetBinContent(i,energy*energy*Flux);
	hEscaledFlux->SetBinError(i,energy*energy*dFlux);
    }

    TCanvas *c40 = new TCanvas();
    c40->SetLogx();
    c40->SetLogy();
    c40->SetGridx();
    c40->SetGridy();

    hEscaledFlux->SetTitle("Escaled Flux"); 
    hEscaledFlux ->SetXTitle("E [GeV]"); 
    hEscaledFlux ->SetYTitle("E^{2}*Flux [TeV s^{-1} cm^{-2}]");
    
    hEscaledFlux->Draw();


//     // -----------------------------------------------------------------
//     //
//     // Graph move a 30 %
//     //
//     TCanvas *c4 = new TCanvas();
//     c4->SetLogx();
//     c4->SetLogy();
//     c4->SetGridx();
//     c4->SetGridy();

//     Int_t nbins = fAverageFlux.GetNbinsX();
// 	TArrayD x(nbins),y(nbins),dx(nbins),dy(nbins);

//     for(int i=1; i<=nbins; i++)
//     {
// 	 x[i-1] = fAverageFlux.GetXaxis()->GetBinCenter(i)*.7;
// 	 y[i-1] = fAverageFlux.GetBinContent(i);
// 	dx[i-1] = fAverageFlux.GetXaxis()->GetBinWidth(i)*0.62;
// 	dy[i-1] = fAverageFlux.GetBinError(i);
//     }

//     TGraphErrors* gr = new TGraphErrors(fAverageFlux.GetNbinsX(), x.GetArray(), y.GetArray(), dx.GetArray(), dy.GetArray());
   
//     gr->SetMarkerStyle(8);
//     gr->Draw("Ap");
//     gr->Print("all");  

//   //
//     TF1* fluxfit2 = new TF1("f2","[0]*pow(x,-[1])",70,2500);
//     fluxfit2->SetParNames("f0","a");
//     fluxfit2->SetParameter(0,5.10986e-05);
//     fluxfit2->SetParameter(1,2.4);
//     fluxfit2->SetTitle("Flux fit");
//     fluxfit2->SetLineColor(27);
//     fluxfit2->SetLineWidth(3);


//     gr->Fit("f2","R");  
   
//     gr->SetTitle("");
//     gr->SetMaximum(1e-3);
//     gr->SetMinimum(1e-12);
 
//     TLegend* leg2 = new TLegend(0.67,0.72,0.89,0.89);
   
//     leg2->AddEntry(gr,"MAGIC Sept. 2004","p");
 
//     gr->SetMarkerStyle(8);
//     //  gr->SetLineColor(6);

// //     //
// //     // Draw the Crab spectrum measured by HEGRA between 500 GeV and 80 TeV
// //     //
// //     TF1* CrabFlux = new TF1("CrabFlux","[0]*pow(x/1000.,-[1])",350,2000);
// //     CrabFlux->SetParameter(0,2.83e-11);
// //     CrabFlux->SetParameter(1,2.62);
// //     CrabFlux->SetLineStyle(2); 
// //     CrabFlux->SetLineColor(4);
// //     CrabFlux->Draw("same");

//     leg2->AddEntry(CrabFlux,"HEGRA ApJ 614","l");
//     leg2->Draw();
//     lab->Draw(); 
 
//     TPaveText* func2 = new TPaveText(0.16, 0.22, 0.67, 0.28,"NDC");
//     func2->AddText(Form("#frac{dF}{dE} = %.2e * E^{-%.2f} [#frac{ph}{cm^{2} s TeV}]",fluxfit2->GetParameter(0),fluxfit2->GetParameter(1)));
//     func2->SetFillStyle(0);
//     func2->SetBorderSize(0);
//     func2->Draw();  


//     gr->GetHistogram()->SetXTitle("E [GeV]"); 
//     gr->GetHistogram()->SetYTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]");

//     TH1F* h = gr->GetHistogram();
//     TCanvas *c12 = new TCanvas();
//     h->Draw();
//     cout << "Integral flux = "<< h->Integral("width") << endl;
}