[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(0x0)
  ,fResultsList(0x0)

{
// Constructor

//  printf("AliAnalysisCentralExtrapolate::AliAnalysisCentralExtrapolate(const char *name)\n");

	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;
    
//     TCanvas *ccorrdata;
    	
	if(fPartType.Contains("kPi")){
		data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("AliAnalysisCentral_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("AliAnalysisCentral_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("AliAnalysisCentral_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");
		exit(1);
	}
	
	if(!eff){
		printf("No Eff Grid found! \n");
		exit(1);
    }
 
// 	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();


//scalez hist corectata la nr de ev

	TH1F *hNoEvt = dynamic_cast<TH1F*>(fInputList->FindObject("AliAnalysisCentral_NoEvt"));
	if(!hNoEvt){
		printf("Unable to get the number of events! \n");
		exit(1);
	}

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

    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
    TH1F *hPtExtTsallis = new TH1F("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]

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

    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.0,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);


    TH1F *hPtExtBoltzmann = new TH1F("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
	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
	25
	26	#include "TH1D.h"
	27	#include "TFile.h"
	28	#include "TList.h"
	29	// #include "TCanvas.h"
	30	// #include "TLegend.h"
	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()
	56	,fInputList(0x0)
	57	,fResultsList(0x0)
	58
	59	{
	60	// Constructor
	61
	62	// printf("AliAnalysisCentralExtrapolate::AliAnalysisCentralExtrapolate(const char *name)\n");
	63
	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;
93
94	TFile *file1 = new TFile("$ALICE_ROOT/PWG2/data/AliAnalysisCentralEfficiency.root", "read");
95
96	AliCFEffGrid *eff = 0;
97
98	// TCanvas *ccorrdata;
99
100	if(fPartType.Contains("kPi")){
101	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("AliAnalysisCentral_CFCont_Pi"));
102	eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_Pi"));
103	mass= 0.13957;//(GeV - pions)
104	// ccorrdata =new TCanvas("Pions ccorrdata","Pions - corrected data",0,0,600,800);
105	printf("\n\n**************************************\n");
106	printf("\tRunning for pions!\n");
107	}
108	else if(fPartType.Contains("kK")){
109	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("AliAnalysisCentral_CFCont_K"));
110	eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_K"));
111	mass= 0.49368;//(GeV - kaons)
112	// ccorrdata =new TCanvas("Kaons ccorrdata","Kaons - corrected data",0,0,600,800);
113	printf("\n\n**************************************\n");
114	printf("\tRunning for kaons!\n");
115	}
116	else if(fPartType.Contains("kProton")){
117	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("AliAnalysisCentral_CFCont_P"));
118	eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_P"));
119	mass = 0.93827;//(GeV - protons)
120	// ccorrdata =new TCanvas("Protons ccorrdata","Protons - corrected data",0,0,600,800);
121	printf("\n\n**************************************\n");
122	printf("\tRunning for protons!\n");
123	}
124	else printf("Unsupported particle type!\n");
125
126	if(!data){
127	printf("Unable to get CFContainer! \n");
128	exit(1);
129	}
130
131	if(!eff){
132	printf("No Eff Grid found! \n");
133	exit(1);
134	}
135
136	// ccorrdata->Divide(1,2);
137	// ccorrdata->cd(1);
138	// ccorrdata->cd(1)->SetLogy();
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
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}}");
152	// hPtMC->Draw("p e1");
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");
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();
175
176	TH1D *efficiency = eff->Project(0); //the efficiency vs pt (ipt = 0)
177	efficiency->GetXaxis()->SetTitle("p_{T}(GeV/c)");
178	// efficiency->Draw();
179
180
181	//scalez hist corectata la nr de ev
182
183	TH1F hNoEvt = dynamic_cast<TH1F>(fInputList->FindObject("AliAnalysisCentral_NoEvt"));
184	if(!hNoEvt){
185	printf("Unable to get the number of events! \n");
186	exit(1);
187	}
188
189	Int_t noEvt = (Int_t)(hNoEvt->GetEntries());
190
191	printf("\n No of processed events = %i \n",noEvt);
192
193	Double_t scale = 1.0/noEvt;
194
195	TH1D hPtESDCorrNorm = (TH1D)hPtESDCorr->Clone("hPtESDCorrNorm");
196	hPtESDCorrNorm->Scale(scale);
197
198	fResultsList->SetName(fPartType);
199	fResultsList->Add(hPtESDCorr);
200	fResultsList->Add(hPtESDCorrNorm);
201
202	}
203
204
205	//*********************************************************************************
206	//_________________________________________________________________________________
207	// fit Tsallis
208
209
210	Double_t Integ1(const Double_t *x){
211	//the function under the integral - TBW
212
213	const Double_t rMax = 11.0;//(fm)
214
215	Double_t betaR = betaS*(x[0]/rMax);
216
217
218	Double_t rho = TMath::ATanH(betaR);
219
220
221	Double_t temp = 1.0+((q-1.0)(mtTMath::CosH(x[2])TMath::CosH(rho) - ptTMath::SinH(rho)*TMath::Cos(x[1]))/T);
222
223	Double_t power = -1.0/(q-1.0);
224
225	Double_t temp1 = pow(temp, power);
226
227	Double_t f = TMath::CosH(x[1])x[0]temp1;
228
229	return f;
230	}
231
232
233	// TF1 requires the function to have the ( )( Double_t , Double_t ) signature
234	Double_t Tsallis(Double_t* const x, Double_t* const par){
235	//computes the triple integral in the TBW formula
236
237	pt = x[0];
238
239	T = par[0];
240	betaS = par[1];
241	q = par[2];
242
243	mt = sqrt(massmass+ ptpt);
244
245
246	ROOT::Math::WrappedMultiFunction<> f1(Integ1,3);
247
248	ROOT::Math::IntegratorMultiDim ig(f1, ROOT::Math::IntegrationMultiDim::kPLAIN);
249
250
251	Double_t xmin[3] = {0.0, -TMath::Pi(), -0.5}; //rmin, phi_min, ymin
252	Double_t xmax[3] = {11.0, TMath::Pi(), 0.5}; //rmax, phi_max, ymax
253
254	Double_t rez = mt*ig.Integral(xmin,xmax);
255
256	return rez;
257	}
258
259
260	void AliAnalysisCentralExtrapolate::TsallisFit(){
261	// fits and extrpolates the pt spectrum using TBW model
262
263	printf("---------------------------------------------\n");
264	printf("*** Tsallis Blast Wave Fit ***\n");
265
266	printf("AliAnalysisCentralExtrapolate::Tsallis mass = %g\n", mass);
267
268
269	TStopwatch timer;
270
271	TH1::SetDefaultSumw2();
272	TH1::AddDirectory(0);
273
274	timer.Start();
275
276	TH1D ptSpectra = dynamic_cast<TH1D>(fResultsList->FindObject("hPtESDCorrNorm"));
277
278	TVirtualFitter::SetDefaultFitter("Minuit2");//options Minuit, Minuit2, Fumili, Fumili2
279
280	TF1 *ptFit = new TF1("ptFit",Tsallis,0.2,4.0,3);
281
282	gStyle->SetOptFit(1112);
283
284	ptFit->SetParName(0,"T");
285	ptFit->SetParName(1,"betaS");
286	ptFit->SetParName(2,"q");
287
288	ptFit->SetParameters(0.1,0.5,1.05); //GeV
289	ptFit->SetParLimits(0,0.0,0.5);//GeV
290	ptFit->SetParLimits(1,0.0,1.0);
291	ptFit->SetParLimits(2,1.00001,1.5);
292
293	ptSpectra->Fit("ptFit","Q R B ME N");
294
295	timer.Stop();
296	timer.Print();
297
298	// ----------- print the fit result ----------
299	Double_t chi2 = ptFit->GetChisquare();
300	Double_t ndf = ptFit->GetNDF();
301	Double_t p1 = ptFit->GetParameter(0);
302	Double_t param1Err = ptFit->GetParError(0);
303	Double_t p2 = ptFit->GetParameter(1);
304	Double_t param2Err = ptFit->GetParError(1);
305	Double_t p3 = ptFit->GetParameter(2);
306	Double_t param3Err = ptFit->GetParError(2);
307
308	printf("chi2/ndf = %f/%f\n", chi2,ndf);
309	printf("p1 = %f\tparam1Err = %f\n", p1, param1Err);
310	printf("p2 = %f\tparam2Err = %f\n", p2, param2Err);
311	printf("p3 = %f\tparam3Err = %f\n", p3, param3Err);
312
313
314	// create a new, "extended" histogram
315	TH1F *hPtExtTsallis = new TH1F("PtExtTsallis","Pt Corr Norm Ext",25,0.0,5.0);
316
317	Double_t bin, binerr, test;
318
319	for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){
320
321	bin = ptSpectra->GetBinContent(i);
322	binerr = ptSpectra->GetBinError(i);
323	test = ptSpectra->GetBinLowEdge(i);
324
325	hPtExtTsallis->SetBinContent(i, bin);
326	hPtExtTsallis->SetBinError(i, binerr);
327	}
328
329	Double_t eval;
330	eval = ptFit->Eval(0.1);
331	printf("extrapolated pt value = %g \n", eval);
332	printf("****************************************************\n\n");
333
334	hPtExtTsallis->SetBinContent(1, eval);
335
336	hPtExtTsallis->SetMarkerStyle(24);
337	hPtExtTsallis->SetMarkerColor(kRed); //ESD corrected
338	hPtExtTsallis->GetXaxis()->SetTitle("Pt");
339
340
341	fResultsList->Add(hPtExtTsallis);
342
343	}
344
345	//*********************************************************************************
346	//_________________________________________________________________________________
347	// fit Boltzmann
348
349	Double_t func(Double_t r){
350	//the function under the integral - BGBW
351
352	const Double_t rMax = 11.0;//(fm)
353
354	Double_t betaR = betaS*(r/rMax);
355
356	Double_t rho = TMath::ATanH(betaR);
357
358	Double_t mt = sqrt(massmass + ptpt);// !!! par[0]= pt !!!!!
359
360	Double_t argI0 = (pt*TMath::SinH(rho))/T; //T = par[1]
361
362	Double_t argK1 = (mt*TMath::CosH(rho))/T;
363
364	Double_t i0 = TMath::BesselI0(argI0);
365
366	Double_t k1 = TMath::BesselK1(argK1);
367
368	Double_t f = rmti0*k1;
369
370	return f;
371	}
372
373
374	Double_t Boltzmann(Double_t* const x, Double_t* const par){
375	//computes the integral in the BGBW formula
376
377	pt = x[0];
378	T = par[0];
379
380	betaS = par[1];
381
382
383	ROOT::Math::WrappedFunction<> f1(func);
384
385	ROOT::Math::Integrator ig(f1, ROOT::Math::IntegrationOneDim::kGAUSS,1.E-12,1.E-12,10000);
386
387	Double_t rez = ig.Integral(0.0,11.0);
388
389	return rez;
390	}
391
392
393	void AliAnalysisCentralExtrapolate::BoltzmannFit(){
394	//fits and extrapoates the pt spectrum using the BGBW model
395
396	printf("---------------------------------------------------\n");
397	printf("* Boltzmann-Gibbs Blast Wave Fit *\n");
398
399	printf("AliAnalysisCentralExtrapolate::Boltzmann mass = %g\n", mass);
400
401
402	TStopwatch timer;
403
404	TH1::SetDefaultSumw2();
405	TH1::AddDirectory(0);
406
407	timer.Start();
408
409
410	TH1D ptSpectra = dynamic_cast<TH1D>(fResultsList->FindObject("hPtESDCorrNorm"));
411
412	gStyle->SetOptStat("neRM");
413
414	TVirtualFitter::SetDefaultFitter("Minuit2");
415
416	TF1 *ptFit = new TF1("ptFit",Boltzmann,0.2,4.0,2);
417
418	gStyle->SetOptFit(1112);
419
420	ptFit->SetParName(0,"T");
421	ptFit->SetParName(1,"betaS");
422
423	ptFit->SetParameters(0.1,0.5); //GeV
424	ptFit->SetParLimits(0,0.0,0.5);//GeV
425	ptFit->SetParLimits(1,0.0,1.0);
426
427	ptSpectra->Fit("ptFit","Q R B ME I N");
428
429	timer.Stop();
430	timer.Print();
431
432	// ----------- print the fit results ----------
433	Double_t chi2 = ptFit->GetChisquare();
434	Double_t ndf = ptFit->GetNDF();
435	Double_t p1 = ptFit->GetParameter(0);
436	Double_t param1Err = ptFit->GetParError(0);
437	Double_t p2 = ptFit->GetParameter(1);
438	Double_t param2Err = ptFit->GetParError(1);
439
440	printf("chi2/ndf = %f/%f = %f\n", chi2,ndf,chi2/ndf);
441	printf("p1 = %f (GeV)\tparam1Err = %f\n", p1, param1Err);
442	printf("p2 = %f\tparam2Err = %f\n", p2, param2Err);
443
444
445	TH1F *hPtExtBoltzmann = new TH1F("PtExtBoltzmann","Pt Corr Norm Ext",25,0.0,5.0);
446
447	Double_t bin, binerr, test;
448
449	for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){
450
451	bin = ptSpectra->GetBinContent(i);
452	binerr = ptSpectra->GetBinError(i);
453	test = ptSpectra->GetBinLowEdge(i);
454
455	hPtExtBoltzmann->SetBinContent(i, bin);
456	hPtExtBoltzmann->SetBinError(i, binerr);
457	}
458
459	Double_t eval;
460	eval = ptFit->Eval(0.1);
461	printf("extrapolated pt value = %g \n", eval);
462	printf("********************************************\n\n");
463
464	hPtExtBoltzmann->SetBinContent(1, eval);
465
466	hPtExtBoltzmann->SetMarkerStyle(24);
467	hPtExtBoltzmann->SetMarkerColor(kRed); //ESD corrected
468	hPtExtBoltzmann->GetXaxis()->SetTitle("Pt");
469
470
471	fResultsList->Add(hPtExtBoltzmann);
472
473	}
474
475