[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliAnalysisCentralExtrapolate.cxx

/*************************************************************************
* 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("AliAnalysisCentralExtrapolate::AliAnalysisCentralExtrapolate(const char *name)\n");
	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);

//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 = 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);


	Double_t 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);

}
Commit	Line	Data
bde49c8a	1	/*************************************************************************
	2	* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
410ff8cc	16	// -------------------------------------------------------
	17	// pt spectra extrapolation in the 0. - 0.2 region using
	18	// Boltzmann-Gibbs Blast Wave model or Tsallis Blast
	19	// Wave model for azimuthal isotropic expansion in
	20	// highly central collisions analysis
	21	// author: Cristian Andrei
	22	// acristian@niham.nipne.ro
	23	// ----------------------------------------------------------
	24
bde49c8a	25
	26	#include "TH1D.h"
	27	#include "TFile.h"
	28	#include "TList.h"
410ff8cc	29	#include "TCanvas.h"
410ff8cc	30	#include "TLegend.h"
bde49c8a	31	// #include "TString.h"
	32	// #include "TPaveStats.h"
	33	#include "TMath.h"
	34	#include "TStopwatch.h"
	35	#include "TStyle.h"
	36	#include "TF1.h"
	37
	38	#include "TVirtualFitter.h"
	39	#include "Math/WrappedFunction.h"
	40	#include "Math/Integrator.h"
	41	#include "Math/IntegratorMultiDim.h"
	42
	43	#include "AliCFContainer.h"
	44	#include "AliCFDataGrid.h"
	45	#include "AliCFEffGrid.h"
	46
	47	#include "AliAnalysisCentralExtrapolate.h"
	48
	49
	50	ClassImp(AliAnalysisCentralExtrapolate)
	51
	52	//________________________________________________________________________
	53	AliAnalysisCentralExtrapolate::AliAnalysisCentralExtrapolate(const char *name)
	54	: TObject()
	55	,fPartType()
410ff8cc	56	,fInputList(0)
bde49c8a	57	,fResultsList(0x0)
	58
	59	{
	60	// Constructor
	61
	62	// printf("AliAnalysisCentralExtrapolate::AliAnalysisCentralExtrapolate(const char *name)\n");
410ff8cc	63	fInputList = new TList();
bde49c8a	64	fResultsList = new TList();
	65
	66	}
	67
	68	AliAnalysisCentralExtrapolate::~AliAnalysisCentralExtrapolate() {
	69
	70	if(fInputList) delete fInputList;
	71	if(fResultsList) delete fResultsList;
	72
	73	}
	74
	75	//************************************************************************************
	76	//____________________________________________________________________________________
	77	// using global variables to avoid the limitations of ROOT::Math::WrappedMultiFunction<>
	78	static Double_t mass = 0.;
	79	static Double_t pt = 0.;
	80	static Double_t T = 0.;
	81	static Double_t betaS = 0.;
	82	static Double_t q = 0.;
	83	static Double_t mt = 0.;
	84
	85	void AliAnalysisCentralExtrapolate::ApplyEff(){
	86	// applies the efficiency map to the pt spectra
	87
	88	TH1::SetDefaultSumw2();
	89	// Int_t stepGen = 0;
	90	Int_t stepRec = 1;
	91
	92	AliCFContainer *data = 0;
410ff8cc	93
410ff8cc	94
bde49c8a	95	TFile *file1 = new TFile("$ALICE_ROOT/PWG2/data/AliAnalysisCentralEfficiency.root", "read");
	96
	97	AliCFEffGrid *eff = 0;
	98
410ff8cc	99	if(fPartType.Contains("kPi")){
410ff8cc	100	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_Pi"));
bde49c8a	101	eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_Pi"));
	102	mass= 0.13957;//(GeV - pions)
	103	// ccorrdata =new TCanvas("Pions ccorrdata","Pions - corrected data",0,0,600,800);
	104	printf("\n\n**************************************\n");
	105	printf("\tRunning for pions!\n");
	106	}
	107	else if(fPartType.Contains("kK")){
410ff8cc	108	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_K"));
bde49c8a	109	eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_K"));
	110	mass= 0.49368;//(GeV - kaons)
	111	// ccorrdata =new TCanvas("Kaons ccorrdata","Kaons - corrected data",0,0,600,800);
	112	printf("\n\n**************************************\n");
	113	printf("\tRunning for kaons!\n");
	114	}
	115	else if(fPartType.Contains("kProton")){
410ff8cc	116	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_P"));
bde49c8a	117	eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_P"));
	118	mass = 0.93827;//(GeV - protons)
	119	// ccorrdata =new TCanvas("Protons ccorrdata","Protons - corrected data",0,0,600,800);
	120	printf("\n\n**************************************\n");
	121	printf("\tRunning for protons!\n");
	122	}
	123	else printf("Unsupported particle type!\n");
	124
	125	if(!data){
	126	printf("Unable to get CFContainer! \n");
8a993eac	127	return;
bde49c8a	128	}
	129
	130	if(!eff){
	131	printf("No Eff Grid found! \n");
8a993eac	132	return;
bde49c8a	133	}
410ff8cc	134
410ff8cc	135	// TCanvas *ccorrdata = new TCanvas();
bde49c8a	136	// ccorrdata->Divide(1,2);
bde49c8a	137	// ccorrdata->cd(1);
410ff8cc	138	// ccorrdata->cd(1)->SetLogy();
bde49c8a	139
	140
	141	AliCFDataGrid corrdata = new AliCFDataGrid("corrdata","corrected data",data);
	142
	143	//correct selection step "reconstructed"
	144	corrdata->SetMeasured(stepRec); //set data to be corrected
	145	corrdata->ApplyEffCorrection(*eff);//apply the correction for efficiency
	146
410ff8cc	147	// TH1D *hPtMC = data->ShowProjection(0,0); //MC distribution ShowProjection(ivar, istep)
	148	// hPtMC->SetMarkerStyle(20);
	149	// hPtMC->SetMarkerColor(kGreen-3);
	150	// hPtMC->GetXaxis()->SetTitle("p_{T}(GeV/c)");
	151	// hPtMC->GetYaxis()->SetTitle("#frac{dN}{p_{T}dp_{T}}");
bde49c8a	152	// hPtMC->Draw("p e1");
410ff8cc	153	//
	154	// TH1D *hPtESDI = corrdata->GetData()->Project(0); //uncorrected ESD Project(ivar)
	155	// hPtESDI->SetMarkerStyle(25);
	156	// hPtESDI->SetMarkerColor(kBlue);
	157	// hPtESDI->GetXaxis()->SetTitle("Pt");
bde49c8a	158	// hPtESDI->Draw("p e1 same");
	159
	160	TH1D *hPtESDCorr = corrdata->Project(0); //corrected data
	161	hPtESDCorr->SetMarkerStyle(26);
	162	hPtESDCorr->SetMarkerColor(kRed); //ESD corrected
	163	hPtESDCorr->SetName("hPtESDCorr");
	164	// hPtESDCorr->Draw("p e1 same");
	165
	166	// TLegend *leg2 = new TLegend(0.40,0.65,0.87,0.8);
	167	// leg2->AddEntry(hPtMC," generated","p");
	168	// leg2->AddEntry(hPtESDI," reconstructed (not corrected)","p");
	169	// leg2->AddEntry(hPtESDCorr," reconstructed & corrected","p");
	170	// leg2->Draw();
	171
	172
	173	// ccorrdata->cd(2);
	174	// ccorrdata->cd(2)->SetLogy();
410ff8cc	175	//
	176	// TH1D *efficiency = eff->Project(0); //the efficiency vs pt (ipt = 0)
	177	// efficiency->GetXaxis()->SetTitle("p_{T}(GeV/c)");
bde49c8a	178	// efficiency->Draw();
	179
	180
0fb1f0cf	181	TH1D hNoEvt = dynamic_cast<TH1D>(fInputList->FindObject("TaskCentral_NoEvt"));
bde49c8a	182	if(!hNoEvt){
bde49c8a	183	printf("Unable to get the number of events! \n");
8a993eac	184	return;
bde49c8a	185	}
bde49c8a	186
410ff8cc	187
bde49c8a	188	Int_t noEvt = (Int_t)(hNoEvt->GetEntries());
	189
	190	printf("\n No of processed events = %i \n",noEvt);
8a993eac	191
bde49c8a	192	Double_t scale = 1.0/noEvt;
8a993eac	193
bde49c8a	194	TH1D hPtESDCorrNorm = (TH1D)hPtESDCorr->Clone("hPtESDCorrNorm");
	195	hPtESDCorrNorm->Scale(scale);
	196
	197	fResultsList->SetName(fPartType);
	198	fResultsList->Add(hPtESDCorr);
	199	fResultsList->Add(hPtESDCorrNorm);
	200
	201	}
	202
	203
	204	//*********************************************************************************
	205	//_________________________________________________________________________________
	206	// fit Tsallis
	207
	208
	209	Double_t Integ1(const Double_t *x){
	210	//the function under the integral - TBW
	211
	212	const Double_t rMax = 11.0;//(fm)
	213
	214	Double_t betaR = betaS*(x[0]/rMax);
	215
	216
	217	Double_t rho = TMath::ATanH(betaR);
	218
	219
	220	Double_t temp = 1.0+((q-1.0)(mtTMath::CosH(x[2])TMath::CosH(rho) - ptTMath::SinH(rho)*TMath::Cos(x[1]))/T);
	221
	222	Double_t power = -1.0/(q-1.0);
	223
	224	Double_t temp1 = pow(temp, power);
	225
	226	Double_t f = TMath::CosH(x[1])x[0]temp1;
	227
	228	return f;
	229	}
	230
	231
	232	// TF1 requires the function to have the ( )( Double_t , Double_t ) signature
	233	Double_t Tsallis(Double_t* const x, Double_t* const par){
	234	//computes the triple integral in the TBW formula
	235
	236	pt = x[0];
	237
	238	T = par[0];
	239	betaS = par[1];
	240	q = par[2];
	241
	242	mt = sqrt(massmass+ ptpt);
	243
	244
	245	ROOT::Math::WrappedMultiFunction<> f1(Integ1,3);
	246
	247	ROOT::Math::IntegratorMultiDim ig(f1, ROOT::Math::IntegrationMultiDim::kPLAIN);
	248
	249
	250	Double_t xmin[3] = {0.0, -TMath::Pi(), -0.5}; //rmin, phi_min, ymin
	251	Double_t xmax[3] = {11.0, TMath::Pi(), 0.5}; //rmax, phi_max, ymax
	252
	253	Double_t rez = mt*ig.Integral(xmin,xmax);
	254
	255	return rez;
	256	}
	257
258
259	void AliAnalysisCentralExtrapolate::TsallisFit(){
260	// fits and extrpolates the pt spectrum using TBW model
261
262	printf("---------------------------------------------\n");
263	printf("*** Tsallis Blast Wave Fit ***\n");
264
265	printf("AliAnalysisCentralExtrapolate::Tsallis mass = %g\n", mass);
266
267
268	TStopwatch timer;
269
270	TH1::SetDefaultSumw2();
271	TH1::AddDirectory(0);
272
273	timer.Start();
274
275	TH1D ptSpectra = dynamic_cast<TH1D>(fResultsList->FindObject("hPtESDCorrNorm"));
8a993eac	276	if(!ptSpectra){
410ff8cc	277	printf("TsallisFit: Can't get the normalized spectrum\n");
8a993eac	278	return;
	279	}
	280
	281	if(!ptSpectra->GetEntries()){
410ff8cc	282	printf("TsallisFit: The fit data is empty!\n");
8a993eac	283	return;
8a993eac	284	}
bde49c8a	285
	286	TVirtualFitter::SetDefaultFitter("Minuit2");//options Minuit, Minuit2, Fumili, Fumili2
	287
	288	TF1 *ptFit = new TF1("ptFit",Tsallis,0.2,4.0,3);
	289
	290	gStyle->SetOptFit(1112);
	291
410ff8cc	292	// ptFit->SetParName(0,"T");
	293	// ptFit->SetParName(1,"betaS");
	294	// ptFit->SetParName(2,"q");
bde49c8a	295
	296	ptFit->SetParameters(0.1,0.5,1.05); //GeV
	297	ptFit->SetParLimits(0,0.0,0.5);//GeV
	298	ptFit->SetParLimits(1,0.0,1.0);
	299	ptFit->SetParLimits(2,1.00001,1.5);
	300
	301	ptSpectra->Fit("ptFit","Q R B ME N");
	302
	303	timer.Stop();
	304	timer.Print();
	305
	306	// ----------- print the fit result ----------
	307	Double_t chi2 = ptFit->GetChisquare();
	308	Double_t ndf = ptFit->GetNDF();
	309	Double_t p1 = ptFit->GetParameter(0);
	310	Double_t param1Err = ptFit->GetParError(0);
	311	Double_t p2 = ptFit->GetParameter(1);
	312	Double_t param2Err = ptFit->GetParError(1);
	313	Double_t p3 = ptFit->GetParameter(2);
	314	Double_t param3Err = ptFit->GetParError(2);
	315
	316	printf("chi2/ndf = %f/%f\n", chi2,ndf);
	317	printf("p1 = %f\tparam1Err = %f\n", p1, param1Err);
	318	printf("p2 = %f\tparam2Err = %f\n", p2, param2Err);
	319	printf("p3 = %f\tparam3Err = %f\n", p3, param3Err);
	320
	321
	322	// create a new, "extended" histogram
0fb1f0cf	323	TH1D *hPtExtTsallis = new TH1D("PtExtTsallis","Pt Corr Norm Ext",25,0.0,5.0);
bde49c8a	324
	325	Double_t bin, binerr, test;
	326
	327	for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){
	328
	329	bin = ptSpectra->GetBinContent(i);
	330	binerr = ptSpectra->GetBinError(i);
	331	test = ptSpectra->GetBinLowEdge(i);
	332
	333	hPtExtTsallis->SetBinContent(i, bin);
	334	hPtExtTsallis->SetBinError(i, binerr);
	335	}
	336
	337	Double_t eval;
	338	eval = ptFit->Eval(0.1);
	339	printf("extrapolated pt value = %g \n", eval);
	340	printf("****************************************************\n\n");
	341
	342	hPtExtTsallis->SetBinContent(1, eval);
	343
	344	hPtExtTsallis->SetMarkerStyle(24);
	345	hPtExtTsallis->SetMarkerColor(kRed); //ESD corrected
	346	hPtExtTsallis->GetXaxis()->SetTitle("Pt");
	347
	348
	349	fResultsList->Add(hPtExtTsallis);
	350
	351	}
	352
	353	//*********************************************************************************
	354	//_________________________________________________________________________________
	355	// fit Boltzmann
	356
	357	Double_t func(Double_t r){
	358	//the function under the integral - BGBW
	359
	360	const Double_t rMax = 11.0;//(fm)
	361
	362	Double_t betaR = betaS*(r/rMax);
	363
410ff8cc	364
bde49c8a	365	Double_t rho = TMath::ATanH(betaR);
bde49c8a	366
410ff8cc	367
bde49c8a	368	Double_t mt = sqrt(massmass + ptpt);// !!! par[0]= pt !!!!!
410ff8cc	369
bde49c8a	370	Double_t argI0 = (pt*TMath::SinH(rho))/T; //T = par[1]
410ff8cc	371	if(argI0>700.0) return 0.0; // !! floating point exception protection
bde49c8a	372
	373	Double_t argK1 = (mt*TMath::CosH(rho))/T;
	374
410ff8cc	375
bde49c8a	376	Double_t i0 = TMath::BesselI0(argI0);
bde49c8a	377
410ff8cc	378
bde49c8a	379	Double_t k1 = TMath::BesselK1(argK1);
bde49c8a	380
bde49c8a	381
410ff8cc	382	Double_t f = rmti0*k1;
410ff8cc	383
bde49c8a	384	return f;
	385	}
	386
	387
	388	Double_t Boltzmann(Double_t* const x, Double_t* const par){
	389	//computes the integral in the BGBW formula
	390
	391	pt = x[0];
	392	T = par[0];
	393
	394	betaS = par[1];
	395
	396
	397	ROOT::Math::WrappedFunction<> f1(func);
	398
	399	ROOT::Math::Integrator ig(f1, ROOT::Math::IntegrationOneDim::kGAUSS,1.E-12,1.E-12,10000);
	400
	401	Double_t rez = ig.Integral(0.0,11.0);
	402
	403	return rez;
	404	}
	405
	406
	407	void AliAnalysisCentralExtrapolate::BoltzmannFit(){
	408	//fits and extrapoates the pt spectrum using the BGBW model
	409
	410	printf("---------------------------------------------------\n");
	411	printf("* Boltzmann-Gibbs Blast Wave Fit *\n");
	412
	413	printf("AliAnalysisCentralExtrapolate::Boltzmann mass = %g\n", mass);
	414
	415
	416	TStopwatch timer;
	417
	418	TH1::SetDefaultSumw2();
	419	TH1::AddDirectory(0);
8a993eac	420
bde49c8a	421	timer.Start();
	422
	423
	424	TH1D ptSpectra = dynamic_cast<TH1D>(fResultsList->FindObject("hPtESDCorrNorm"));
8a993eac	425	if(!ptSpectra){
	426	printf("BoltzmannFit: Can't get the normalized spectrum\n");
	427	return;
	428	}
	429
410ff8cc	430	printf("pt spectra get entries: %f\n",ptSpectra->GetEntries());
410ff8cc	431
8a993eac	432	if(!ptSpectra->GetEntries()){
	433	printf("BoltzmannFit: The fit data is empty!\n");
	434	return;
	435	}
bde49c8a	436
	437	gStyle->SetOptStat("neRM");
	438
	439	TVirtualFitter::SetDefaultFitter("Minuit2");
	440
	441	TF1 *ptFit = new TF1("ptFit",Boltzmann,0.2,4.0,2);
	442
	443	gStyle->SetOptFit(1112);
	444
410ff8cc	445	// ptFit->SetParName(0,"T");
410ff8cc	446	// ptFit->SetParName(1,"betaS");
bde49c8a	447
bde49c8a	448	ptFit->SetParameters(0.1,0.5); //GeV
410ff8cc	449	ptFit->SetParLimits(0,0.001,0.5);//GeV
bde49c8a	450	ptFit->SetParLimits(1,0.0,1.0);
410ff8cc	451
bde49c8a	452	ptSpectra->Fit("ptFit","Q R B ME I N");
bde49c8a	453
410ff8cc	454	timer.Stop();
bde49c8a	455	timer.Print();
	456
	457	// ----------- print the fit results ----------
	458	Double_t chi2 = ptFit->GetChisquare();
	459	Double_t ndf = ptFit->GetNDF();
	460	Double_t p1 = ptFit->GetParameter(0);
	461	Double_t param1Err = ptFit->GetParError(0);
	462	Double_t p2 = ptFit->GetParameter(1);
	463	Double_t param2Err = ptFit->GetParError(1);
	464
	465	printf("chi2/ndf = %f/%f = %f\n", chi2,ndf,chi2/ndf);
	466	printf("p1 = %f (GeV)\tparam1Err = %f\n", p1, param1Err);
	467	printf("p2 = %f\tparam2Err = %f\n", p2, param2Err);
	468
	469
0fb1f0cf	470	TH1D *hPtExtBoltzmann = new TH1D("PtExtBoltzmann","Pt Corr Norm Ext",25,0.0,5.0);
bde49c8a	471
	472	Double_t bin, binerr, test;
	473
	474	for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){
	475
	476	bin = ptSpectra->GetBinContent(i);
	477	binerr = ptSpectra->GetBinError(i);
	478	test = ptSpectra->GetBinLowEdge(i);
	479
	480	hPtExtBoltzmann->SetBinContent(i, bin);
	481	hPtExtBoltzmann->SetBinError(i, binerr);
	482	}
	483
	484	Double_t eval;
	485	eval = ptFit->Eval(0.1);
	486	printf("extrapolated pt value = %g \n", eval);
	487	printf("********************************************\n\n");
	488
	489	hPtExtBoltzmann->SetBinContent(1, eval);
	490
	491	hPtExtBoltzmann->SetMarkerStyle(24);
	492	hPtExtBoltzmann->SetMarkerColor(kRed); //ESD corrected
	493	hPtExtBoltzmann->GetXaxis()->SetTitle("Pt");
	494
	495
	496	fResultsList->Add(hPtExtBoltzmann);
	497
	498	}
	499
	500