[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / MultEvShape / 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("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);

}
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
011e988c	62	Printf("Running %s!\n", name);
011e988c	63
410ff8cc	64	fInputList = new TList();
bde49c8a	65	fResultsList = new TList();
011e988c	66
bde49c8a	67	}
011e988c	68
bde49c8a	69	AliAnalysisCentralExtrapolate::~AliAnalysisCentralExtrapolate() {
011e988c	70
bde49c8a	71	if(fInputList) delete fInputList;
	72	if(fResultsList) delete fResultsList;
	73
	74	}
	75
	76	//************************************************************************************
	77	//____________________________________________________________________________________
	78	// using global variables to avoid the limitations of ROOT::Math::WrappedMultiFunction<>
	79	static Double_t mass = 0.;
	80	static Double_t pt = 0.;
	81	static Double_t T = 0.;
	82	static Double_t betaS = 0.;
	83	static Double_t q = 0.;
	84	static Double_t mt = 0.;
	85
	86	void AliAnalysisCentralExtrapolate::ApplyEff(){
	87	// applies the efficiency map to the pt spectra
	88
	89	TH1::SetDefaultSumw2();
	90	// Int_t stepGen = 0;
	91	Int_t stepRec = 1;
	92
	93	AliCFContainer *data = 0;
410ff8cc	94
410ff8cc	95
bde49c8a	96	TFile *file1 = new TFile("$ALICE_ROOT/PWG2/data/AliAnalysisCentralEfficiency.root", "read");
	97
	98	AliCFEffGrid *eff = 0;
	99
410ff8cc	100	if(fPartType.Contains("kPi")){
410ff8cc	101	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_Pi"));
bde49c8a	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")){
410ff8cc	109	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_K"));
bde49c8a	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")){
410ff8cc	117	data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_P"));
bde49c8a	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");
8a993eac	128	return;
bde49c8a	129	}
	130
	131	if(!eff){
	132	printf("No Eff Grid found! \n");
8a993eac	133	return;
bde49c8a	134	}
410ff8cc	135
410ff8cc	136	// TCanvas *ccorrdata = new TCanvas();
bde49c8a	137	// ccorrdata->Divide(1,2);
bde49c8a	138	// ccorrdata->cd(1);
410ff8cc	139	// ccorrdata->cd(1)->SetLogy();
bde49c8a	140
bde49c8a	141
67a4103d	142	AliCFDataGrid corrdata = new AliCFDataGrid("corrdata","corrected data",data, stepRec);
bde49c8a	143
bde49c8a	144	//correct selection step "reconstructed"
67a4103d	145	// corrdata->SetMeasured(stepRec); //set data to be corrected
bde49c8a	146	corrdata->ApplyEffCorrection(*eff);//apply the correction for efficiency
bde49c8a	147
410ff8cc	148	// TH1D *hPtMC = data->ShowProjection(0,0); //MC distribution ShowProjection(ivar, istep)
	149	// hPtMC->SetMarkerStyle(20);
	150	// hPtMC->SetMarkerColor(kGreen-3);
	151	// hPtMC->GetXaxis()->SetTitle("p_{T}(GeV/c)");
	152	// hPtMC->GetYaxis()->SetTitle("#frac{dN}{p_{T}dp_{T}}");
bde49c8a	153	// hPtMC->Draw("p e1");
410ff8cc	154	//
	155	// TH1D *hPtESDI = corrdata->GetData()->Project(0); //uncorrected ESD Project(ivar)
	156	// hPtESDI->SetMarkerStyle(25);
	157	// hPtESDI->SetMarkerColor(kBlue);
	158	// hPtESDI->GetXaxis()->SetTitle("Pt");
bde49c8a	159	// hPtESDI->Draw("p e1 same");
bde49c8a	160
7faab71b	161	TH1D hPtESDCorr = (TH1D)corrdata->Project(0); //corrected data
bde49c8a	162	hPtESDCorr->SetMarkerStyle(26);
	163	hPtESDCorr->SetMarkerColor(kRed); //ESD corrected
	164	hPtESDCorr->SetName("hPtESDCorr");
	165	// hPtESDCorr->Draw("p e1 same");
	166
	167	// TLegend *leg2 = new TLegend(0.40,0.65,0.87,0.8);
	168	// leg2->AddEntry(hPtMC," generated","p");
	169	// leg2->AddEntry(hPtESDI," reconstructed (not corrected)","p");
	170	// leg2->AddEntry(hPtESDCorr," reconstructed & corrected","p");
	171	// leg2->Draw();
	172
	173
	174	// ccorrdata->cd(2);
	175	// ccorrdata->cd(2)->SetLogy();
410ff8cc	176	//
	177	// TH1D *efficiency = eff->Project(0); //the efficiency vs pt (ipt = 0)
	178	// efficiency->GetXaxis()->SetTitle("p_{T}(GeV/c)");
bde49c8a	179	// efficiency->Draw();
	180
	181
0fb1f0cf	182	TH1D hNoEvt = dynamic_cast<TH1D>(fInputList->FindObject("TaskCentral_NoEvt"));
bde49c8a	183	if(!hNoEvt){
bde49c8a	184	printf("Unable to get the number of events! \n");
8a993eac	185	return;
bde49c8a	186	}
bde49c8a	187
410ff8cc	188
bde49c8a	189	Int_t noEvt = (Int_t)(hNoEvt->GetEntries());
	190
	191	printf("\n No of processed events = %i \n",noEvt);
8a993eac	192
bde49c8a	193	Double_t scale = 1.0/noEvt;
8a993eac	194
bde49c8a	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"));
8a993eac	277	if(!ptSpectra){
410ff8cc	278	printf("TsallisFit: Can't get the normalized spectrum\n");
8a993eac	279	return;
	280	}
	281
	282	if(!ptSpectra->GetEntries()){
410ff8cc	283	printf("TsallisFit: The fit data is empty!\n");
8a993eac	284	return;
8a993eac	285	}
bde49c8a	286
	287	TVirtualFitter::SetDefaultFitter("Minuit2");//options Minuit, Minuit2, Fumili, Fumili2
	288
	289	TF1 *ptFit = new TF1("ptFit",Tsallis,0.2,4.0,3);
	290
	291	gStyle->SetOptFit(1112);
	292
410ff8cc	293	// ptFit->SetParName(0,"T");
	294	// ptFit->SetParName(1,"betaS");
	295	// ptFit->SetParName(2,"q");
bde49c8a	296
	297	ptFit->SetParameters(0.1,0.5,1.05); //GeV
	298	ptFit->SetParLimits(0,0.0,0.5);//GeV
	299	ptFit->SetParLimits(1,0.0,1.0);
	300	ptFit->SetParLimits(2,1.00001,1.5);
	301
	302	ptSpectra->Fit("ptFit","Q R B ME N");
	303
	304	timer.Stop();
	305	timer.Print();
	306
	307	// ----------- print the fit result ----------
	308	Double_t chi2 = ptFit->GetChisquare();
	309	Double_t ndf = ptFit->GetNDF();
	310	Double_t p1 = ptFit->GetParameter(0);
	311	Double_t param1Err = ptFit->GetParError(0);
	312	Double_t p2 = ptFit->GetParameter(1);
	313	Double_t param2Err = ptFit->GetParError(1);
	314	Double_t p3 = ptFit->GetParameter(2);
	315	Double_t param3Err = ptFit->GetParError(2);
	316
	317	printf("chi2/ndf = %f/%f\n", chi2,ndf);
	318	printf("p1 = %f\tparam1Err = %f\n", p1, param1Err);
	319	printf("p2 = %f\tparam2Err = %f\n", p2, param2Err);
	320	printf("p3 = %f\tparam3Err = %f\n", p3, param3Err);
	321
	322
	323	// create a new, "extended" histogram
0fb1f0cf	324	TH1D *hPtExtTsallis = new TH1D("PtExtTsallis","Pt Corr Norm Ext",25,0.0,5.0);
bde49c8a	325
	326	Double_t bin, binerr, test;
	327
	328	for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){
	329
	330	bin = ptSpectra->GetBinContent(i);
	331	binerr = ptSpectra->GetBinError(i);
	332	test = ptSpectra->GetBinLowEdge(i);
	333
	334	hPtExtTsallis->SetBinContent(i, bin);
	335	hPtExtTsallis->SetBinError(i, binerr);
	336	}
	337
	338	Double_t eval;
	339	eval = ptFit->Eval(0.1);
	340	printf("extrapolated pt value = %g \n", eval);
	341	printf("****************************************************\n\n");
	342
	343	hPtExtTsallis->SetBinContent(1, eval);
	344
	345	hPtExtTsallis->SetMarkerStyle(24);
	346	hPtExtTsallis->SetMarkerColor(kRed); //ESD corrected
	347	hPtExtTsallis->GetXaxis()->SetTitle("Pt");
	348
	349
	350	fResultsList->Add(hPtExtTsallis);
	351
	352	}
	353
	354	//*********************************************************************************
	355	//_________________________________________________________________________________
	356	// fit Boltzmann
	357
	358	Double_t func(Double_t r){
	359	//the function under the integral - BGBW
	360
	361	const Double_t rMax = 11.0;//(fm)
	362
	363	Double_t betaR = betaS*(r/rMax);
	364
410ff8cc	365
bde49c8a	366	Double_t rho = TMath::ATanH(betaR);
bde49c8a	367
410ff8cc	368
011e988c	369	mt = sqrt(massmass + ptpt);// !!! par[0]= pt !!!!!
410ff8cc	370
bde49c8a	371	Double_t argI0 = (pt*TMath::SinH(rho))/T; //T = par[1]
410ff8cc	372	if(argI0>700.0) return 0.0; // !! floating point exception protection
bde49c8a	373
	374	Double_t argK1 = (mt*TMath::CosH(rho))/T;
	375
410ff8cc	376
bde49c8a	377	Double_t i0 = TMath::BesselI0(argI0);
bde49c8a	378
410ff8cc	379
bde49c8a	380	Double_t k1 = TMath::BesselK1(argK1);
bde49c8a	381
bde49c8a	382
410ff8cc	383	Double_t f = rmti0*k1;
410ff8cc	384
bde49c8a	385	return f;
	386	}
	387
	388
	389	Double_t Boltzmann(Double_t* const x, Double_t* const par){
	390	//computes the integral in the BGBW formula
	391
	392	pt = x[0];
	393	T = par[0];
	394
	395	betaS = par[1];
	396
	397
	398	ROOT::Math::WrappedFunction<> f1(func);
	399
	400	ROOT::Math::Integrator ig(f1, ROOT::Math::IntegrationOneDim::kGAUSS,1.E-12,1.E-12,10000);
	401
	402	Double_t rez = ig.Integral(0.0,11.0);
	403
	404	return rez;
	405	}
	406
	407
	408	void AliAnalysisCentralExtrapolate::BoltzmannFit(){
	409	//fits and extrapoates the pt spectrum using the BGBW model
	410
	411	printf("---------------------------------------------------\n");
	412	printf("* Boltzmann-Gibbs Blast Wave Fit *\n");
	413
	414	printf("AliAnalysisCentralExtrapolate::Boltzmann mass = %g\n", mass);
	415
	416
	417	TStopwatch timer;
	418
	419	TH1::SetDefaultSumw2();
	420	TH1::AddDirectory(0);
8a993eac	421
bde49c8a	422	timer.Start();
	423
	424
	425	TH1D ptSpectra = dynamic_cast<TH1D>(fResultsList->FindObject("hPtESDCorrNorm"));
8a993eac	426	if(!ptSpectra){
	427	printf("BoltzmannFit: Can't get the normalized spectrum\n");
	428	return;
	429	}
	430
410ff8cc	431	printf("pt spectra get entries: %f\n",ptSpectra->GetEntries());
410ff8cc	432
8a993eac	433	if(!ptSpectra->GetEntries()){
	434	printf("BoltzmannFit: The fit data is empty!\n");
	435	return;
	436	}
bde49c8a	437
	438	gStyle->SetOptStat("neRM");
	439
	440	TVirtualFitter::SetDefaultFitter("Minuit2");
	441
	442	TF1 *ptFit = new TF1("ptFit",Boltzmann,0.2,4.0,2);
	443
	444	gStyle->SetOptFit(1112);
	445
410ff8cc	446	// ptFit->SetParName(0,"T");
410ff8cc	447	// ptFit->SetParName(1,"betaS");
bde49c8a	448
bde49c8a	449	ptFit->SetParameters(0.1,0.5); //GeV
410ff8cc	450	ptFit->SetParLimits(0,0.001,0.5);//GeV
bde49c8a	451	ptFit->SetParLimits(1,0.0,1.0);
410ff8cc	452
bde49c8a	453	ptSpectra->Fit("ptFit","Q R B ME I N");
bde49c8a	454
410ff8cc	455	timer.Stop();
bde49c8a	456	timer.Print();
	457
	458	// ----------- print the fit results ----------
	459	Double_t chi2 = ptFit->GetChisquare();
	460	Double_t ndf = ptFit->GetNDF();
	461	Double_t p1 = ptFit->GetParameter(0);
	462	Double_t param1Err = ptFit->GetParError(0);
	463	Double_t p2 = ptFit->GetParameter(1);
	464	Double_t param2Err = ptFit->GetParError(1);
	465
	466	printf("chi2/ndf = %f/%f = %f\n", chi2,ndf,chi2/ndf);
	467	printf("p1 = %f (GeV)\tparam1Err = %f\n", p1, param1Err);
	468	printf("p2 = %f\tparam2Err = %f\n", p2, param2Err);
	469
	470
0fb1f0cf	471	TH1D *hPtExtBoltzmann = new TH1D("PtExtBoltzmann","Pt Corr Norm Ext",25,0.0,5.0);
bde49c8a	472
	473	Double_t bin, binerr, test;
	474
	475	for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){
	476
	477	bin = ptSpectra->GetBinContent(i);
	478	binerr = ptSpectra->GetBinError(i);
	479	test = ptSpectra->GetBinLowEdge(i);
	480
	481	hPtExtBoltzmann->SetBinContent(i, bin);
	482	hPtExtBoltzmann->SetBinError(i, binerr);
	483	}
	484
	485	Double_t eval;
	486	eval = ptFit->Eval(0.1);
	487	printf("extrapolated pt value = %g \n", eval);
	488	printf("********************************************\n\n");
	489
	490	hPtExtBoltzmann->SetBinContent(1, eval);
	491
	492	hPtExtBoltzmann->SetMarkerStyle(24);
	493	hPtExtBoltzmann->SetMarkerColor(kRed); //ESD corrected
	494	hPtExtBoltzmann->GetXaxis()->SetTitle("Pt");
	495
	496
	497	fResultsList->Add(hPtExtBoltzmann);
	498
	499	}
	500
	501