/*************************************************************************
* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

//  -------------------------------------------------------
//  pt spectra extrapolation in the 0. - 0.2 region using
//  Boltzmann-Gibbs Blast Wave model or Tsallis Blast
//  Wave model for azimuthal isotropic  expansion in
//  highly central collisions analysis
//  author: Cristian Andrei
//          acristian@niham.nipne.ro
//  ----------------------------------------------------------


#include "TH1D.h"
#include "TFile.h"
#include "TList.h"
#include "TCanvas.h"
#include "TLegend.h"
// #include "TString.h"
// #include "TPaveStats.h"
#include "TMath.h"
#include "TStopwatch.h"
#include "TStyle.h"
#include "TF1.h"

#include "TVirtualFitter.h"
#include "Math/WrappedFunction.h"
#include "Math/Integrator.h"
#include "Math/IntegratorMultiDim.h"

#include "AliCFContainer.h"
#include "AliCFDataGrid.h"
#include "AliCFEffGrid.h"

#include "AliAnalysisCentralExtrapolate.h"


ClassImp(AliAnalysisCentralExtrapolate)

//________________________________________________________________________
AliAnalysisCentralExtrapolate::AliAnalysisCentralExtrapolate(const char *name) 
  : TObject()
  ,fPartType()
  ,fInputList(0)
  ,fResultsList(0x0)

{
// Constructor

	Printf("Running %s!\n", name);

	fInputList = new TList();
	fResultsList = new TList();

}

AliAnalysisCentralExtrapolate::~AliAnalysisCentralExtrapolate() {

	if(fInputList) delete fInputList;
	if(fResultsList) delete fResultsList;

}

//************************************************************************************
//____________________________________________________________________________________
// using global variables to avoid the limitations of ROOT::Math::WrappedMultiFunction<>
static Double_t mass = 0.;
static Double_t pt = 0.;
static Double_t T = 0.;
static Double_t betaS = 0.;
static Double_t q = 0.;
static Double_t mt = 0.;

void AliAnalysisCentralExtrapolate::ApplyEff(){
// applies the efficiency map to the pt spectra

    TH1::SetDefaultSumw2();
//     Int_t stepGen = 0;
    Int_t stepRec = 1;

	AliCFContainer *data = 0;
	
	
	TFile *file1 = new TFile("$ALICE_ROOT/PWG2/data/AliAnalysisCentralEfficiency.root", "read");
    
    AliCFEffGrid *eff = 0;
    
    if(fPartType.Contains("kPi")){
		data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_Pi"));
		eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_Pi"));
		mass= 0.13957;//(GeV - pions)
// 		ccorrdata =new TCanvas("Pions ccorrdata","Pions - corrected data",0,0,600,800);
		printf("\n\n**************************************\n");
		printf("\tRunning for pions!\n");
	}
	else if(fPartType.Contains("kK")){
		data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_K"));
		eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_K"));
		mass= 0.49368;//(GeV - kaons)
// 		ccorrdata =new TCanvas("Kaons ccorrdata","Kaons - corrected data",0,0,600,800);
		printf("\n\n**************************************\n");
 		printf("\tRunning for kaons!\n");
	}
	else if(fPartType.Contains("kProton")){
		data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_P"));
		eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_P"));
		mass = 0.93827;//(GeV - protons)
// 		ccorrdata =new TCanvas("Protons ccorrdata","Protons - corrected data",0,0,600,800);
		printf("\n\n**************************************\n");
 		printf("\tRunning for protons!\n");
	}
	else printf("Unsupported particle type!\n");

    if(!data){
		printf("Unable to get CFContainer! \n");
		return;
	}
	
	if(!eff){
		printf("No Eff Grid found! \n");
		return;
    }
	
// 	TCanvas *ccorrdata = new TCanvas();
// 	ccorrdata->Divide(1,2);
// 	ccorrdata->cd(1);
// 	ccorrdata->cd(1)->SetLogy();

 
    AliCFDataGrid *corrdata = new AliCFDataGrid("corrdata","corrected data",*data, stepRec);

//correct selection step "reconstructed"
//    corrdata->SetMeasured(stepRec); //set data to be corrected
    corrdata->ApplyEffCorrection(*eff);//apply the correction for efficiency

//     TH1D *hPtMC = data->ShowProjection(0,0); //MC distribution ShowProjection(ivar, istep)
//     hPtMC->SetMarkerStyle(20);
//     hPtMC->SetMarkerColor(kGreen-3);
//     hPtMC->GetXaxis()->SetTitle("p_{T}(GeV/c)");
//     hPtMC->GetYaxis()->SetTitle("#frac{dN}{p_{T}dp_{T}}");
//     hPtMC->Draw("p e1");
// 
//     TH1D *hPtESDI = corrdata->GetData()->Project(0); //uncorrected ESD  Project(ivar)
//     hPtESDI->SetMarkerStyle(25);
//     hPtESDI->SetMarkerColor(kBlue);
//     hPtESDI->GetXaxis()->SetTitle("Pt");
//     hPtESDI->Draw("p e1 same");

    TH1D *hPtESDCorr = (TH1D*)corrdata->Project(0); //corrected data
    hPtESDCorr->SetMarkerStyle(26);
    hPtESDCorr->SetMarkerColor(kRed);        //ESD corrected 
	hPtESDCorr->SetName("hPtESDCorr");
// 	hPtESDCorr->Draw("p e1 same");

// 	TLegend *leg2 = new TLegend(0.40,0.65,0.87,0.8);
//     leg2->AddEntry(hPtMC," generated","p");
//     leg2->AddEntry(hPtESDI," reconstructed (not corrected)","p");    
//     leg2->AddEntry(hPtESDCorr," reconstructed & corrected","p");    
// 	leg2->Draw();


// 	ccorrdata->cd(2);
// 	ccorrdata->cd(2)->SetLogy();
// 	
// 	TH1D *efficiency = eff->Project(0); //the efficiency vs pt (ipt = 0)
// 	efficiency->GetXaxis()->SetTitle("p_{T}(GeV/c)");
// 	efficiency->Draw();


	TH1D *hNoEvt = dynamic_cast<TH1D*>(fInputList->FindObject("TaskCentral_NoEvt"));
	if(!hNoEvt){
		printf("Unable to get the number of events! \n");
		return;
	}

	
    Int_t noEvt = (Int_t)(hNoEvt->GetEntries());

    printf("\n** No of processed events = %i **\n",noEvt);

    Double_t scale = 1.0/noEvt;

    TH1D *hPtESDCorrNorm = (TH1D*)hPtESDCorr->Clone("hPtESDCorrNorm");
    hPtESDCorrNorm->Scale(scale);

	fResultsList->SetName(fPartType);
	fResultsList->Add(hPtESDCorr);
	fResultsList->Add(hPtESDCorrNorm);
	
}


//*********************************************************************************
//_________________________________________________________________________________
// fit Tsallis


Double_t Integ1(const Double_t *x){
//the function under the integral - TBW

    const Double_t rMax = 11.0;//(fm)

    Double_t betaR = betaS*(x[0]/rMax);


    Double_t rho = TMath::ATanH(betaR);

    
    Double_t temp = 1.0+((q-1.0)*(mt*TMath::CosH(x[2])*TMath::CosH(rho) - pt*TMath::SinH(rho)*TMath::Cos(x[1]))/T);

    Double_t power = -1.0/(q-1.0);

    Double_t temp1 = pow(temp, power);

    Double_t f = TMath::CosH(x[1])*x[0]*temp1;

    return f;
}


// TF1 requires the function to have the ( )( Double_t *, Double_t *) signature 
Double_t Tsallis(Double_t* const x, Double_t* const par){ 
//computes the triple integral in the TBW formula

    pt = x[0]; 

    T = par[0];
    betaS = par[1];
    q = par[2];

    mt = sqrt(mass*mass+ pt*pt);
	

   ROOT::Math::WrappedMultiFunction<> f1(Integ1,3);

    ROOT::Math::IntegratorMultiDim ig(f1, ROOT::Math::IntegrationMultiDim::kPLAIN);


    Double_t xmin[3] = {0.0, -TMath::Pi(), -0.5}; //rmin, phi_min, ymin
    Double_t xmax[3] = {11.0, TMath::Pi(), 0.5}; //rmax, phi_max, ymax

    Double_t rez = mt*ig.Integral(xmin,xmax);

    return rez;
}


void AliAnalysisCentralExtrapolate::TsallisFit(){
// fits and extrpolates the pt spectrum using TBW model

	printf("---------------------------------------------\n");
	printf("*****   Tsallis Blast Wave Fit    *****\n");

	printf("AliAnalysisCentralExtrapolate::Tsallis mass = %g\n", mass);


    TStopwatch timer; 

    TH1::SetDefaultSumw2();
    TH1::AddDirectory(0);

    timer.Start(); 

	TH1D *ptSpectra  = dynamic_cast<TH1D*>(fResultsList->FindObject("hPtESDCorrNorm"));
	if(!ptSpectra){
		printf("TsallisFit: Can't get the normalized spectrum\n");
		return;	
	}

	if(!ptSpectra->GetEntries()){
		printf("TsallisFit: The fit data is empty!\n");
		return;
	}

    TVirtualFitter::SetDefaultFitter("Minuit2");//options Minuit, Minuit2, Fumili, Fumili2

    TF1 *ptFit = new TF1("ptFit",Tsallis,0.2,4.0,3);

    gStyle->SetOptFit(1112);

//     ptFit->SetParName(0,"T");
//     ptFit->SetParName(1,"betaS");
//     ptFit->SetParName(2,"q");

    ptFit->SetParameters(0.1,0.5,1.05); //GeV
    ptFit->SetParLimits(0,0.0,0.5);//GeV
    ptFit->SetParLimits(1,0.0,1.0);
    ptFit->SetParLimits(2,1.00001,1.5);

	ptSpectra->Fit("ptFit","Q R B ME N");

    timer.Stop(); 
    timer.Print();

//  ----------- print the fit result ----------
    Double_t chi2 = ptFit->GetChisquare();
    Double_t ndf = ptFit->GetNDF();
    Double_t p1 = ptFit->GetParameter(0);
    Double_t param1Err = ptFit->GetParError(0);
    Double_t p2 = ptFit->GetParameter(1);
    Double_t param2Err = ptFit->GetParError(1);
    Double_t p3 = ptFit->GetParameter(2);
    Double_t param3Err = ptFit->GetParError(2);

    printf("chi2/ndf = %f/%f\n", chi2,ndf);
    printf("p1 = %f\tparam1Err = %f\n", p1, param1Err);
    printf("p2 = %f\tparam2Err = %f\n", p2, param2Err);
    printf("p3 = %f\tparam3Err = %f\n", p3, param3Err);


// create a new, "extended" histogram
    TH1D *hPtExtTsallis = new TH1D("PtExtTsallis","Pt Corr Norm Ext",25,0.0,5.0);
    
    Double_t bin, binerr, test;
    
    for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){
	
		bin = ptSpectra->GetBinContent(i);
		binerr = ptSpectra->GetBinError(i);
		test = ptSpectra->GetBinLowEdge(i);

		hPtExtTsallis->SetBinContent(i, bin);
		hPtExtTsallis->SetBinError(i, binerr);
    }
    
    Double_t eval;
    eval = ptFit->Eval(0.1);
    printf("extrapolated pt value = %g \n", eval);
    printf("****************************************************\n\n");
	
    hPtExtTsallis->SetBinContent(1, eval);

    hPtExtTsallis->SetMarkerStyle(24);
    hPtExtTsallis->SetMarkerColor(kRed);        //ESD corrected 
    hPtExtTsallis->GetXaxis()->SetTitle("Pt");


	fResultsList->Add(hPtExtTsallis);

}

//*********************************************************************************
//_________________________________________________________________________________
// fit Boltzmann

Double_t func(Double_t r){ 
//the function under the integral - BGBW

    const Double_t rMax = 11.0;//(fm)

    Double_t betaR = betaS*(r/rMax);


    Double_t rho = TMath::ATanH(betaR);


	mt = sqrt(mass*mass + pt*pt);// !!! par[0]= pt !!!!!

    Double_t argI0 = (pt*TMath::SinH(rho))/T; //T = par[1]
	if(argI0>700.0) return 0.0; // !! floating point exception protection

    Double_t argK1 = (mt*TMath::CosH(rho))/T;

	
    Double_t i0 = TMath::BesselI0(argI0);


    Double_t k1 = TMath::BesselK1(argK1);


    Double_t f = r*mt*i0*k1;
	
    return f;
}


Double_t Boltzmann(Double_t* const x, Double_t* const par){
//computes the integral in the BGBW formula

	pt = x[0];
    T = par[0];
	
    betaS = par[1];


    ROOT::Math::WrappedFunction<> f1(func);

    ROOT::Math::Integrator ig(f1, ROOT::Math::IntegrationOneDim::kGAUSS,1.E-12,1.E-12,10000);

    Double_t rez = ig.Integral(0.0,11.0);

    return rez;
}


void AliAnalysisCentralExtrapolate::BoltzmannFit(){
//fits and extrapoates the pt spectrum using the BGBW model

	printf("---------------------------------------------------\n");
	printf("*** Boltzmann-Gibbs Blast Wave Fit  ***\n");

	printf("AliAnalysisCentralExtrapolate::Boltzmann mass = %g\n", mass);


    TStopwatch timer; 

    TH1::SetDefaultSumw2();
    TH1::AddDirectory(0);

    timer.Start(); 


	TH1D *ptSpectra  = dynamic_cast<TH1D*>(fResultsList->FindObject("hPtESDCorrNorm"));
	if(!ptSpectra){
		printf("BoltzmannFit: Can't get the normalized spectrum\n");
		return;	
	}

	printf("pt spectra get entries: %f\n",ptSpectra->GetEntries());

	if(!ptSpectra->GetEntries()){
		printf("BoltzmannFit: The fit data is empty!\n");
		return;
	}

    gStyle->SetOptStat("neRM");

    TVirtualFitter::SetDefaultFitter("Minuit2");

    TF1 *ptFit = new TF1("ptFit",Boltzmann,0.2,4.0,2); 

    gStyle->SetOptFit(1112);

//     ptFit->SetParName(0,"T");
//     ptFit->SetParName(1,"betaS");

    ptFit->SetParameters(0.1,0.5); //GeV
    ptFit->SetParLimits(0,0.001,0.5);//GeV
    ptFit->SetParLimits(1,0.0,1.0);

	ptSpectra->Fit("ptFit","Q R B ME I N");

    timer.Stop();
    timer.Print();

//  ----------- print the fit results ----------
    Double_t chi2 = ptFit->GetChisquare();
    Double_t ndf = ptFit->GetNDF();
    Double_t p1 = ptFit->GetParameter(0);
    Double_t param1Err = ptFit->GetParError(0);
    Double_t p2 = ptFit->GetParameter(1);
    Double_t param2Err = ptFit->GetParError(1);

    printf("chi2/ndf = %f/%f = %f\n", chi2,ndf,chi2/ndf);
    printf("p1 = %f (GeV)\tparam1Err = %f\n", p1, param1Err);
    printf("p2 = %f\tparam2Err = %f\n", p2, param2Err);


    TH1D *hPtExtBoltzmann = new TH1D("PtExtBoltzmann","Pt Corr Norm Ext",25,0.0,5.0);

    Double_t bin, binerr, test;

    for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){
	
		bin = ptSpectra->GetBinContent(i);
		binerr = ptSpectra->GetBinError(i);
		test = ptSpectra->GetBinLowEdge(i);
	
		hPtExtBoltzmann->SetBinContent(i, bin);
		hPtExtBoltzmann->SetBinError(i, binerr);
    }

    Double_t eval;
    eval = ptFit->Eval(0.1);
    printf("extrapolated pt value = %g \n", eval);
    printf("********************************************\n\n");

    hPtExtBoltzmann->SetBinContent(1, eval);

    hPtExtBoltzmann->SetMarkerStyle(24);
    hPtExtBoltzmann->SetMarkerColor(kRed);        //ESD corrected 
    hPtExtBoltzmann->GetXaxis()->SetTitle("Pt");


	fResultsList->Add(hPtExtBoltzmann);

}