#pragma link C++ class AliAnalysisCentralCutEvtESD+;
#pragma link C++ class AliAnalysisCentralExtrapolate+;
#pragma link C++ class AliAnalysisTaskCentral+;
+#pragma link C++ class AliBWTools+;
+#pragma link C++ class AliBWFunc+;
+#pragma link C++ class AliLatexTable+;
#endif
--- /dev/null
+void ALICEWorkInProgress(TCanvas *c,TString today="11/05/2010", TString label = "ALICE performance"){
+
+ TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.75,0.82,0.89);
+ myPadLogo->SetFillColor(0);
+ myPadLogo->SetBorderMode(0);
+ myPadLogo->SetBorderSize(2);
+ myPadLogo->SetFrameBorderMode(0);
+ myPadLogo->SetLeftMargin(0.0);
+ myPadLogo->SetTopMargin(0.0);
+ myPadLogo->SetBottomMargin(0.0);
+ myPadLogo->SetRightMargin(0.0);
+ myPadLogo->SetFillStyle(0);
+ myPadLogo->Draw();
+ myPadLogo->cd();
+ TASImage *myAliceLogo = new TASImage("~/WORK/ALICE/ANALYSIS/macros/alice_logo.png");
+ myAliceLogo->Draw();
+ c->cd();
+ TPaveText* t1=new TPaveText(0.65,0.7,0.89,0.75,"NDC");
+ t1->SetFillStyle(0);
+ t1->SetBorderSize(0);
+ t1->AddText(0.,0.,label);
+ t1->SetTextColor(kRed);
+ t1->SetTextFont(42);
+ t1->Draw();
+ TPaveText* t2=new TPaveText(0.65,0.65,0.89,0.7,"NDC");
+ t2->SetFillStyle(0);
+ t2->SetBorderSize(0);
+ t2->SetTextColor(kRed);
+ t2->SetTextFont(52);
+ t2->AddText(0.,0.,today.Data());
+ t2->Draw();
+}
--- /dev/null
+//----------------------------------------------------------------------
+// AliBWFunc
+//
+// This class implements several function useful to fit pt spectra,
+// including but not limited to blast wave models.
+//
+// It can return the same functional for as a function of different
+// variables: dNdpt vs pt, 1/pt dNdpt vs pt, 1/mt dNdmt vs mt.
+//
+// Before getting the function you need, you have to chose the
+// variable you want to use calling AliBWFunc::SetVarType with one of
+// the elements of the VarType_t enum.
+//
+// Warning: not all variables are iplemented for a
+//
+// Author: M. Floris, CERN
+//----------------------------------------------------------------------
+
+#include "AliBWFunc.h"
+#include "TMath.h"
+#include "TF1.h"
+#include "TF3.h"
+#include "TH1.h"
+#include "TSpline.h"
+
+ClassImp(AliBWFunc)
+
+AliBWFunc::AliBWFunc () : fLastFunc(0), fLineWidth(1), fVarType(kdNdpt) {
+
+ // ctor
+ fLineWidth = 1;
+}
+AliBWFunc::~AliBWFunc(){
+
+ // dtor
+ if (fLastFunc) delete fLastFunc;
+
+}
+
+
+TF1 * AliBWFunc::GetHistoFunc(TH1 * h, const char * name) {
+
+ // Regardless of the variable type, this returns a function made
+ // from the histo * a multiplicative normalization.
+
+ fLastFunc = new TF1 (name, StaticHistoFunc, 0.0, 10, 2);
+ fLastFunc->SetParameter(0,1);
+ fLastFunc->FixParameter(1,Double_t(Long64_t(h)));
+ fLastFunc->SetParNames("norm", "pointer to histo");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+
+TF1 * AliBWFunc::GetBGBW(Double_t mass, Double_t beta, Double_t T,
+ Double_t norm, const char * name){
+
+ // Boltzmann-Gibbs blast wave
+
+ switch (fVarType) {
+ case kdNdpt:
+ return GetBGBWdNdptTimesPt(mass,beta,T,norm,name);
+ break;
+ case kOneOverPtdNdpt:
+ return GetBGBWdNdpt(mass,beta,T,norm,name);
+ break;
+ case kOneOverMtdNdmt:
+ AliFatal("Not implemented");
+ break;
+ default:
+ AliFatal("Not implemented");
+ }
+
+ return 0;
+
+}
+
+
+TF1 * AliBWFunc::GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char * name){
+ // Boltzmann
+ switch (fVarType) {
+ case kdNdpt:
+ return GetBoltzmanndNdptTimesPt(mass, T, norm, name);
+ case kOneOverPtdNdpt:
+ AliFatal("Not implemented");
+ break;
+ case kOneOverMtdNdmt:
+ AliFatal("Not implemented");
+ break;
+ default:
+ AliFatal("Not implemented");
+ }
+
+ return 0;
+
+}
+
+
+TF1 * AliBWFunc::GetTsallisBW(Double_t mass, Double_t beta, Double_t T, Double_t q,
+ Double_t norm, Double_t ymax, const char * name){
+
+ // Tsallis blast wave
+ switch (fVarType) {
+ case kdNdpt:
+ return GetTsallisBWdNdptTimesPt(mass,beta,T,q,norm,ymax,name);
+ break;
+ case kOneOverPtdNdpt:
+ return GetTsallisBWdNdpt(mass,beta,T,q,norm,ymax,name);
+ break;
+ case kOneOverMtdNdmt:
+ AliFatal("Not implemented");
+ break;
+ default:
+ AliFatal("Not implemented");
+ }
+
+ return 0;
+
+}
+
+
+TF1 * AliBWFunc::GetMTExp(Double_t mass, Double_t T, Double_t norm, const char * name){
+
+ // Simple exponential in 1/mt*MT
+ switch (fVarType) {
+ case kdNdpt:
+ return GetMTExpdNdptTimesPt(mass,T,norm,name);
+ break;
+ case kOneOverPtdNdpt:
+ return GetMTExpdNdpt(mass,T,norm,name);
+ break;
+ case kOneOverMtdNdmt:
+ AliFatal("Not implemented");
+ break;
+ default:
+ AliFatal("Not implemented");
+ }
+
+ return 0;
+
+
+}
+
+TF1 * AliBWFunc::GetPTExp(Double_t T, Double_t norm, const char * name){
+
+ // Simple exponential in 1/mt*MT
+ switch (fVarType) {
+ case kdNdpt:
+ return GetPTExpdNdptTimesPt(T,norm,name);
+ break;
+ case kOneOverPtdNdpt:
+ AliFatal("Not implemented");
+ break;
+ case kOneOverMtdNdmt:
+ AliFatal("Not implemented");
+ break;
+ default:
+ AliFatal("Not implemented");
+ }
+
+ return 0;
+
+
+}
+
+
+TF1 * AliBWFunc::GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
+ // Levi function (aka Tsallis)
+ switch (fVarType) {
+ case kdNdpt:
+ return GetLevidNdptTimesPt(mass,T,n,norm,name);
+ break;
+ case kOneOverPtdNdpt:
+ return GetLevidNdpt(mass,T,n,norm,name);
+ break;
+ case kOneOverMtdNdmt:
+ return GetLevidNdmt(mass,T,n,norm,name);
+ break;
+ default:
+ AliFatal("Not implemented");
+ }
+
+ return 0;
+
+
+}
+
+TF1 * AliBWFunc::GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name){
+ // power law Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
+ // This is sometimes also called Hagedorn or modified Hagedorn
+
+ switch (fVarType) {
+ case kdNdpt:
+ return GetPowerLawdNdptTimesPt(pt0,n,norm,name);
+ break;
+ case kOneOverPtdNdpt:
+ return GetPowerLawdNdpt(pt0,n,norm,name);
+ break;
+ case kOneOverMtdNdmt:
+ AliFatal("Not Implemented");
+ // return GetUA1dNdmt(mass,T,n,norm,name);
+ break;
+ default:
+ AliFatal("Not implemented");
+ }
+
+ return 0;
+
+
+}
+
+TF1 * AliBWFunc::GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char * name) {
+ // UA1 parametrization Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
+
+ switch (fVarType) {
+ case kdNdpt:
+
+ fLastFunc = new TF1 (name, StaticUA1Func, 0.0, 10, 6);
+ fLastFunc->FixParameter(0,mass);
+ fLastFunc->SetParameter(1,p0star);
+ fLastFunc->SetParameter(2,pt0);
+ fLastFunc->SetParameter(3,n);
+ fLastFunc->SetParameter(4,T);
+ fLastFunc->SetParameter(5,norm);
+ fLastFunc->SetParLimits(1,0.01,1);
+ fLastFunc->SetParLimits(2,0.01,100);
+ fLastFunc->SetParLimits(3,0.01,100);
+ fLastFunc->SetParLimits(4,0.01,100);
+ fLastFunc->SetParNames("mass","p0star","pt0","n","T","norm");
+ fLastFunc->SetNpx(5000);
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+ break;
+ case kOneOverPtdNdpt:
+ AliFatal("Not Implemented");
+ break;
+ case kOneOverMtdNdmt:
+ AliFatal("Not Implemented");
+ // return GetUA1dNdmt(mass,T,n,norm,name);
+ break;
+ default:
+ AliFatal("Not implemented");
+ }
+
+ return 0;
+}
+
+
+
+
+// ________________________________________________________________________
+
+// Backend (private functions and support functions for numerical integration)
+
+Double_t AliBWFunc::StaticHistoFunc(double * x, double* p){
+
+ // provides a function interpolating a histo with a spline;
+ // using double to store a pointer... This is a bad hack. To be replaced
+
+ double norm = p[0];
+
+ TH1 * h = (TH1*) Long64_t(p[1]);
+
+// Int_t bin = h->FindBin(x[0]);
+// double value = h->GetBinContent(bin);
+
+ if (h->FindBin(x[0]) > h->GetNbinsX()) return 0;
+
+ static TH1 * oldptr = 0;
+ TSpline3 * spl = 0;
+ if (h!=oldptr) {
+ spl = new TSpline3(h);
+ }
+ double value = spl->Eval(x[0]);
+ delete spl;
+
+ return value*norm;
+
+}
+
+Double_t AliBWFunc::StaticUA1Func(double * x, double* p) {
+
+
+ // "mass","p0star","pt0","n","T","norm"
+ Double_t mass = p[0];
+ Double_t p0star = p[1];
+ Double_t pt0 = p[2];
+ Double_t n = p[3];
+ Double_t T = p[4];
+ Double_t norm = p[5];
+
+ Double_t xx = x[0];
+
+ static AliBWFunc * self = new AliBWFunc;
+ static TF1 * fPLaw = self->GetPowerLawdNdptTimesPt(pt0, n, norm, "fLocalPLawUA1");
+ static TF1 * fPMTExp = self->GetMTExpdNdptTimesPt (mass, T, norm, "fLocalMTexpUA1");
+
+ fPLaw->SetParameters(norm,pt0,n);
+ fPMTExp->SetParameters(1,T);
+
+
+ Double_t normMT =fPMTExp->Eval(p0star) >0 ? fPLaw->Eval(p0star) / fPMTExp->Eval(p0star) * fPMTExp->GetParameter(0) : 1;
+ fPMTExp->SetParameter(0,normMT);
+
+
+ if (TMath::Abs(fPMTExp->Eval(p0star) - fPLaw->Eval(p0star)) > 0.0001 ) {
+ Printf("AliBWFunc::StaticUA1Func - Wrong norm") ;
+ Printf(" p0* %f NMT: %f N: %f PL: %f MT: %f", p0star, normMT, norm, fPLaw->Eval(p0star), fPMTExp->Eval(p0star));
+ }
+
+ if (xx > p0star) return fPLaw->Eval(xx);
+ return fPMTExp->Eval(xx);
+
+
+
+}
+
+
+Double_t AliBWFunc::IntegrandBG(double * x, double* p){
+ // integrand for boltzman-gibbs blast wave
+
+ double x0 = x[0];
+
+ double mass = p[0];
+ double pT = p[1];
+ double beta = p[2];
+ double T = p[3];
+
+ double mT = TMath::Sqrt(mass*mass+pT*pT);
+
+ double rho0 = TMath::ATanH(beta*x0);
+ double arg0_0 = pT*TMath::SinH(rho0)/T;
+ double arg0_1 = mT*TMath::CosH(rho0)/T;
+ double f0 = x0*mT*TMath::BesselI0(arg0_0)*TMath::BesselK1(arg0_1);
+
+ return f0;
+}
+
+
+
+Double_t AliBWFunc::StaticBGdNdPt(double * x, double* p) {
+
+ // implementation of BGBW (1/pt dNdpt)
+
+ double pT = x[0];;
+
+
+ double mass = p[0];
+ double beta = p[1];
+ double T = p[2];
+
+ static TF1 * fIntBG = 0;
+ if(!fIntBG)
+ fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 4);
+
+ fIntBG->SetParameters(mass, pT, beta, T);
+ double result = fIntBG->Integral(0,1);
+ return result*p[3];//*1e30;;
+
+}
+
+Double_t AliBWFunc::StaticBGdNdPtTimesPt(double * x, double* p) {
+ // BGBW dNdpt implementation
+ return x[0]*StaticBGdNdPt(x,p);
+}
+
+
+TF1 * AliBWFunc::GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t T,
+ Double_t norm, const char * name){
+
+ // BGBW 1/pt dNdpt
+
+ fLastFunc = new TF1 (name, StaticBGdNdPt, 0.0, 10, 4);
+ fLastFunc->SetParameters(mass,beta,T,norm);
+ fLastFunc->SetParNames("mass", "#beta", "T", "norm");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+}
+
+
+//_____________________________________________________________________
+// Tsallis
+
+Double_t AliBWFunc::IntegrandTsallis(double * x, double* p){
+
+ // integrand for numerical integration (tsallis)
+
+ Double_t r = x[0];
+ Double_t phi = x[1];
+ Double_t y = x[2];
+
+ Double_t mass = p[0];
+ Double_t pt = p[1];
+ Double_t beta = p[2];
+ Double_t T = p[3];
+ Double_t q = p[4];
+
+ Double_t mt = TMath::Sqrt(mass*mass+pt*pt);
+
+ Double_t rho = TMath::ATanH(beta*r); // TODO: implement different velocity profiles
+
+ Double_t res = mt*
+ r*TMath::CosH(y) *TMath::Power( (
+ 1+(q-1)/T * (
+ mt*TMath::CosH(y)*TMath::CosH(rho) -
+ pt*TMath::SinH(rho)*TMath::Cos(phi)
+ )
+ ),
+ -1/(q-1)
+ );
+
+
+ return res;
+}
+
+
+
+Double_t AliBWFunc::StaticTsallisdNdPt(double * x, double* p) {
+
+ // tsallis BW implementation 1/pt dNdpt
+
+ double pT = x[0];;
+
+
+ double mass = p[0];
+ double beta = p[1];
+ double T = p[2];
+ double q = p[3];
+
+ Double_t ymax = p[5];
+
+
+ static TF3 * fInt = 0;
+ if(!fInt){
+ fInt = new TF3 ("fIntTsa", IntegrandTsallis, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
+// fInt->SetNpx(10000);
+// fInt->SetNpy(10000);
+// fInt->SetNpz(10000);
+ }
+
+ fInt->SetParameters(mass, pT, beta, T, q);
+ double result = fInt->Integral(0,1, -TMath::Pi(), TMath::Pi(), -ymax, ymax);
+ // double result = fInt->Integral(0,1, -2, 2, -ymax, ymax);
+
+ return result*p[4];//*1e30;;
+
+}
+
+Double_t AliBWFunc::StaticTsallisdNdPtTimesPt(double * x, double* p) {
+
+ // tsallis BW , implementatio of dNdpt
+ return x[0]*StaticTsallisdNdPt(x,p);
+
+}
+
+TF1 * AliBWFunc::GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t T, Double_t q,
+ Double_t norm, Double_t ymax,const char * name){
+
+
+ // tsallis BW, 1/pt dNdpt
+
+ fLastFunc = new TF1 (name, StaticTsallisdNdPt, 0.0, 10, 6);
+ fLastFunc->SetParameters(mass,beta,T,q,norm,ymax);
+ fLastFunc->SetParLimits(1,0.0,0.99);
+ fLastFunc->SetParLimits(2,0.01,0.99);
+ fLastFunc->SetParLimits(3,1.0001,1.9);
+ fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+}
+
+// Times Pt funcs
+// Boltzmann-Gibbs Blast Wave
+TF1 * AliBWFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T,
+ Double_t norm, const char * name){
+
+ // BGBW, dNdpt
+
+ fLastFunc = new TF1 (name, StaticBGdNdPtTimesPt, 0.0, 10, 4);
+ fLastFunc->SetParameters(mass,beta,T,norm);
+ fLastFunc->SetParNames("mass", "#beta", "T", "norm");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+
+
+TF1 * AliBWFunc::GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T, Double_t q,
+ Double_t norm, Double_t ymax, const char * name){
+
+// Tsallis blast wave, dNdpt
+
+ fLastFunc = new TF1 (name, StaticTsallisdNdPtTimesPt, 0.0, 10, 6);
+ fLastFunc->SetParameters(mass,beta,T,q,norm,ymax);
+ fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+
+
+TF1 * AliBWFunc::GetMTExpdNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name){
+
+ // Simple exponential in 1/mt*MT, as a function of dNdpt
+ char formula[500];
+ sprintf(formula,"[0]*x*exp(-sqrt(x**2+%f**2)/[1])", mass);
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, T);
+ fLastFunc->SetParLimits(1, 0.01, 10);
+ fLastFunc->SetParNames("norm", "T");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+
+TF1 * AliBWFunc::GetPTExpdNdptTimesPt(Double_t T, Double_t norm, const char * name){
+
+ // Simple exponential in 1/pt*dNdpT, as a function of dNdpt
+ char formula[500];
+ sprintf(formula,"[0]*x*exp(-x/[1])");
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, T);
+ fLastFunc->SetParLimits(1, 0.01, 10);
+ fLastFunc->SetParNames("norm", "T");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+
+TF1 * AliBWFunc::GetBoltzmanndNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name){
+ // Boltzmann (exp in 1/mt*dNdmT times mt) as a function of dNdpt
+ char formula[500];
+ sprintf(formula,"[0]*x*sqrt(x**2+%f**2)*exp(-sqrt(x**2+%f**2)/[1])", mass,mass);
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, T);
+ fLastFunc->SetParLimits(1, 0.01, 10);
+ fLastFunc->SetParNames("norm", "T");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+
+// Tsallis (no BW, a la CMS)
+// TF1 * AliBWFunc::GetTsallisdNdptTimesPt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name){
+
+// char formula[500];
+// // sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2)-%f)),(-1/([2]-1)))", mass, mass); //CMS
+// sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass); // STAR
+// //sprintf(formula,"[0]*x*sqrt(x**2+%f**2)*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass,mass); // STAR * mt
+// fLastFunc=new TF1(name,formula,0,10);
+// fLastFunc->SetParameters(norm, T, q);
+// fLastFunc->SetParLimits(1, 0.001, 10);
+// fLastFunc->SetParNames("norm", "T", "q");
+// fLastFunc->SetLineWidth(fLineWidth);
+// return fLastFunc;
+
+
+// }
+
+
+TF1 * AliBWFunc::GetLevidNdptTimesPt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
+
+ // Levi function, dNdpt
+ char formula[500];
+
+ sprintf(formula,"( x*[0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2]) )^(-[1])");
+ // sprintf(formula,"( x*[0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (sqrt([3]*[3]+x*x))/([1]*[2]) )^(-[1])");
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, n, T,mass);
+ fLastFunc->SetParLimits(2, 0.01, 10);
+ fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
+ fLastFunc->FixParameter(3,mass);
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+TF1 * AliBWFunc::GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name){
+
+ // PowerLaw function, dNdpt
+ char formula[500];
+
+ sprintf(formula,"x*[0]*( 1 + x/[1] )^(-[2])");
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, pt0, n);
+ fLastFunc->SetParLimits(1, 0.01, 10);
+ //fLastFunc->SetParLimits(2, 0.01, 50);
+ fLastFunc->SetParNames("norm", "pt0", "n");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+TF1 * AliBWFunc::GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name){
+
+ // PowerLaw function, 1/pt dNdpt
+ char formula[500];
+
+ sprintf(formula," [0]*( 1 + x/[1] )^(-[2])");
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, pt0, n);
+ // fLastFunc->SetParLimits(2, 0.01, 10);
+ fLastFunc->SetParNames("norm", "pt0", "n");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+
+TF1 * AliBWFunc::GetLevidNdpt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
+
+ // Levi function, dNdpt
+ char formula[500];
+
+ sprintf(formula,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2]) )^(-[1])");
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, n, T,mass);
+ fLastFunc->SetParLimits(2, 0.01, 10);
+ fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
+ fLastFunc->FixParameter(3,mass);
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+TF1 * AliBWFunc::GetLevidNdmt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
+
+ // Levi function, dNdmt
+ char formula[500];
+
+ // sprintf(formula,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x -[3])/([1]*[2]) )^(-[1])");
+ sprintf(formula,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + x/([1]*[2]) )^(-[1])");
+ // sprintf(formula,"[0] * ( 1 + x/([1]*[2]) )^(-[1])");
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, n, T,mass);
+ fLastFunc->SetParLimits(2, 0.01, 10);
+ fLastFunc->SetParNames("norm", "n", "T", "mass");
+ fLastFunc->FixParameter(3,mass);
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+
+}
+
+
+
+// Test Function
+Double_t AliBWFunc::IntegrandTest(double * x, double* p){
+
+ // test function
+
+ Double_t y = x[0];
+
+ Double_t mass = p[0];
+ Double_t pt = p[1];
+ Double_t T = p[2];
+
+ Double_t mt = TMath::Sqrt(mass*mass+pt*pt);
+
+ return mt*TMath::CosH(y)*TMath::Exp(-mt*TMath::CosH(y)/T);
+
+}
+
+Double_t AliBWFunc::StaticTest(double * x, double* p) {
+
+ // test function
+
+ double pT = x[0];;
+
+
+ double mass = p[0];
+ double T = p[1];
+ Double_t ymax = p[3];
+
+
+ static TF3 * fIntTest = 0;
+ if(!fIntTest){
+ fIntTest = new TF3 ("fIntTest", IntegrandTest, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
+ // fInt->SetNpx(10000);
+ }
+
+ fIntTest->SetParameters(mass, pT, T);
+ double result = fIntTest->Integral(-ymax, ymax);
+
+ return result*p[2];//*1e30;;
+
+}
+
+TF1 * AliBWFunc::GetTestFunc(Double_t mass, Double_t T, Double_t norm, Double_t ymax, const char * name){
+
+ // test function
+
+ fLastFunc = new TF1 (name, StaticTest, 0.0, 10, 4);
+ fLastFunc->SetParameters(mass,T,norm,ymax);
+ fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+
+}
+
+
+//___________________________________________________________
+
+
+TF1 * AliBWFunc::GetMTExpdNdpt(Double_t mass, Double_t T, Double_t norm, const char * name){
+ // Simple exp in 1/mt dNdmt, as a function of dNdpt
+ // mt scaling
+ char formula[500];
+ sprintf(formula,"[0]*exp(-sqrt(x**2+%f**2)/[1])", mass);
+ fLastFunc=new TF1(name,formula,0,10);
+ fLastFunc->SetParameters(norm, T);
+ fLastFunc->SetParLimits(1, 0.01, 10);
+ fLastFunc->SetParNames("norm", "T");
+ fLastFunc->SetLineWidth(fLineWidth);
+ return fLastFunc;
+}
+
+
+// // Simple tsallis (a la CMS)
+// TF1 * AliBWFunc::GetTsallisdNdpt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name){
+
+// char formula[500];
+// sprintf(formula,"[0]*sqrt(x**2+%f**2)*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass,mass);
+// fLastFunc=new TF1(name,formula,0,10);
+// fLastFunc->SetParameters(norm, T, q);
+// fLastFunc->SetParLimits(1, 0.01, 10);
+// fLastFunc->SetParNames("norm", "T", "q");
+// fLastFunc->SetLineWidth(fLineWidth);
+// return fLastFunc;
+// }
--- /dev/null
+//----------------------------------------------------------------------
+// AliBWFunc
+//
+// This class implements several function useful to fit pt spectra,
+// including but not limited to blast wave models.
+//
+// It can return the same functional for as a function of different
+// variables: dNdpt vs pt, 1/pt dNdpt vs pt, 1/mt dNdmt vs mt.
+//
+// Before getting the function you need, you have to chose the
+// variable you want to use calling AliBWFunc::SetVarType with one of
+// the elements of the VarType_t enum.
+//
+// Author: M. Floris, CERN
+//----------------------------------------------------------------------
+
+#ifndef ALIBWFUNC_H
+#define ALIBWFUNC_H
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TObject.h"
+#include "AliLog.h"
+
+class TF1;
+class TH1;
+
+#endif
+
+
+
+class AliBWFunc : public TObject {
+
+
+public:
+ // define the variables used for the function
+ typedef enum {kdNdpt,kOneOverPtdNdpt,kOneOverMtdNdmt} VarType_t;
+
+ AliBWFunc();
+ ~AliBWFunc();
+
+ // Boltzmann-Gibbs blast wave
+ TF1 * GetBGBW(Double_t mass, Double_t beta, Double_t T,
+ Double_t norm, const char * name = "fBGBW");
+
+ // Boltzmann
+ TF1 * GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char * name ="fBoltzmann");
+
+ // Tsallis blast wave
+ TF1 * GetTsallisBW(Double_t mass, Double_t beta, Double_t T, Double_t q,
+ Double_t norm, Double_t ymax = 0.5, const char * name = "fTsallisBW");
+
+ // Simple exponential in 1/mt*dNdmt
+ TF1 * GetMTExp(Double_t mass, Double_t T, Double_t norm, const char * name ="fExp");
+
+ // Simple exponential in 1/pt*dNdpt
+ TF1 * GetPTExp(Double_t T, Double_t norm, const char * name ="fExp");
+
+ // Tsallis (no BW, a la CMS)
+ TF1 * GetTsallis(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name="fTsallis")
+ {return GetLevi (mass,T,1/(q-1),norm,name);}
+
+ // Levi function (aka Tsallis)
+ TF1 * GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name="fLevi");
+
+ // UA1 function
+ TF1 * GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char * name="fUA1");
+
+ // Function derived from a histo
+ TF1 * GetHistoFunc(TH1 * h, const char * name = "fHisto");
+
+ // Power law
+ TF1 * GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name="fPowerLaw");
+
+
+ void SetVarType(VarType_t tp) {fVarType=tp;}
+
+protected:
+
+ // dNdpt here means 1/pt dN/dpt
+
+ // Boltzmann-Gibbs Blast Wave
+ TF1 * GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t T,
+ Double_t norm, const char * name = "fBGBW");
+
+ // Tsallis blast wave
+ TF1 * GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t T, Double_t q,
+ Double_t norm, Double_t ymax = 0.5, const char * name = "fTsallisBW");
+
+ // Simple exponential in 1/mt*MT
+ TF1 * GetMTExpdNdpt(Double_t mass, Double_t T, Double_t norm, const char * name ="fExp");
+
+ // Tsallis (no BW, a la CMS)
+ TF1 * GetTsallisdNdpt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name="fTsallis");
+
+ // Levi function
+ TF1 * GetLevidNdpt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name="fLevi");
+
+ // Power Law function
+ TF1 * GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name="fLevi");
+
+ // UA1 function
+ TF1 * GetUA1dNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name="fLevi");
+
+ // TimesPt means dNdpt
+
+ // Boltzmann-Gibbs Blast Wave
+ TF1 * GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T,
+ Double_t norm, const char * name = "fBGBWTimesPt");
+
+ // Tsallis blast wave
+ TF1 * GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T, Double_t q,
+ Double_t norm, Double_t ymax = 0.5, const char * name = "fTsallisBWTimesPt");
+ // Levi function
+ TF1 * GetLevidNdptTimesPt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name="fLevi");
+
+ // UA1 function
+ TF1 * GetUA1dNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name="fLevi");
+
+ // PowerLaw function
+ TF1 * GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name="fLevi");
+
+ // Simple exponential in 1/mt*dNdmT
+ TF1 * GetMTExpdNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name ="fMtExpTimesPt");
+
+ // Simple exponential in 1/mp*dNdpT
+ TF1 * GetPTExpdNdptTimesPt(Double_t T, Double_t norm, const char * name ="fPtExpTimesPt");
+
+ // Boltzmann (exp in 1/mt*dNdmT times mt)
+ TF1 * GetBoltzmanndNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name ="fBoltzmannTimesPt");
+
+ // Tsallis (no BW, a la CMS)
+ TF1 * GetTsallisdNdptTimesPt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name="fTsallisTimesPt");
+
+ // 1/mt dNdmt
+
+ // Levi
+ TF1 * GetLevidNdmt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name = "fLeviMt");
+
+
+ // gereral setters
+ void SetLineWidth(Width_t width) { fLineWidth = width;}
+
+ TF1 * GetTestFunc(Double_t mass, Double_t T, Double_t norm, Double_t ymax, const char * name ="fTest") ;
+
+ // static functions for TF1
+ // Boltzmann-Gibbs Blast Wave
+ static Double_t StaticBGdNdPt(double * x, double* p);
+ static Double_t StaticBGdNdPtTimesPt(double * x, double* p);
+ // Tsallis blast wave
+ static Double_t StaticTsallisdNdPt(double * x, double* p);
+ static Double_t StaticTsallisdNdPtTimesPt(double * x, double* p);
+ // Helper funcs for numeric integration
+ static Double_t IntegrandBG(double * x, double* p);
+ static Double_t IntegrandTsallis(double * x, double* p);
+
+ // Test func
+ static Double_t StaticTest(double * x, double* p);
+ static Double_t IntegrandTest(double * x, double* p);
+
+ // histo func
+ static Double_t StaticHistoFunc(double * x, double* p);
+
+ // UA1 parametrization
+ static Double_t StaticUA1Func(double * x, double* p);
+
+
+private:
+
+
+ TF1 * fLastFunc;
+ Width_t fLineWidth;
+ VarType_t fVarType;
+
+ AliBWFunc(const AliBWFunc&); // not implemented
+ AliBWFunc& operator=(const AliBWFunc&); // not implemented
+
+
+ ClassDef(AliBWFunc, 1)
+
+
+};
+
+#endif
--- /dev/null
+// ----------------------------------------------------------------------
+// AliBWTools
+//
+// This class provides some tools which can be useful in the analsis
+// of spectra, to fit or transform histograms. See the comments of the
+// individual methods for details
+//
+// Author: M. Floris (CERN)
+// ----------------------------------------------------------------------
+
+#include "AliBWTools.h"
+#include "TH1D.h"
+#include "TF1.h"
+#include "TH1.h"
+#include "TMath.h"
+#include "TGraphErrors.h"
+#include "TVirtualFitter.h"
+#include "TMinuit.h"
+#include "AliLog.h"
+#include <iostream>
+
+using namespace std;
+
+TF1 * AliBWTools::fFuncForNormalized = 0;
+
+ClassImp(AliBWTools)
+
+AliBWTools::AliBWTools() {
+ // ctor
+}
+
+AliBWTools::~AliBWTools(){
+ // dtor
+}
+
+TH1 * AliBWTools::GetdNdmtFromdNdpt(TH1 * hpt, Double_t mass) {
+ // convert the x axis from pt to mt. Assumes you have 1/pt dNdpt in the histo you start with
+
+ Int_t nbins = hpt->GetNbinsX();
+ Float_t * xbins = new Float_t[nbins+1];
+ for(Int_t ibins = 0; ibins <= nbins; ibins++){
+ xbins[ibins] = TMath::Sqrt(hpt->GetBinLowEdge(ibins+1)*hpt->GetBinLowEdge(ibins+1) +
+ mass *mass) - mass;
+// xbins[ibins] = TMath::Sqrt(hpt->GetBinLowEdge(ibins+1)*hpt->GetBinLowEdge(ibins+1) +
+// mass *mass);
+ // cout << ibins << " "<< xbins[ibins] << endl;
+
+ }
+
+ TH1D * hmt = new TH1D(TString(hpt->GetName())+"_mt",
+ TString(hpt->GetName())+"_mt",
+ nbins, xbins);
+ for(Int_t ibins = 1; ibins <= nbins; ibins++){
+ hmt->SetBinContent(ibins, hpt->GetBinContent(ibins));
+ hmt->SetBinError(ibins, hpt->GetBinError(ibins));
+
+ }
+
+ hmt->SetXTitle("m_{t} - m_{0} (GeV/c^{2})");
+ hmt->SetYTitle("1/m_{t} dN/dm_{t} (a.u.)");
+
+ return hmt;
+
+}
+
+TH1 * AliBWTools::GetdNdPtFromOneOverPt(TH1 * h1Pt) {
+
+ // convert an histo from 1/pt dNdpt to dNdpt, using the pt at the center of the bin
+
+
+ TH1 * hPt = (TH1 *) h1Pt->Clone((TString(h1Pt->GetName()) + "_inv").Data());
+ hPt->Reset();
+
+ Int_t nbinx = hPt->GetNbinsX();
+
+ for(Int_t ibinx = 1; ibinx <= nbinx; ibinx++){
+
+ Double_t cont = h1Pt->GetBinContent(ibinx);
+ Double_t err = h1Pt->GetBinError(ibinx);
+
+ Double_t pt = h1Pt->GetBinCenter(ibinx);
+
+ if(pt != 0) {
+ cont *= pt;
+ err *= pt;
+ } else {
+ cont = 0;
+ err = 0;
+ }
+
+ hPt->SetBinContent(ibinx, cont);
+ hPt->SetBinError(ibinx, err);
+
+ }
+
+ hPt->SetXTitle("p_{t} (GeV)");
+ hPt->SetYTitle("dN/dp_{t} (GeV^{-2})");
+
+ return hPt;
+
+}
+
+
+
+
+TH1 * AliBWTools::GetOneOverPtdNdPt(TH1 * hPt) {
+
+ // convert an histo from dNdpt to 1/pt dNdpt, using the pt at the center of the bin
+
+ TH1 * h1Pt = (TH1 *) hPt->Clone((TString(hPt->GetName()) + "_inv").Data());
+ h1Pt->Reset();
+
+ Int_t nbinx = h1Pt->GetNbinsX();
+
+ for(Int_t ibinx = 1; ibinx <= nbinx; ibinx++){
+
+ Double_t cont = hPt->GetBinContent(ibinx);
+ Double_t err = hPt->GetBinError(ibinx);
+
+ Double_t pt = hPt->GetBinCenter(ibinx);
+
+ if(pt != 0) {
+ cont /= pt;
+ err /= pt;
+ } else {
+ cont = 0;
+ err = 0;
+ }
+
+ h1Pt->SetBinContent(ibinx, cont);
+ h1Pt->SetBinError(ibinx, err);
+
+ }
+
+ h1Pt->SetXTitle("p_{t} (GeV)");
+ h1Pt->SetYTitle("1/p_{t} dN/dp_{t} (GeV^{-2})");
+
+ return h1Pt;
+
+}
+
+
+TGraphErrors * AliBWTools::GetGraphFromHisto(TH1F * h, Bool_t binWidth) {
+ // Convert a histo to a graph
+ // if binWidth is true ex is set to the bin width of the histos, otherwise it is set to zero
+ Int_t nbin = h->GetNbinsX();
+
+ TGraphErrors * g = new TGraphErrors();
+ Int_t ipoint = 0;
+ for(Int_t ibin = 1; ibin <= nbin; ibin++){
+ Double_t xerr = binWidth ? h->GetBinWidth(ibin)/2 : 0;
+ if (h->GetBinContent(ibin)!=0) {
+ g->SetPoint (ipoint, h->GetBinCenter(ibin), h->GetBinContent(ibin));
+ g->SetPointError(ipoint, xerr, h->GetBinError(ibin));
+ ipoint++;
+ }
+ }
+
+ g->SetMarkerStyle(h->GetMarkerStyle());
+ g->SetMarkerColor(h->GetMarkerColor());
+ g->SetLineColor(h->GetLineColor());
+
+ g->SetTitle(h->GetTitle());
+ g->SetName(TString("g_")+h->GetName());
+
+ return g;
+
+}
+
+TH1F * AliBWTools::GetHistoFromGraph(TGraphErrors * g, TH1F* hTemplate) {
+
+ // convert a graph to histo with the binning of hTemplate.
+ // Warning: the template should be chosen
+ // properly: if you have two graph points in the same histo bin this
+ // won't work!
+
+ TH1F * h = (TH1F*) hTemplate->Clone(TString("h_")+g->GetName());
+ h->Reset();
+ Int_t npoint = g->GetN();
+ // g->Print();
+ for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+ Float_t x = g->GetX() [ipoint];
+ Float_t y = g->GetY() [ipoint];
+ Float_t ey = g->GetEY()[ipoint];
+ Int_t bin = h->FindBin(x);
+ // cout << "bin: "<< bin << " -> " << x << ", "<< y <<", " << ey << endl;
+
+ h->SetBinContent(bin,y);
+ h->SetBinError (bin,ey);
+ }
+
+ h->SetMarkerStyle(g->GetMarkerStyle());
+ h->SetMarkerColor(g->GetMarkerColor());
+ h->SetLineColor (g->GetLineColor());
+
+
+ return h;
+}
+
+TGraphErrors * AliBWTools::ConcatenateGraphs(TGraphErrors * g1,TGraphErrors * g2){
+
+ // concatenates two graphs
+
+ Int_t npoint1=g1->GetN();
+ Int_t npoint2=g2->GetN();
+
+ TGraphErrors * gClone = (TGraphErrors*) g1->Clone();
+
+// for(Int_t ipoint = 0; ipoint < npoint1; ipoint++){
+// gClone->SetPointError(ipoint,0,g1->GetEY()[ipoint]);
+
+// }
+ for(Int_t ipoint = 0; ipoint < npoint2; ipoint++){
+ gClone->SetPoint(ipoint+npoint1,g2->GetX()[ipoint],g2->GetY()[ipoint]);
+ gClone->SetPointError(ipoint+npoint1,g2->GetEX()[ipoint],g2->GetEY()[ipoint]);
+ // gClone->SetPointError(ipoint+npoint1,0,g2->GetEY()[ipoint]);
+ }
+
+ gClone->GetHistogram()->GetXaxis()->SetTimeDisplay(1);
+ gClone->SetTitle(TString(gClone->GetTitle())+" + "+g2->GetTitle());
+ gClone->SetName(TString(gClone->GetName())+"_"+g2->GetName());
+
+ return gClone;
+}
+
+
+TH1F * AliBWTools::CombineHistos(TH1 * h1, TH1 * h2, TH1* htemplate, Float_t renorm1){
+
+ // Combine two histos. This assumes the histos have the same binning
+ // in the overlapping region. It computes the arithmetic mean in the
+ // overlapping region and assigns as an error the relative error
+ // h2. TO BE IMPROVED
+
+ TH1F * hcomb = (TH1F*) htemplate->Clone(TString(h1->GetName())+"_"+h2->GetName());
+
+ TH1F * h1local = (TH1F*) h1->Clone();
+ h1local->Scale(renorm1);
+
+ Int_t nBinComb = hcomb->GetNbinsX();
+ for(Int_t ibin = 1; ibin <= nBinComb; ibin++){
+ Int_t ibin1 = h1local->FindBin(hcomb->GetBinCenter(ibin));
+ Int_t ibin2 = h2->FindBin(hcomb->GetBinCenter(ibin));
+
+ if (h1local->GetBinContent(ibin1) == 0 && h2->GetBinContent(ibin2) == 0) {
+ // None has data: go to next bin
+ hcomb->SetBinContent(ibin,0);
+ hcomb->SetBinError (ibin,0);
+ } else if(h1local->GetBinContent(ibin1) != 0 && h2->GetBinContent(ibin2) == 0) {
+ // take data from h1local:
+ hcomb->SetBinContent(ibin,h1local->GetBinContent(ibin1));
+ hcomb->SetBinError (ibin,h1local->GetBinError(ibin1));
+ } else if(h1local->GetBinContent(ibin1) == 0 && h2->GetBinContent(ibin2) != 0) {
+ // take data from h2:
+ hcomb->SetBinContent(ibin,h2->GetBinContent(ibin2));
+ hcomb->SetBinError (ibin,h2->GetBinError(ibin2));
+ } else {
+ hcomb->SetBinContent(ibin,(h1local->GetBinContent(ibin1) +h2->GetBinContent(ibin2))/2);
+ // hcomb->SetBinError (ibin,h1local->GetBinError(ibin1)/h1local->GetBinContent(ibin1)*hcomb->GetBinContent(ibin));
+ hcomb->SetBinError (ibin,h2->GetBinError(ibin2)/h2->GetBinContent(ibin2)*hcomb->GetBinContent(ibin));
+ }
+
+
+ }
+
+
+ hcomb->SetMarkerStyle(h1local->GetMarkerStyle());
+ hcomb->SetMarkerColor(h1local->GetMarkerColor());
+ hcomb->SetLineColor (h1local->GetLineColor());
+
+ hcomb->SetXTitle(h1local->GetXaxis()->GetTitle());
+ hcomb->SetYTitle(h1local->GetYaxis()->GetTitle());
+ delete h1local;
+ return hcomb;
+
+}
+
+
+void AliBWTools::GetFromHistoGraphDifferentX(TH1F * h, TF1 * f, TGraphErrors ** gBarycentre, TGraphErrors ** gXlw) {
+
+ // Computes the baycentre in each bin and the xlw as defined in NIMA
+ // 355 - 541 using f. Returs 2 graphs with the same y content of h
+ // but with a different x (barycentre and xlw)
+
+ Int_t nbin = h->GetNbinsX();
+
+ (*gBarycentre) = new TGraphErrors();
+ (*gXlw) = new TGraphErrors();
+
+ Int_t ipoint = 0;
+ for(Int_t ibin = 1; ibin <= nbin; ibin++){
+ Float_t min = h->GetBinLowEdge(ibin);
+ Float_t max = h->GetBinLowEdge(ibin+1);
+ Double_t xerr = 0;
+ Double_t xbar = f->Mean(min,max);
+ // compute x_LW
+ Double_t temp = 1./(max-min) * f->Integral(min,max);
+ Double_t epsilon = 0.000000001;
+ Double_t increment = 0.0000000001;
+ Double_t x_lw = min;
+
+ while ((f->Eval(x_lw)- temp) > epsilon) {
+ x_lw += increment;
+ if(x_lw > max) {
+ Printf("Cannot find x_lw");
+ break;
+ }
+ }
+
+ if (h->GetBinContent(ibin)!=0) {
+ (*gBarycentre)->SetPoint (ipoint, xbar, h->GetBinContent(ibin));
+ (*gBarycentre)->SetPointError(ipoint, xerr, h->GetBinError(ibin));
+ (*gXlw) ->SetPoint (ipoint, x_lw, h->GetBinContent(ibin));
+ (*gXlw) ->SetPointError(ipoint, xerr, h->GetBinError(ibin));
+ ipoint++;
+ }
+ }
+
+ (*gBarycentre)->SetMarkerStyle(h->GetMarkerStyle());
+ (*gBarycentre)->SetMarkerColor(h->GetMarkerColor());
+ (*gBarycentre)->SetLineColor(h->GetLineColor());
+
+ (*gBarycentre)->SetTitle(h->GetTitle());
+ (*gBarycentre)->SetName(TString("g_")+h->GetName());
+
+ (*gXlw)->SetMarkerStyle(h->GetMarkerStyle());
+ (*gXlw)->SetMarkerColor(h->GetMarkerColor());
+ (*gXlw)->SetLineColor(h->GetLineColor());
+ (*gXlw)->SetTitle(h->GetTitle());
+ (*gXlw)->SetName(TString("g_")+h->GetName());
+
+
+}
+
+
+Float_t AliBWTools::GetMean(TH1F * h, Float_t min, Float_t max) {
+
+ // Get the mean of histo in a range; root is not reliable in sub ranges with variable binning.
+
+ Int_t minBin = h->FindBin(min);
+ Int_t maxBin = h->FindBin(max);
+
+ Float_t mean = 0 ;
+ Float_t integral = 0;
+ for(Int_t ibin = minBin; ibin < maxBin; ibin++){
+ mean = mean + (h->GetBinCenter(ibin) * h->GetBinWidth(ibin)* h->GetBinContent(ibin));
+ integral = integral + h->GetBinContent(ibin) * h->GetBinWidth(ibin);
+ }
+
+ return mean/integral;
+
+
+}
+
+void AliBWTools::GetMean(TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) {
+
+ // Get the mean of function in a range;
+
+ return GetMoment("fMean", "x*", func, mean, error, min,max);
+
+}
+
+void AliBWTools::GetMeanSquare(TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) {
+
+ // Get the mean2 of function in a range;
+
+ return GetMoment("fMean2", "x*x*", func, mean, error, min,max);
+
+
+}
+
+void AliBWTools::GetMoment(TString name, TString var, TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) {
+
+ // returns the integral of a function derived from func by prefixing some operation.
+ // Used as a base method for mean and mean2
+ Printf("AliBWTools::GetMoment: Error on <pt> is not correct!!! It is overestimated, fix required");
+ Int_t npar = func->GetNpar();
+
+ TF1 * f = new TF1(name, var+func->GetName()); // FIXME
+// fFuncForNormalized = func;// TMP: computing mean pt
+// TF1 * f = new TF1(name,GetNormalizedFunc,0,10,npar);// FIXME
+// for(Int_t ipar = 0; ipar < npar; ipar++){ // FIXME
+// f->SetParameter(ipar,func->GetParameter(ipar)); // FIXME
+// }
+
+
+ // The parameter of the function used to compute the mean should be
+ // the same as the parent function: fixed if needed and they should
+ // also have the same errors.
+
+ // cout << "npar :" << npar << endl;
+
+ for(Int_t ipar = 0; ipar < npar; ipar++){
+ Double_t parmin, parmax;
+ Double_t value = func->GetParameter(ipar);
+ func->GetParLimits(ipar, parmin, parmax);
+ if ( parmin == parmax ) {
+ if ( parmin != 0 || (parmin == 1 && value == 0) ) {
+ f->FixParameter(ipar,func->GetParameter(ipar));
+ // cout << "Fixing " << ipar << "("<<value<<","<<parmin<<","<<parmax<<")"<<endl;
+ }
+ else {
+ f->SetParError (ipar,func->GetParError(ipar) );
+ // cout << "Setting Err" << ipar << "("<<func->GetParError(ipar)<<")"<<endl;
+ }
+ }
+ else {
+ f->SetParError (ipar,func->GetParError(ipar) );
+ // cout << "Setting Err" << ipar << "("<<func->GetParError(ipar)<<")"<<endl;
+ }
+ }
+ // f->Print();
+// mean = f->Integral (min,max)/func->Integral(min,max);
+// error = f->IntegralError(min,max)/func->Integral(min,max);
+ mean = f->Integral (min,max);
+ error = f->IntegralError(min,max);
+// cout << "Mean " << mean <<"+-"<< error<< endl;
+// cout << "Integral Error " << error << endl;
+
+}
+
+Double_t AliBWTools::GetNormalizedFunc(double * x, double* p){
+
+ // Static function used to provide normalized pointer to a function
+
+ Int_t npar = fFuncForNormalized->GetNpar();
+ for(Int_t ipar = 0; ipar < npar; ipar++){ // FIXME
+ fFuncForNormalized->SetParameter(ipar,p[ipar]); // FIXME
+ }
+
+ return x[0]*fFuncForNormalized->Eval(x[0])/fFuncForNormalized->Integral(0,100);
+
+}
+
+
+Bool_t AliBWTools::Fit (TH1 * h1, TF1* func, Float_t min, Float_t max) {
+
+ // Fits h1 with func, preforms several checks on the quality of the
+ // fit and tries to improve it. If the fit is not good enough, it
+ // returs false.
+
+ Double_t amin; Double_t edm; Double_t errdef; Int_t nvpar; Int_t nparx;
+ TVirtualFitter *fitter;
+ cout << "--- Fitting : " << h1->GetName() << " ["<< h1->GetTitle() <<"] ---"<< endl;
+
+ h1-> Fit(func,"IME0","",min,max);
+ Int_t fitResult = h1-> Fit(func,"IME0","",min,max);
+// h1-> Fit(func,"0","",min,max);
+// Int_t fitResult = h1-> Fit(func,"0IE","",min,max);
+
+
+// From TH1:
+// The fitStatus is 0 if the fit is OK (i.e no error occurred). The
+// value of the fit status code is negative in case of an error not
+// connected with the minimization procedure, for example when a wrong
+// function is used. Otherwise the return value is the one returned
+// from the minimization procedure. When TMinuit (default case) or
+// Minuit2 are used as minimizer the status returned is : fitStatus =
+// migradResult + 10*minosResult + 100*hesseResult +
+// 1000*improveResult. TMinuit will return 0 (for migrad, minos,
+// hesse or improve) in case of success and 4 in case of error (see
+// the documentation of TMinuit::mnexcm). So for example, for an error
+// only in Minos but not in Migrad a fitStatus of 40 will be returned.
+// Minuit2 will return also 0 in case of success and different values
+// in migrad minos or hesse depending on the error. See in this case
+// the documentation of Minuit2Minimizer::Minimize for the
+// migradResult, Minuit2Minimizer::GetMinosError for the minosResult
+// and Minuit2Minimizer::Hesse for the hesseResult. If other
+// minimizers are used see their specific documentation for the status
+// code returned. For example in the case of Fumili, for the status
+// returned see TFumili::Minimize.
+
+
+ if( gMinuit->fLimset ) {
+ Printf("ERROR: AliBWTools: Parameters at limits");
+ return kFALSE;
+ }
+
+
+ ///
+ fitter = TVirtualFitter::GetFitter();
+ Int_t fitStat = fitter->GetStats(amin, edm, errdef, nvpar, nparx);
+
+ if( ( (fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ")|| (edm > 1e6) || (fitResult !=0 && fitResult < 4000) ) &&
+ TString(gMinuit->fCstatu) != "SUCCESSFUL" &&
+ TString(gMinuit->fCstatu) != "CONVERGED " ) {
+ if(fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ") {
+ Printf("WARNING: AliBWTools: Cannot properly compute errors");
+ if (fitStat == 0) Printf(" not calculated at all");
+ if (fitStat == 1) Printf(" approximation only, not accurate");
+ if (fitStat == 2) Printf(" full matrix, but forced positive-definite");
+ }
+
+ if (edm > 1e6) {
+
+ Printf("WARNING: AliBWTools: Huge EDM (%f): Fit probably failed!", edm);
+ }
+ if (fitResult != 0) {
+ Printf("WARNING: AliBWTools: Fit Result (%d)", fitResult);
+ }
+
+ Printf ("AliBWTools: Trying Again with Strategy = 2");
+
+ gMinuit->Command("SET STRATEGY 2"); // more effort
+ fitResult = 0;
+ fitResult = h1-> Fit(func,"0","",min,max);
+ fitResult = h1-> Fit(func,"IME0","",min,max);
+ fitResult = h1-> Fit(func,"IME0","",min,max);
+
+ fitter = TVirtualFitter::GetFitter();
+
+ fitStat = fitter->GetStats(amin, edm, errdef, nvpar, nparx);
+
+ if(fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ") {
+ Printf("ERROR: AliBWTools: Cannot properly compute errors");
+ if (fitStat == 0) Printf(" not calculated at all");
+ if (fitStat == 1) Printf(" approximation only, not accurate");
+ if (fitStat == 2) Printf(" full matrix, but forced positive-definite");
+ cout << "[" <<gMinuit->fCstatu<<"]" << endl;
+ return kFALSE;
+ }
+
+ if (edm > 1e6) {
+ Printf("ERROR: AliBWTools: Huge EDM (%f): Fit probably failed!", edm);
+
+ return kFALSE;
+ }
+ if (fitResult != 0) {
+ Printf("ERROR: AliBWTools: Fit Result (%d)", fitResult);
+
+ return kFALSE;
+ }
+
+ gMinuit->Command("SET STRATEGY 1"); // back to normal value
+
+ }
+
+ cout << "---- FIT OK ----" << endl;
+
+ return kTRUE;
+
+}
+
+Int_t AliBWTools::GetLowestNotEmptyBin(TH1*h) {
+
+ // Return the index of the lowest non empty bin in the histo h
+
+ Int_t nbin = h->GetNbinsX();
+ for(Int_t ibin = 1; ibin <= nbin; ibin++){
+ if(h->GetBinContent(ibin)>0) return ibin;
+ }
+
+ return -1;
+
+}
+
+Int_t AliBWTools::GetHighestNotEmptyBin(TH1*h) {
+
+ // Return the index of the highest non empty bin in the histo h
+
+ Int_t nbin = h->GetNbinsX();
+ for(Int_t ibin = nbin; ibin > 0; ibin--){
+ if(h->GetBinContent(ibin)>0) return ibin;
+ }
+
+ return -1;
+
+}
+
+void AliBWTools::GetResiduals(TGraphErrors * gdata, TF1 * func, TH1F ** hres, TGraphErrors ** gres) {
+
+ // Returns a graph of residuals vs point and the res/err distribution
+
+ Int_t npoint = gdata->GetN();
+
+ (*gres) =new TGraphErrors();
+ (*hres) = new TH1F(TString("hres_")+gdata->GetName()+"-"+func->GetName(),
+ TString("hres_")+gdata->GetName()+"-"+func->GetName(),
+ 20,-5,5);
+
+
+ for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+ Float_t x = gdata->GetX()[ipoint];
+ Float_t res = (gdata->GetY()[ipoint] - func->Eval(x))/func->Eval(x);
+ Float_t err = gdata->GetEY()[ipoint]/func->Eval(x);
+ (*hres)->Fill(res/err);
+ (*gres)->SetPoint(ipoint, x, res/err);
+ // (*gres)->SetPointError(ipoint, 0, err);
+
+ }
+
+ (*gres)->SetMarkerStyle(gdata->GetMarkerStyle());
+ (*gres)->SetMarkerColor(gdata->GetMarkerColor());
+ (*gres)->SetLineColor (gdata->GetLineColor());
+ (*gres)->GetHistogram()->GetYaxis()->SetTitle("(data-function)/function");
+ (*hres)->SetMarkerStyle(gdata->GetMarkerStyle());
+ (*hres)->SetMarkerColor(gdata->GetMarkerColor());
+ (*hres)->SetLineColor (gdata->GetLineColor());
+
+
+
+}
+
+void AliBWTools::GetResiduals(TH1F* hdata, TF1 * func, TH1F ** hres, TH1F ** hresVsBin) {
+
+ // Returns an histo of residuals bin by bin and the res/err distribution
+
+ if (!func) {
+ Printf("AliBWTools::GetResiduals: No function provided");
+ return;
+ }
+ if (!hdata) {
+ Printf("AliBWTools::GetResiduals: No data provided");
+ return;
+ }
+
+ (*hresVsBin) = (TH1F*) hdata->Clone(TString("hres_")+hdata->GetName());
+ (*hresVsBin)->Reset();
+ (*hres) = new TH1F(TString("hres_")+hdata->GetName()+"-"+func->GetName(),
+ TString("hres_")+hdata->GetName()+"-"+func->GetName(),
+ 20,-5,5);
+
+ Int_t nbin = hdata->GetNbinsX();
+ for(Int_t ibin = 1; ibin <= nbin; ibin++){
+ if(hdata->GetBinContent(ibin)==0) continue;
+ Float_t res = (hdata->GetBinContent(ibin) - func->Eval(hdata->GetBinCenter(ibin)) ) /
+ func->Eval(hdata->GetBinCenter(ibin));
+ Float_t err = hdata->GetBinError (ibin) / func->Eval(hdata->GetBinCenter(ibin));
+ (*hresVsBin)->SetBinContent(ibin,res);
+ (*hresVsBin)->SetBinError (ibin,err);
+ (*hres)->Fill(res/err);
+
+ }
+
+ (*hresVsBin)->SetMarkerStyle(hdata->GetMarkerStyle());
+ (*hresVsBin)->SetMarkerColor(hdata->GetMarkerColor());
+ (*hresVsBin)->SetLineColor (hdata->GetLineColor() );
+ (*hresVsBin)->GetYaxis()->SetTitle("(data-function)/function");
+ (*hres)->SetMarkerStyle(hdata->GetMarkerStyle());
+ (*hres)->SetMarkerColor(hdata->GetMarkerColor());
+ (*hres)->SetLineColor (hdata->GetLineColor() );
+
+}
+
+void AliBWTools::GetYield(TH1* h, TF1 * f, Double_t &yield, Double_t &yieldError, Float_t min, Float_t max,
+ Double_t *partialYields, Double_t *partialYieldsErrors){
+
+ // Returns the yield extracted from the data in the histo where
+ // there are points and from the fit to extrapolate, in the given
+ // range.
+
+ // Partial yields are also returned if the corresponding pointers are non null
+
+ Int_t bin1 = h->FindBin(min);
+ Int_t bin2 = h->FindBin(max);
+ Float_t bin1Edge = GetLowestNotEmptyBinEdge (h);
+ Float_t bin2Edge = GetHighestNotEmptyBinEdge(h);
+
+ Double_t integralFromHistoError ;
+ Double_t integralFromHisto = h->IntegralAndError(bin1,bin2,integralFromHistoError,"width");
+
+ Double_t integralBelow = min < bin1Edge ? f->Integral(min,bin1Edge) : 0;
+ Double_t integralBelowError = min < bin1Edge ? f->IntegralError(min,bin1Edge) : 0;
+ Double_t integralAbove = max > bin2Edge ? f->Integral(bin2Edge,max) : 0;
+ Double_t integralAboveError = max > bin2Edge ? f->Integral(bin2Edge,max) : 0;
+
+ cout << "GetYield INFO" << endl;
+ cout << " " << bin1Edge << " " << bin2Edge << endl;
+ cout << " " << integralFromHisto << " " << integralBelow << " " << integralAbove << endl;
+ cout << " " << integralFromHistoError << " " << integralBelowError << " " << integralAboveError << endl;
+
+ if(partialYields) {
+ partialYields[0] = integralFromHisto;
+ partialYields[1] = integralBelow;
+ partialYields[2] = integralAbove;
+ }
+ if(partialYieldsErrors) {
+ partialYieldsErrors[0] = integralFromHistoError;
+ partialYieldsErrors[1] = integralBelowError;
+ partialYieldsErrors[2] = integralAboveError;
+ }
+ yield = integralFromHisto+integralBelow+integralAbove;
+ yieldError = TMath::Sqrt(integralFromHistoError*integralFromHistoError+
+ integralBelowError*integralBelowError+
+ integralAboveError*integralAboveError);
+
+}
+
+TGraphErrors * AliBWTools::DivideGraphByFunc(TGraphErrors * g, TF1 * f, Bool_t invert){
+
+ // Divides g/f. If invert == true => f/g
+
+ TGraphErrors * gRatio = new TGraphErrors();
+ Int_t npoint = g->GetN();
+ for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+ Double_t x = g->GetX()[ipoint];
+ Double_t ratio = invert ? f->Eval(x)/g->GetY()[ipoint] :g->GetY()[ipoint]/f->Eval(x);
+ gRatio->SetPoint (ipoint, x, ratio);
+ gRatio->SetPointError(ipoint, 0, g->GetEY()[ipoint]/f->Eval(x));
+ // cout << x << " " << g->GetY()[ipoint] << " " << f->Eval(x) << endl;
+
+ }
+ gRatio->SetMarkerStyle(20);
+ //gRatio->Print();
+ return gRatio;
+
+}
+
+TGraphErrors * AliBWTools::DivideGraphByHisto(TGraphErrors * g, TH1 * h, Bool_t invert){
+
+ // Divides g/h. If invert == true => h/g
+
+
+ TGraphErrors * gRatio = new TGraphErrors();
+ Int_t npoint = g->GetN();
+ for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+ Double_t xj = g->GetX()[ipoint];
+ Double_t yj = g->GetY()[ipoint];
+ Double_t yje = g->GetEY()[ipoint];
+
+ Int_t binData = h->FindBin(xj);
+ Double_t yd = h->GetBinContent(binData);
+ Double_t yde = h->GetBinError(binData);
+ Double_t xd = h->GetBinCenter(binData);
+
+ // cout << TMath::Abs((xd-xj)/xd) << endl;
+
+
+
+ if (yd == 0) continue;
+ // if (binData == 28 ) cout << TMath::Abs(xd-xj)/TMath::Abs(xd) << endl;
+
+ if (TMath::Abs(xd-xj)/TMath::Abs(xd) > 0.01){
+ Printf( "WARNING: bin center (%f) and x graph (%f) are more than 1 %% away, skipping",xd,xj );
+ continue;
+
+ }
+
+ Double_t ratio = invert ? yd/yj : yj/yd;
+
+ gRatio->SetPoint(ipoint, xj, ratio);
+ gRatio->SetPointError(ipoint, 0, TMath::Sqrt(yde*yde/yd/yd + yje*yje/yj/yj)*ratio);
+ // gRatio->SetPointError(ipoint, 0, yje/yj * ratio);
+ }
+
+ return gRatio;
+
+
+}
+
+TH1F * AliBWTools::DivideHistoByFunc(TH1F * h, TF1 * f, Bool_t invert){
+
+ // Divides h/f. If invert == true => f/g
+ // Performs the integral of f on the bin range to perform the ratio
+ // Returns a histo with the same binnig as h
+
+ // Prepare histo for ratio
+ TH1F * hRatio = (TH1F*) h->Clone(TString("hRatio_")+h->GetName()+"_"+f->GetName());
+ hRatio->Reset();
+ // Set y title
+ if(!invert) hRatio->SetYTitle(TString(h->GetName())+"/"+f->GetName());
+ else hRatio->SetYTitle(TString(f->GetName())+"/"+h->GetName());
+
+ // Loop over all bins
+ Int_t nbin = hRatio->GetNbinsX();
+
+ for(Int_t ibin = 1; ibin <= nbin; ibin++){
+ Double_t yhisto = h->GetBinContent(ibin);
+ Double_t yerror = h->GetBinError(ibin);
+ Double_t xmin = h->GetBinLowEdge(ibin);
+ Double_t xmax = h->GetBinLowEdge(ibin+1);
+ Double_t yfunc = f->Integral(xmin,xmax)/(xmax-xmin);
+ Double_t ratio = invert ? yfunc/yhisto : yhisto/yfunc ;
+ Double_t error = yerror/yfunc ;
+ hRatio->SetBinContent(ibin,ratio);
+ hRatio->SetBinError (ibin,error);
+ }
+
+ return hRatio;
+
+}
--- /dev/null
+// ----------------------------------------------------------------------
+// AliBWTools
+//
+// This class provides some tools which can be useful in the analsis
+// of spectra, to fit or transform histograms. See the comments of the
+// individual methods for details
+//
+// Author: M. Floris (CERN)
+// ----------------------------------------------------------------------
+
+#ifndef ALIBWTOOLS_H
+#define ALIBWTOOLS_H
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TObject.h"
+#include "TH1.h"
+
+class TF1;
+class TH1D;
+class TH1F;
+class TGraphErrors;
+#endif
+
+
+
+class AliBWTools : public TObject {
+
+public:
+
+ AliBWTools();
+ ~AliBWTools();
+
+ static TH1 * GetOneOverPtdNdPt(TH1 * hPt) ;
+ static TH1 * GetdNdmtFromdNdpt(TH1 * hpt, Double_t mass);
+ static TH1 * GetdNdPtFromOneOverPt(TH1 * h1Pt) ;
+
+ static TGraphErrors * ConcatenateGraphs(TGraphErrors * g1,TGraphErrors * g2);
+ static TH1F * GetHistoFromGraph(TGraphErrors * g, TH1F* hTemplate) ;
+ static TGraphErrors * GetGraphFromHisto(TH1F * h, Bool_t binWidth = kTRUE) ;
+ static TH1F * CombineHistos(TH1 * h1, TH1 * h2, TH1* htemplate, Float_t renorm1=1.);
+ static void GetFromHistoGraphDifferentX(TH1F * h, TF1 * f, TGraphErrors ** gBarycentre, TGraphErrors ** gXlw);
+ static Float_t GetMean(TH1F * h, Float_t min, Float_t max) ;
+
+ static void GetMean(TF1 * func, Float_t &mean, Float_t &error, Float_t min=0, Float_t max=100);
+ static void GetMeanSquare(TF1 * func, Float_t &mean, Float_t &error, Float_t min=0, Float_t max=100) ;
+
+ static Bool_t Fit (TH1 * h, TF1* f, Float_t min, Float_t max) ;
+
+ static Int_t GetLowestNotEmptyBin(TH1*h);
+ static Int_t GetHighestNotEmptyBin(TH1*h);
+ static Float_t GetLowestNotEmptyBinEdge(TH1*h) { return h->GetBinLowEdge(GetLowestNotEmptyBin(h));}
+ static Float_t GetHighestNotEmptyBinEdge(TH1*h) { return h->GetBinLowEdge(GetHighestNotEmptyBin(h)+1);}
+
+ static void GetResiduals(TGraphErrors * gdata, TF1 * func, TH1F ** hres, TGraphErrors ** gres) ;
+ static void GetResiduals(TH1F* hdata, TF1 * func, TH1F ** hres, TH1F ** hresVsBin) ;
+
+ static void GetYield(TH1* h, TF1 * f, Double_t &yield, Double_t &yieldError, Float_t min = 0,
+ Float_t max = 100, Double_t *partialYields=0, Double_t *partialYieldsErrors=0);
+
+ static TGraphErrors * DivideGraphByFunc (TGraphErrors * g, TF1 * f, Bool_t invert = kFALSE);
+ static TGraphErrors * DivideGraphByHisto(TGraphErrors * g, TH1 * h, Bool_t invert = kFALSE);
+ static TH1F * DivideHistoByFunc (TH1F * h, TF1 * f, Bool_t invert = kFALSE);
+
+private:
+
+ static void GetMoment(TString name, TString var, TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) ;
+ static Double_t GetNormalizedFunc(double * x, double* p);
+
+ static TF1 * fFuncForNormalized; // Function used in GetNormalizedFunc
+
+ AliBWTools(const AliBWTools&); // not implemented
+ AliBWTools& operator=(const AliBWTools&); // not implemented
+
+ ClassDef(AliBWTools,1);
+
+};
+
+#endif
--- /dev/null
+// ------------------------------------------------------------------
+//
+// AliLatexTable
+//
+// This class is used to produce tables using latex sintax. It can
+// than print the table as Latex source code (to be pasted on a paper)
+// or ASCII.
+//
+// The basic idea is that you add columns one after the other with the
+// SetNextCol method and than you insert this row. The SetNextCol
+// method comes in different flavours to alow the insertion of
+// different type of values (numbers with or without errors,
+// strings...).
+//
+// TODO:
+// 1. Make the class drawable
+// 2. Implement vertical lines in ascii print
+// 3. Print output in HTML format
+//
+// Author: Michele Floris, CERN
+// ------------------------------------------------------------------
+
+#include "AliLatexTable.h"
+#include "TString.h"
+#include "TRegexp.h"
+#include "TObjString.h"
+#include "TObjArray.h"
+#include <stdarg.h>
+#include "snprintf.h"
+#include "Varargs.h"
+#include "TPRegexp.h"
+#include "TMath.h"
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliLatexTable);
+
+AliLatexTable::AliLatexTable(Int_t ncol, TString format) : fNcol(0), fFormat(""), fRows(0), fCols(0), fNcolReady(0){
+ // constructor, specify number of cols
+ fNcol = ncol;
+ fFormat = format;
+ fNcolReady = 0;
+ fRows = new TObjArray;
+ fCols = new TObjArray;
+
+ fRows->SetOwner();
+ fCols->SetOwner();
+}
+
+
+AliLatexTable::~AliLatexTable() {
+
+ // dtor
+ if (fRows) delete fRows;
+ if (fCols) delete fCols;
+
+}
+
+void AliLatexTable::SetNextCol(Int_t val){
+ // Set next column in current row - integer
+ char col[200];
+ sprintf(col, " %d ", val);
+ SetNextCol(col);
+
+}
+
+void AliLatexTable::SetNextCol(Int_t val, Int_t err){
+
+ // Set next column in current row - int +- int
+ char col[200];
+ sprintf(col, " $%d \\pm %d$ ", val, err);
+ SetNextCol(col);
+
+}
+
+void AliLatexTable::SetNextCol(Double_t val, Int_t scientificNotation){
+
+ // Set next column in current row - double, specify resolution
+
+ char col[200];
+ if(scientificNotation >= 0) {
+ char format[100];
+
+ // cout << format << endl;
+ Double_t mantissa, exp;
+ GetMantissaAndExpBase10(val, mantissa, exp);
+
+ if (exp == 0) {
+ sprintf(format, " $%%%d.%df $", scientificNotation, scientificNotation);
+ sprintf(col, format, mantissa);
+ }
+ else {
+ sprintf(format, " $%%%d.%df \\cdot 10^{%%0.0f} $", scientificNotation, scientificNotation);
+ sprintf(col, format, mantissa, exp);
+ }
+
+
+
+
+ SetNextCol(col);
+
+ } else if (scientificNotation < -1) {
+ char format[100];
+ sprintf(format, " $%%%d.%df $", -scientificNotation,-scientificNotation);
+ sprintf(col, format , val);
+ SetNextCol(col);
+ } else {
+ sprintf(col, " %f ", val);
+ SetNextCol(col);
+ }
+
+}
+void AliLatexTable::SetNextCol(Double_t val, Double_t err, Int_t scientificNotation){
+
+ // Set next column in current row - double +- double, specify resolution
+
+ // if scientific notation is != -1 is used to determine number of
+ // digits in error
+
+ //if it is > 0 exp notation is used
+
+ //if it is < -1 number is truncated to the number of digits after
+ //point dictated by scientificNotation,
+
+ char col[200];
+ if(scientificNotation >=0 ) {
+
+ Double_t mantissa, exp;
+ GetMantissaAndExpBase10(val, mantissa, exp);
+
+ Double_t mantissa_err, exp_err;
+ GetMantissaAndExpBase10(err, mantissa_err, exp_err);
+
+ Int_t nSigDigits =TMath::Nint(exp - exp_err);
+ if(scientificNotation != 0) nSigDigits = nSigDigits + scientificNotation - 1;
+
+
+ char format[100];
+ if (exp == 0) {
+ sprintf(format, " $%%%d.%df \\pm %%%d.%df $", nSigDigits, nSigDigits, nSigDigits, nSigDigits);
+ sprintf(col, format, mantissa, mantissa_err/TMath::Power(10,exp - exp_err));
+ }
+ else {
+ sprintf(format, " $%%%d.%df \\pm %%%d.%df \\cdot 10^{%%0.0f}$", nSigDigits, nSigDigits, nSigDigits, nSigDigits);
+ sprintf(col, format, mantissa, mantissa_err/TMath::Power(10,exp - exp_err), exp);
+ }
+
+
+ //cout << format << endl;
+
+ SetNextCol(col);
+
+ } else if (scientificNotation < -1) {
+ char format[100];
+ sprintf(format, " $%%%d.%df \\pm %%%d.%df $", -scientificNotation,-scientificNotation,-scientificNotation,-scientificNotation);
+ sprintf(col, format , val, err);
+ SetNextCol(col);
+ }
+ else {
+ sprintf(col, " $%f \\pm %f$ ", val, err);
+ SetNextCol(col);
+ }
+}
+
+void AliLatexTable::SetNextCol(Double_t val, Double_t err, Double_t errSyst, Int_t scientificNotation){
+
+ // Set next column in current row - double +- double +- double, specify resolution
+
+ // if scientific notation is != -1 is used to determine number of
+ // digits in error
+
+ //if it is > 0 exp notation is used
+
+ //if it is < -1 number is truncated to the number of digits after
+ //point dictated by scientificNotation,
+
+ char col[200];
+ if(scientificNotation >=0 ) {
+
+ Double_t mantissa, exp;
+ GetMantissaAndExpBase10(val, mantissa, exp);
+
+ Double_t mantissa_err, exp_err;
+ GetMantissaAndExpBase10(err, mantissa_err, exp_err);
+
+ Double_t mantissa_errsyst, exp_errsyst;
+ GetMantissaAndExpBase10(errSyst, mantissa_errsyst, exp_errsyst);
+
+ Int_t nSigDigits =TMath::Nint(exp - exp_err);
+ if(scientificNotation != 0) nSigDigits = nSigDigits + scientificNotation - 1;
+
+
+ char format[100];
+ if (exp == 0) {
+ sprintf(format, " $%%%d.%df \\pm %%%d.%df \\pm %%%d.%df $",
+ nSigDigits, nSigDigits, nSigDigits, nSigDigits, nSigDigits, nSigDigits);
+ sprintf(col, format, mantissa,
+ mantissa_err/TMath::Power(10,exp - exp_err),
+ mantissa_errsyst/TMath::Power(10,exp - exp_errsyst));
+ }
+ else {
+ sprintf(format, " $%%%d.%df \\pm %%%d.%df \\pm %%%d.%df \\cdot 10^{%%0.0f}$",
+ nSigDigits, nSigDigits, nSigDigits, nSigDigits, nSigDigits, nSigDigits);
+ sprintf(col, format, mantissa,
+ mantissa_err/TMath::Power(10,exp - exp_err),
+ mantissa_errsyst/TMath::Power(10,exp - exp_errsyst),
+ exp);
+ }
+
+
+ //cout << format << endl;
+
+ SetNextCol(col);
+
+ } else if (scientificNotation < -1) {
+ char format[100];
+ sprintf(format, " $%%%d.%df \\pm %%%d.%df \\pm %%%d.%df $", -scientificNotation,-scientificNotation,-scientificNotation,-scientificNotation, -scientificNotation,-scientificNotation);
+ sprintf(col, format , val, err);
+ SetNextCol(col);
+ }
+ else {
+ sprintf(col, " $%f \\pm %f \\pm %f$ ", val, err, errSyst);
+ SetNextCol(col);
+ }
+}
+
+// void AliLatexTable::SetNextColPrintf(const char *va_(fmt), ...){
+
+
+// const Int_t buf_size = 1024; // hope it's long enough
+// char colBuffer[buf_size];
+
+// va_list ap;
+// va_start(ap, va_(fmt));
+// Int_t n = vsnprintf(colBuffer, buf_size, fmt, ap);
+
+// if (n == -1 || n >= buf_size) {
+// Error("SetNextCol", "vsnprintf error!");
+// }
+
+// va_end(ap);
+
+// // for(int iobj=0; iobj<num; iobj++){ //Loop until all objects are added
+// // TH1 * obj = (TH1 *) va_arg(arguments, void *);
+// // returnValue = returnValue && AddObject(obj,(char*)fOptMany.Data(), EColor(iobj+1));
+// // }
+// SetNextCol(colBuffer);
+// }
+
+
+void AliLatexTable::SetNextCol(TString val){
+
+ // Set next column in current row - tstring
+
+ fCols->Add(new TObjString(val));
+ fNcolReady++;
+
+}
+
+void AliLatexTable::InsertCustomRow(TString row){
+ // insert a full row from string. Make sure you have all the columns
+ fRows->Add(new TObjString(row));
+}
+
+
+
+void AliLatexTable::InsertRow(){
+
+ // Insert the row, based on all the individual columns
+
+ if ( fNcolReady != fNcol) {
+ Warning("InsertRow", "Wrong number of cols: %d (!= %d)", fNcolReady, fNcol);
+ }
+
+ TString row = "";
+ for(Int_t icol = 0; icol < fNcol; icol++){
+ row = row + ((TObjString*) fCols->At(icol))->String();
+ if(icol != (fNcol-1)) row = row + " & ";
+ }
+ row+="\\\\";
+
+ fRows->Add(new TObjString(row));
+
+ fNcolReady = 0;
+ fCols->Clear();
+
+}
+
+void AliLatexTable::InsertHline(){
+
+ // insert an horizontal line
+ fRows->Add(new TObjString("\\hline"));
+
+
+}
+
+Int_t * AliLatexTable::GetColWidths() {
+
+ // Compute the width of columns, for nice ascii printing
+
+ Int_t * col_widths= new Int_t[fNcol];
+ Int_t nrow = fRows->GetEntriesFast();
+
+ for(Int_t icol = 0; icol < fNcol; icol++){
+ col_widths[icol] = 0;
+ }
+
+
+ for(Int_t irow = 0; irow < nrow; irow++){
+ TString row(((TObjString*) fRows->At(irow))->String());
+ if(row.Contains("\\hline"))continue;
+
+ StripLatex(row, "ASCII");
+
+ TObjArray * cols = row.Tokenize("&");
+ if(cols->GetEntries() != fNcol) {
+ Error("GetColWidths", "Wrong number of cols in row %s: %d - %d", row.Data(), cols->GetEntries(), fNcol);
+ }
+ for(Int_t icol = 0; icol < fNcol; icol++){
+ Int_t w = ((TObjString *) cols->At(icol))->String().Length();
+ if (w>col_widths[icol]) col_widths[icol] = w;
+ }
+
+ }
+
+ return col_widths;
+}
+
+void AliLatexTable::PrintTable(Option_t * opt){
+
+ // Print the table on screen. You can specify the format with opt.
+ // Currently supported:
+ // "" -> LaTeX
+ // "ASCII" -> plain ASCII
+ // "HTML" -> HTML, to be improved
+
+ if(TString(opt) == "ASCII" || TString(opt)=="HTML") {
+
+ Int_t nrow = fRows->GetEntriesFast();
+
+ Int_t * col_widths = GetColWidths();
+
+ Int_t total_lenght = 0;
+ for(Int_t icol = 0; icol < fNcol; icol++) total_lenght = total_lenght + col_widths[icol] + 2 ;
+
+
+ for(Int_t irow = 0; irow < nrow; irow++){
+ TString row = ((TObjString*) fRows->At(irow))->String();
+ if (row.Contains("\\hline") ){
+ for(Int_t il = 0; il < total_lenght; il++) printf("-");
+ printf("\n");
+ continue;
+ }
+ StripLatex(row, opt);
+ TObjArray * cols = row.Tokenize("&");
+ for(Int_t icol = 0; icol < fNcol; icol++){
+ const char * colstr = ((TObjString *) cols->At(icol))->String().Data();
+ char format [200];
+ sprintf(format, "%%%ds", col_widths[icol] + 2);
+ printf(format, colstr);
+
+ }
+ printf ("\n");
+ }
+
+
+ return;
+ }
+
+
+ cout << "\\begin{tabular}{"<<fFormat<<"}"<<endl;
+
+ Int_t nrow = fRows->GetEntriesFast();
+
+ for(Int_t irow = 0; irow < nrow; irow++){
+ cout << ((TObjString*) fRows->At(irow))->String() << endl;
+ }
+
+ cout << "\\end{tabular}" << endl;
+
+
+}
+
+// void AliLatexTable::ParseExponent(TString &expo){
+// // TString parseExponent = col;
+// TRegexp re = "e[+-][0-9][0-9]";
+// // cout << col << endl;
+
+// if(expo.Contains(re)){
+// Int_t index = expo.Index(re);
+// TString replacement = "\\cdot 10^{";
+// replacement = replacement + expo(index+1, 3) +"}";
+// // cout << replacement <<" --- "<< endl;
+// replacement.ReplaceAll("+","");
+// // cout << "B: " << expo << endl;
+// // cout << "RE " << expo(re) << endl;
+
+// expo.ReplaceAll(expo(re), replacement);
+// // cout << "A: " << expo << endl;
+// // recursion to parse all exponents
+// if (expo.Contains(re)) ParseExponent(expo);
+// }
+// else Warning("", "Error parsing exponent");
+// }
+
+void AliLatexTable::GetMantissaAndExpBase10(Double_t num, Double_t &man, Double_t &exp) {
+
+ // Helper used to get mantissa and exponent
+
+ exp = TMath::Floor(TMath::Log10(TMath::Abs(num)));
+ man = num / TMath::Power(10, exp);
+
+// cout << "" << endl;
+// cout << num << " = " << man << " x10^{"<<exp<<"} " << endl;
+
+
+}
+
+void AliLatexTable::StripLatex(TString &text, TString format) {
+
+ // Strip latex away for ascii and html printing. Replaces latex
+ // command with corresponding text/tags
+
+ text.ReplaceAll("\\cdot", "x");
+ text.ReplaceAll("\\pm", "+-");
+ text.ReplaceAll("$", "");
+ if (format == "ASCII") {
+ text.ReplaceAll("\\right>", ">");
+ text.ReplaceAll("\\left<", "<");
+ text.ReplaceAll("\\rangle", ">");
+ text.ReplaceAll("\\langle", "<");
+ } else if (format == "HTML") {
+ text.ReplaceAll("\\right>", "⟩");
+ text.ReplaceAll("\\left<", "⟨");
+ text.ReplaceAll("\\rangle", "⟩");
+ text.ReplaceAll("\\langle", "⟨");
+ }
+ if(text.Contains("multicolumn")) {
+ // cout << "col " << text.Data() << endl;
+ // in case of multicol span, replace first column with content and
+ // add n empty cols
+ TObjArray * array = TPRegexp("multicolumn\{(.*)\\}\{.*\\}\{(.*)\\}").MatchS(text); // Added double \\ because gcc 4 triggers hard error on unknown escape sequence. Hope it still works...
+ const TString content = ((TObjString *)array->At(2))->GetString();
+ Int_t nspan = ((TObjString *)array->At(1))->GetString().Atoi();
+ text = content;
+// cout << "ns " << nspan << ((TObjString *)array->At(1))->GetString() << endl;
+// cout << "t " << text << endl;
+
+ for(Int_t ispan = 1; ispan < nspan; ispan++){
+ text+=" & ";
+ }
+ // cout << "t " << text << endl;
+
+
+ }
+
+ text.ReplaceAll("\\mathrm", "");
+
+ text.ReplaceAll("\\", "");
+ text.Strip(TString::EStripType(1), ' ');
+
+}
--- /dev/null
+// This class is used to print a table which can be pasted on a latex
+// document
+
+#ifndef ALILATEXTABLE_H
+#define ALILATEXTABLE_H
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TObject.h"
+#include "TString.h"
+class TObjArray;
+
+#endif
+
+
+using namespace std;
+
+class AliLatexTable : public TObject {
+
+public:
+ AliLatexTable() : fNcol(1), fFormat("c"), fRows(0), fCols(0), fNcolReady(0) {;}
+ AliLatexTable(Int_t ncol, TString format);
+ ~AliLatexTable();
+
+ // first you set the value of each column
+ void SetNextCol(Int_t val);
+ void SetNextCol(Int_t val, Int_t err);
+
+ void SetNextCol(Double_t val, Int_t scientificNotation = -1); // if different from -1 gives significant digits
+ void SetNextCol(Double_t val, Double_t err, Int_t scientificNotation = -1);
+ void SetNextCol(Double_t val, Double_t err, Double_t errSyst, Int_t scientificNotation = -1);
+
+ void SetNextCol(TString val);
+// // allows to use printf syntax
+// void SetNextColPrintf(const char *va_(fmt), ...);
+
+ // Then you add the row (it's up to user make sure all columns are
+ // there)
+ void InsertRow();
+
+
+ // insert a row without building it up with the methods. May be
+ // usefull for header or multcol rows
+ void InsertCustomRow(TString row);
+
+
+ void InsertHline();
+
+ void PrintTable(Option_t * opt = "");
+ // void ParseExponent(TString &expo);
+ void GetMantissaAndExpBase10(Double_t num, Double_t &man, Double_t &exp) ;
+
+ // used to print columns with correct width
+ Int_t * GetColWidths();
+ void StripLatex(TString &row, TString format) ;
+
+
+private:
+
+ Int_t fNcol; // number of columns
+ TString fFormat; // latex format (es "c|ccc")
+
+ TObjArray * fRows; // rows
+ TObjArray * fCols; // columns
+
+ Int_t fNcolReady; // number of cols ready to be insert
+
+
+ AliLatexTable(const AliLatexTable&);
+ AliLatexTable& operator=(const AliLatexTable&);
+
+ ClassDef(AliLatexTable, 1)
+
+
+};
+
+#endif
+
--- /dev/null
+#include <iostream>
+#include "TH1F.h"
+#include "TGraphErrors.h"
+#include "AliBWFunc.h"
+#include "AliBWTools.h"
+#include "TF1.h"
+#include "TFile.h"
+#include "TDatabasePDG.h"
+#include "TROOT.h"
+#include "TCanvas.h"
+#include "TFolder.h"
+#include "TStyle.h"
+#include "AliLatexTable.h"
+#include "TLegend.h"
+#include "TVirtualFitter.h"
+#include "TMath.h"
+#include "TH2F.h"
+#include "TSystem.h"
+#include "TLine.h"
+#include "TLatex.h"
+#include "TMath.h"
+
+#include "TASImage.h"
+#include "TPaveText.h"
+#include "ALICEWorkInProgress.C"
+#include "StarPPSpectra.C"
+#include "GetE735Ratios.C"
+
+using namespace std;
+
+
+// A bunch of useful enums and constants
+enum {kPion=0,kKaon,kProton,kNPart};
+enum {kTPC=0,kTOF,kITS,kITSTPC,kK0,kKinks,kCombTOFTPC,kNHist};// "k0" listed here as a kind of PID method...
+const Int_t kNDet = kITS+2;
+enum {kPos=0,kNeg,kNCharge};
+enum {kPhojet=0,kPyTuneAtlasCSC, kPyTuneCMS6D6T, kPyTunePerugia0, kNTunes} ;
+enum {kFitLevi=0, kFitUA1, kFitPowerLaw,
+ kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0,
+ kNFit};
+
+
+// flags, labels and names
+const char * partFlag[] = {"Pion", "Kaon", "Proton"};
+const char * detFlag[] = {"TPC", "TOF", "ITS", "ITS Global", "K0", "Kinks", "Combined TOF + TPC"};
+const char * chargeFlag[] = {"Pos", "Neg"};
+const char * chargeLabel[] = {"Positive", "Negative"};
+const char * partLabel[kNPart][kNCharge] = {{"#pi^{+}", "#pi^{-}"},
+ // {"K^{+} (#times 2)", "K^{-} (#times 2)"},
+ {"K^{+}", "K^{-}"},
+ {"p" , "#bar{p}"}};
+const char * mcTuneName[] = {"Phojet",
+ "Pythia - CSC 306",
+ "Pythia - D6T 109",
+ "Pythia - Perugia0 - 320", };
+const char * funcName[] = { "Levi", "UA1", "PowerLaw", "Phojet", "AtlasCSC", "CMS6D6T", "Perugia0"};
+
+// Style
+//const Int_t marker[] = {25,24,28,20,21}; // set for highlithining marek
+const Int_t marker[] = {20,24,25,28,21}; // standard set
+const Int_t color [] = {1,2,4};
+
+const Int_t mcLineColor[] = {kGreen,kRed,kBlue,kBlack};
+const Int_t mcLineStyle[] = {1,2,3,4};
+
+
+// Template needed to combine different detectors
+const Float_t templBins[] = {0.05,0.1,0.15,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2,2.2,2.4,2.6};
+Int_t nbinsTempl=31;
+
+TH1F * htemplate = new TH1F("htemplate", "htemplate",nbinsTempl, templBins );
+
+// Globals
+TH1F * hSpectra[kNHist][kNPart][kNCharge];
+TH1F * hSpectraMC[kNTunes][kNPart][kNCharge];
+Double_t mass[kNPart];
+
+// Functions:
+// Loading
+void LoadSpectra() ;
+void LoadMC() ;
+
+// Additional tasks (may be uncommented below)
+void DrawStar(Int_t icharge);
+void GetITSResiduals();
+void DrawWithModels() ;
+void DrawAllAndKaons();
+void DrawWithJacek();
+void DrawRatioToStar();
+void DrawRatios();
+
+// External stuff
+//void ALICEWorkInProgress(TCanvas *c,TString today, TString label);
+
+// Used to tag plots
+TDatime dt;
+TString today = "";
+
+
+
+// Switches
+Bool_t convertToMT = 0;
+Bool_t doPrint = 0;
+Int_t fitFuncID = kFitLevi;
+Bool_t scaleKaons = kFALSE;
+Bool_t correctSecondaries = 1;
+Bool_t correctGeantFlukaXS = 1;
+
+void CombineSpectra() {
+
+ // This macro is used to combine the 900 GeV spectra from 2009
+
+ // Load stuff
+ gSystem->Load("libTree.so");
+ gSystem->Load("libVMC.so");
+ gSystem->Load("libMinuit.so");
+ gSystem->Load("libSTEERBase.so");
+ gSystem->Load("libESD.so");
+ gSystem->Load("libAOD.so");
+ gSystem->Load("libANALYSIS.so");
+ gSystem->Load("libANALYSISalice.so");
+ gSystem->Load("libCORRFW.so");
+ gSystem->Load("libPWG2spectra.so");
+
+ // Set date
+ today = today + long(dt.GetDay()) +"/" + long(dt.GetMonth()) +"/"+ long(dt.GetYear());
+
+
+ // Set Masses
+ mass[kPion] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+ mass[kKaon] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
+ mass[kProton] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
+
+ // Load histos
+ LoadSpectra();
+ LoadMC();
+
+ // Additional tasks
+ //DrawStar(icharge);
+ // GetITSResiduals();
+ //DrawAllAndKaons();
+ // DrawWithModels() ;
+ //DrawWithJacek();
+ //DrawRatioToStar();
+ //DrawRatios();
+ //return;
+
+
+ // Draw combined & Fit
+ AliBWFunc * fm = new AliBWFunc;
+ fm->SetVarType(AliBWFunc::kdNdpt);
+ if (convertToMT) fm->SetVarType(AliBWFunc::kOneOverMtdNdmt);
+
+ // table to print results
+ AliLatexTable table(10,"c|cccccccc");
+ if (fitFuncID == kFitLevi) {
+ table.InsertCustomRow("Part & Yield & Yield (FIT) & T Slope & n & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle & \\langle p_{t}^{2} \\rangle");
+ } else if (fitFuncID == kFitPowerLaw) {
+ table.InsertCustomRow("Part & Yield & Norm & n & pt0 & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle & \\langle p_{t}^{2} \\rangle");
+ } else {
+ table.InsertCustomRow("Part & Yield & Par0 & Par2 & Par1 & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle & \\langle p_{t}^{2} \\rangle");
+
+ }
+
+ AliLatexTable tempTable(4,"c|ccc");
+ tempTable.InsertCustomRow("Part & Yield & Yield Below & Frac Above\\\\");
+
+
+ // Fit all
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+
+ TCanvas * c2 = new TCanvas(TString("cCombined")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombined")+chargeFlag[icharge]);
+ TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.001,5);
+ hempty->SetXTitle("p_{t} (GeV/c)");
+ hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
+ hempty->Draw();
+ // DrawStar(icharge);
+ c2->SetLogy();
+ TLegend * l = new TLegend(0.516779, 0.729021, 0.89094 ,0.916084, chargeLabel[icharge]);
+ l->SetFillColor(kWhite);
+
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ Float_t fitmin = 0;
+ Float_t fitmax = 3;
+
+ // Get functions
+ TF1 * func = 0;
+ if(fitFuncID == kFitLevi) {
+ if (ipart == kPion)
+ func = fm->GetLevi(mass[ipart], 0.12, 7, 1.5);
+ if (ipart == kKaon)
+ func = fm->GetLevi(mass[ipart], 0.17, 7, 0.17);
+ if (ipart == kProton)
+ func = fm->GetLevi(mass[ipart], 0.15, 8, 0.08);
+ }
+ else if(fitFuncID == kFitUA1) func = fm->GetUA1(mass[ipart],0.2,1.25,1.25,0.2,1.5);
+ else if(fitFuncID == kFitPowerLaw) {
+ if (ipart == kPion)
+ func = fm->GetPowerLaw(1.0, 8.6, 7);
+ if (ipart == kKaon)
+ func = fm->GetPowerLaw(3.0, 12, 2.6);
+ if (ipart == kProton)
+ func = fm->GetPowerLaw(24, 72, 0.8);
+ }
+ else if(fitFuncID == kFitPhojet) func = fm->GetHistoFunc(hSpectraMC[kPhojet] [ipart][icharge]);
+ else if(fitFuncID == kFitAtlasCSC) func = fm->GetHistoFunc(hSpectraMC[kPyTuneAtlasCSC][ipart][icharge]);
+ else if(fitFuncID == kFitPerugia0) func = fm->GetHistoFunc(hSpectraMC[kPyTunePerugia0][ipart][icharge]);
+ else if(fitFuncID == kFitCMS6D6T) func = fm->GetHistoFunc(hSpectraMC[kPyTuneCMS6D6T] [ipart][icharge]);
+ else {
+ cout << "Unknown fit ID " << fitFuncID << endl;
+ return;
+ }
+
+ if(fitFuncID >= kFitPhojet){
+ fitmin = 0.0;
+ fitmax = 1.0;
+ }
+
+ if(!AliBWTools::Fit(hSpectra[kCombTOFTPC][ipart][icharge],func,fitmin,fitmax)) {
+ cout << " FIT ERROR " << endl;
+ return;
+ }
+ hSpectra[kCombTOFTPC][ipart][icharge]->Draw("same");
+ hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
+ ((TF1*)hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0))->SetRange(0,4);
+ l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge],
+ scaleKaons && ipart == kKaon ?
+ (TString(partLabel[ipart][icharge])+" #times 2").Data()
+ : partLabel[ipart][icharge]);
+// TF1 * fClone = (TF1*) hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Clone();
+// hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->Add(fClone);
+// fClone->SetLineStyle(kDashed);
+// fClone->SetRange(0,100);
+// fClone->Draw("same");
+
+ // populate table
+ // Float_t yield = func->Integral(0.45,1.05);
+ // Float_t yieldE = func->IntegralError(0.45,1.05);
+
+ Float_t yield = func->Integral(0.,100);
+ //Float_t yieldE = func->IntegralError(0.,100);
+
+ Double_t yieldTools, yieldETools;
+ Double_t partialYields[3],partialYieldsErrors[3];
+ AliBWTools::GetYield(hSpectra[kCombTOFTPC][ipart][icharge], func, yieldTools, yieldETools,
+ 0, 100, partialYields,partialYieldsErrors);
+ Double_t tslope = func->GetParameter(2);
+ Double_t tslopeE = func->GetParError(2);
+
+ table.SetNextCol(partLabel[ipart][icharge]);
+ //table.SetNextCol(yield,yieldE,-4);
+ table.SetNextCol(yieldTools, yieldETools,-4);
+ table.SetNextCol(func->GetParameter(0));
+ table.SetNextCol(tslope,tslopeE,-4);
+ table.SetNextCol(func->GetParameter(1));
+ table.SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
+ Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
+ //Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kITS][ipart][icharge]);
+ Float_t yieldAbove = func->Integral(lowestPoint,100);
+ table.SetNextCol(lowestPoint,-2);
+ table.SetNextCol(yieldAbove/yield,-2);
+ Float_t mean, meane;
+ Float_t mean2, mean2e;
+ AliBWTools::GetMean (func, mean, meane );
+ AliBWTools::GetMeanSquare(func, mean2, mean2e);
+ table.SetNextCol(mean, meane ,-4);
+ table.SetNextCol(mean2, mean2e,-4);
+
+ // fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
+ table.InsertRow();
+
+
+ /// TEMP TABLE
+ tempTable.SetNextCol(partLabel[ipart][icharge]);
+ tempTable.SetNextCol(yieldTools, yieldETools, -4);
+ tempTable.SetNextCol(partialYields[1], partialYieldsErrors[1], -4);
+ tempTable.SetNextCol(yieldAbove/yield,-2);
+ tempTable.InsertRow();
+ }
+ l->Draw();
+ if (doPrint) {
+ c2->Update();
+ gSystem->ProcessEvents();
+ c2->Print(TString(c2->GetName()) + ".eps");
+ ALICEWorkInProgress(c2,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
+ c2->Print(TString(c2->GetName()) + ".png");
+ }
+
+ }
+
+
+ table.PrintTable("ASCII");
+
+ cout << "" << endl;
+ tempTable.PrintTable("ASCII");
+
+}
+
+void LoadSpectra() {
+
+ TFile * f=0;
+
+ // Load
+
+
+ // TOF
+ // f = new TFile("./Files/spectra-pos-y.root");
+ f = new TFile("./Files/spectra-pos-y_20100615.root");
+ hSpectra[kTOF][kPion] [kPos]= (TH1F*) f->Get("hpi");
+ hSpectra[kTOF][kProton][kPos]= (TH1F*) f->Get("hpr");
+ hSpectra[kTOF][kKaon] [kPos]= (TH1F*) f->Get("hka");
+ hSpectra[kTOF][kPion] [kPos]->SetName("hpiPos");
+ hSpectra[kTOF][kProton][kPos]->SetName("hprPos");
+ hSpectra[kTOF][kKaon] [kPos]->SetName("hkaPos");
+ //f = new TFile("./Files/spectra-neg-y.root");
+ f = new TFile("./Files/spectra-neg-y_20100615.root");
+ hSpectra[kTOF][kPion] [kNeg]= (TH1F*) f->Get("hpi");
+ hSpectra[kTOF][kProton][kNeg]= (TH1F*) f->Get("hpr");
+ hSpectra[kTOF][kKaon] [kNeg]= (TH1F*) f->Get("hka");
+ hSpectra[kTOF][kPion] [kNeg]->SetName("hpiNeg");
+ hSpectra[kTOF][kProton][kNeg]->SetName("hprNeg");
+ hSpectra[kTOF][kKaon] [kNeg]->SetName("hkaNeg");
+
+
+ // Clean UP TPC spectra, removing unwanted points
+ cout << "Cleaning Up TOF spectra" << endl;
+ Int_t nbin = hSpectra[kTOF][kKaon][kPos]->GetNbinsX();
+ for(Int_t ibin = 1; ibin <= nbin; ibin++){
+ Float_t pt = hSpectra[kTOF][kKaon][kPos]->GetBinCenter(ibin);
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ if (pt > 2.35) {
+ hSpectra[kTOF][kKaon][icharge]->SetBinContent(ibin,0);
+ hSpectra[kTOF][kKaon][icharge]->SetBinError (ibin,0);
+ hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
+ hSpectra[kTOF][kProton][icharge]->SetBinError (ibin,0);
+ }
+ }
+ }
+// cout << "Scaling TOF to TPC" << endl;
+// // Scale protons to TPC
+// hSpectra[kTOF][kProton][kPos]->Scale(1./1.05);
+// // Scale pbar so that pbar/p is compatible with Panos
+// hSpectra[kTOF][kProton][kNeg]->Scale(1./1.1);
+
+
+
+
+
+ // ITS SA (Emanuele)
+ // f = new TFile("./Files/ITSsaSPECTRA_3clusters20100619.root");
+ f = new TFile("./Files/ITSsaSPECTRAperMICHELE_20100703.root");
+ hSpectra[kITS][kPion] [kPos]= (TH1F*) f->Get("hSpectraPos0");
+ hSpectra[kITS][kKaon] [kPos]= (TH1F*) f->Get("hSpectraPos1");
+ hSpectra[kITS][kProton][kPos]= (TH1F*) f->Get("hSpectraPos2");
+ hSpectra[kITS][kPion] [kNeg]= (TH1F*) f->Get("hSpectraNeg0");
+ hSpectra[kITS][kKaon] [kNeg]= (TH1F*) f->Get("hSpectraNeg1");
+ hSpectra[kITS][kProton][kNeg]= (TH1F*) f->Get("hSpectraNeg2");
+
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ Int_t nbin = hSpectra[kITS][ipart][icharge]->GetNbinsX();
+ // hSpectra[kITS][ipart][icharge]->Scale(276004.);// Emanule divided his spectra...
+ for(Int_t ibin = 1; ibin <= nbin; ibin++){
+ if(hSpectra[kITS][ipart][icharge]->GetBinContent(ibin) < 0 ) {
+ hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
+ hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
+ }
+// if ((ipart == kKaon && ibin >= 12) || (ipart == kProton && ibin >= 20)) {
+// hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
+// hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
+// }
+ }
+
+ }
+ }
+
+
+ // ITS + TPC (Marek)
+ f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons_20100618.root");
+ TList * list = (TList*) gDirectory->Get("output");
+ hSpectra[kITSTPC][kPion] [kPos]= (TH1F*) list->FindObject("Pions");
+ hSpectra[kITSTPC][kKaon] [kPos]= (TH1F*) list->FindObject("Kaons");
+ hSpectra[kITSTPC][kProton][kPos]= (TH1F*) list->FindObject("Protons");
+ hSpectra[kITSTPC][kPion] [kNeg]= (TH1F*) list->FindObject("AntiPions");
+ hSpectra[kITSTPC][kKaon] [kNeg]= (TH1F*) list->FindObject("AntiKaons");
+ hSpectra[kITSTPC][kProton][kNeg]= (TH1F*) list->FindObject("AntiProtons");
+
+ // TPC
+ // htemplate = hSpectra[kITS][kProton][kNeg]; //FIXME
+ f = new TFile("./Files/protonSpectra_20100615.root");
+ hSpectra[kTPC][kProton][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonPosClassic"),htemplate);
+ hSpectra[kTPC][kProton][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonNegClassic"),htemplate);
+ f = new TFile("./Files/pionSpectra_20100615.root");
+ hSpectra[kTPC][kPion][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionPosClassic"),htemplate);
+ hSpectra[kTPC][kPion][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionNegClassic"),htemplate);
+ f = new TFile("./Files/kaonSpectra_20100615.root");
+ hSpectra[kTPC][kKaon][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonPosClassic"),htemplate);
+ hSpectra[kTPC][kKaon][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonNegClassic"),htemplate);
+
+ // Clean UP TPC spectra, removing unwanted points
+ cout << "Cleaning Up TPC spectra" << endl;
+ nbin = hSpectra[kTPC][kKaon][kPos]->GetNbinsX();
+ for(Int_t ibin = 0; ibin < nbin; ibin++){
+ Float_t pt = hSpectra[kTPC][kKaon][kPos]->GetBinCenter(ibin);
+ if (pt > 0.45) {
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ hSpectra[kTPC][kKaon][icharge]->SetBinContent(ibin,0);
+ hSpectra[kTPC][kKaon][icharge]->SetBinError (ibin,0);
+ }
+ }
+ if (pt < 0.45) {
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ hSpectra[kTPC][kProton][icharge]->SetBinContent(ibin,0);
+ hSpectra[kTPC][kProton][icharge]->SetBinError (ibin,0);
+ }
+ }
+ }
+
+
+
+
+
+ // K0s
+ f = new TFile ("./Files/PtSpectraCorrectedK0sOff-dNdy-Jun02.root");
+ // hSpectra[kK0][kKaon][kPos] = (TH1F*) AliBWTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get("hSpectraOff"));
+ hSpectra[kK0][kKaon][kPos] = (TH1F*) gDirectory->Get("hSpectraOff");
+ // hSpectra[kK0][kKaon][kPos]->Scale(2*TMath::Pi());
+ // hSpectra[kK0][kKaon][kPos]->Scale(1./272463);
+ hSpectra[kK0][kKaon][kNeg] = hSpectra[kK0][kKaon][kPos];
+
+ // Kinks: TO BE FIXED WITH POSITIVES AND NEGATIVES
+ // f = new TFile ("./Files/PtAllKaonKinkRap6Apr24.root");
+ f = new TFile ("./Files/PtKaonKinkJune13AllPN_20100615.root");
+ hSpectra[kKinks][kKaon][kPos] = (TH1F*)gDirectory->Get("fptallKPA");
+ hSpectra[kKinks][kKaon][kNeg] = (TH1F*)gDirectory->Get("fptallKNA");
+
+
+ // Apply correction factrs
+ // Secondaries for protons
+ f = new TFile ("./Files/corrFactorProtons_20100615.root");
+ TH1F * hCorrSecondaries = (TH1F*) gDirectory->Get("corrFactorProtons");
+ if(correctSecondaries) {
+ cout << "CORRECTING SECONDARIES" << endl;
+
+ for(Int_t idet = 0; idet <= kTOF; idet++){ // TPC and TOF only
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ Int_t ipart = kProton;
+ TH1* h = hSpectra[idet][ipart][icharge]; // lighter notation below
+ if (h){
+ Int_t nbins = h->GetNbinsX();
+ for(Int_t ibin = 0; ibin < nbins; ibin++){
+ // Float_t pt = convertToMT ? TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) : h->GetBinCenter(ibin);
+ Float_t pt = h->GetBinCenter(ibin);
+ if (icharge == kNeg) pt = -pt;
+ Int_t binCorrection = hCorrSecondaries->FindBin(pt);
+ Float_t correction = hCorrSecondaries->GetBinContent(binCorrection);
+ // cout << pt << " " << correction << endl;
+
+ if (correction != 0) {// If the bin is empty this is a 0
+ h->SetBinContent(ibin,h->GetBinContent(ibin)/correction);
+ h->SetBinError (ibin,h->GetBinError (ibin)/correction);
+ } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+ cout << "Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
+ cout << " Bin content: " << h->GetBinContent(ibin) << endl;
+
+ }
+ }
+ }
+ }
+ }
+ }
+ // geant/fluka absorption
+ if(correctGeantFlukaXS) {
+ cout << "CORRECTING GEANT3/FLUKA" << endl;
+
+ f = new TFile ("./Files/correctionForCrossSection_20100615.root");
+ TH2D * hCorrFluka[kNCharge];
+ hCorrFluka[kPos] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionProtons");
+ hCorrFluka[kNeg] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionAntiProtons");
+ for(Int_t idet = 0; idet < kNDet; idet++){
+ if (idet == kITS) continue; // skip ITS sa
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ Int_t ipart = kProton;
+ TH1 * h = hSpectra[idet][ipart][icharge]; // lighter notation below
+ if (h){
+ Int_t nbins = h->GetNbinsX();
+ for(Int_t ibin = 0; ibin < nbins; ibin++){
+// Float_t pt = convertToMT ?
+// TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) :
+// h->GetBinCenter(ibin);
+ Float_t pt = h->GetBinCenter(ibin);
+ Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
+ Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(h->GetNbinsY()+1);
+ if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
+ if (pt > maxPtCorrection) pt = maxPtCorrection;
+ Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
+ if (idet == kTOF && icharge == kNeg) {
+
+ // Apply parametrized correction computed by francesco
+ // Fitted panos correction and using momentum at the outer radius of the TPC
+
+ Float_t ptav = pt; // Just to use the same name francesco uses...
+
+ // from pT constrained at P.V. (ptav) to pT TPC outer (ptTPCout)
+ Float_t ptTPCout=ptav*(1-6.81059e-01*TMath::Exp(-ptav*4.20094));
+
+ // traking correction (fit to Panos)
+ Float_t antiprotonEC = 1 - 0.129758 *TMath::Exp(-ptav*0.679612);
+
+ // TOF matching efficiency correction (derived from Panos one scaled for M.B.(TOF)/M.B.(TPC)).
+ Float_t antiprotonEC2 = TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCout*0.679612),0.07162/0.03471);
+ correction = antiprotonEC * antiprotonEC2;
+ }
+ // cout << icharge<< " " << h->GetBinCenter(ibin) << " " << pt << " " << correction << endl;
+ if (correction != 0) {// If the bin is empty this is a 0
+ h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
+ h->SetBinError (ibin,h->GetBinError (ibin)*correction);
+// if (idet == kTOF) {
+// cout << "CORRECTING TOF TWICE" << endl;
+// h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
+// h->SetBinError (ibin,h->GetBinError (ibin)*correction);
+// }
+ } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+ cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
+ cout << " Bin content: " << h->GetBinContent(ibin) << endl;
+ }
+
+ }
+
+ }
+ }
+ }
+ }
+
+
+ // Set style and scale
+ for(Int_t idet = 0; idet < kNDet; idet++){
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ if (hSpectra[idet][ipart][icharge]){
+ hSpectra[idet][ipart][icharge]->SetStats(0); // disable stats
+ hSpectra[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
+ hSpectra[idet][ipart][icharge]->SetLineColor (color[ipart] );
+ hSpectra[idet][ipart][icharge]->SetMarkerStyle (marker[ipart]);
+ hSpectra[idet][ipart][icharge]->SetXTitle("p_{t} (GeV/c)");
+ hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} dN/dp_{t} (GeV/c)^{-1}");
+ if (convertToMT) {
+ TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[idet][ipart][icharge]);
+ hSpectra[idet][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
+ hSpectra[idet][ipart][icharge]->SetXTitle("m_{t} (GeV)");
+ hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
+ }
+// if (idet == kTOF || idet == kTPC) {
+// hSpectra[idet][ipart][icharge]->Scale(1./236763);
+// }
+// if (idet == kITS){
+// hSpectra[idet][ipart][icharge]->Scale(202945./236763);
+// }
+ if (scaleKaons && ipart == kKaon) {
+ hSpectra[idet][ipart][icharge]->Scale(2.);
+ }
+ } else {
+ printf("Cannot load %s,%s,%s\n",detFlag[idet],partFlag[ipart],chargeFlag[icharge]);
+ }
+ }
+ }
+ }
+
+
+
+
+ // Combine detectors
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ TH1F * htemplLocal = htemplate; // If we are converting to 1/mt dNdmt we need to convert the template as well...
+
+ if (convertToMT) {
+ TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(htemplate);
+ htemplLocal = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
+
+ }
+ hSpectra[kCombTOFTPC][ipart][icharge] = AliBWTools::CombineHistos(hSpectra[kTOF][ipart][icharge],
+ hSpectra[kTPC][ipart][icharge],
+ htemplLocal,1.);;
+// if (convertToMT) {
+// TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[kCombTOFTPC][ipart][icharge]);
+// hSpectra[kCombTOFTPC][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
+// hSpectra[kCombTOFTPC][ipart][icharge]->SetXTitle("m_{t} (GeV)");
+// hSpectra[kCombTOFTPC][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
+// }
+ }
+ }
+
+
+ // Scale for the number of inelastic collisions and correct for
+ // efficiency losses due to physics selection:
+
+ Double_t effPhysSel[kNPart];
+ effPhysSel[kPion] = 1.012;
+ effPhysSel[kKaon] = 1.013;
+ effPhysSel[kProton] = 1.014;
+
+
+ for(Int_t idet = 0; idet < kNHist; idet++){
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ if(hSpectra[idet][ipart][icharge]) {
+ // cout << "Scaling!" << endl;
+ hSpectra[idet][ipart][icharge]->Scale(1.*effPhysSel[ipart]/278366.15); // Scale PhysSel tutti? // FIXME
+ }
+ }
+ }
+ }
+
+
+}
+
+void LoadMC() {
+
+ TFile * f = 0;
+ const char * evClass= "INEL";
+ const char * files[] = {"./Files/dndeta_Phojet_10M_900GeV.root",
+ "./Files/dndeta_AtlasCSC306_10M_900GeV.root",
+ "./Files/dndeta_CMSD6T109_10M_900GeV.root",
+ "./Files/dndeta_Perugia0320_10M_900GeV.root", };
+
+ // Phojet
+ for(Int_t itune = 0; itune < kNTunes; itune++){
+ f = new TFile(files[itune]);
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ hSpectraMC[itune][ipart][icharge] = (TH1F*) f->Get(Form("fHistPtID_%s_%s%s",evClass,partFlag[ipart],icharge==0 ? "Pos" : "Neg"));
+ }
+ }
+ }
+
+
+ // Set style
+ for(Int_t itune = 0; itune < kNTunes; itune++){
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ if (hSpectraMC[itune][ipart][icharge]){
+ hSpectraMC[itune][ipart][icharge]->SetName(TString(hSpectraMC[itune][ipart][icharge]->GetName())+mcTuneName[itune]);
+ hSpectraMC[itune][ipart][icharge]->SetMarkerColor (mcLineColor[itune] );
+ hSpectraMC[itune][ipart][icharge]->SetLineColor (mcLineColor[itune] );
+ hSpectraMC[itune][ipart][icharge]->SetLineStyle (mcLineStyle[itune] );
+ hSpectraMC[itune][ipart][icharge]->SetMarkerStyle (1);
+ if (convertToMT) {
+ TH1F * htmp = (TH1F*)AliBWTools::GetOneOverPtdNdPt(hSpectraMC[itune][ipart][icharge]);
+ hSpectraMC[itune][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
+ hSpectraMC[itune][ipart][icharge]->SetXTitle("m_{t} (GeV)");
+ hSpectraMC[itune][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
+ }
+
+ } else {
+ printf("Cannot load MC # %d,%s,%s\n",itune,partFlag[ipart],chargeFlag[icharge]);
+ }
+ }
+ }
+ }
+
+}
+
+void DrawStar(Int_t icharge) {
+ // cout << icharge << endl;
+
+ // gROOT->LoadMacro("StarPPSpectra.C");
+ TGraphErrors ** gStar = StarPPSpectra();
+
+ for(Int_t istar = 0; istar < 6; istar++){
+ gStar[istar]->SetMarkerStyle(kOpenStar);
+ if (icharge==kPos && (istar%2) ) gStar[istar]->Draw("P");
+ else if (icharge==kNeg && !(istar%2) ) gStar[istar]->Draw("P");
+ else cout << "Skipping Star " << istar << endl;
+ }
+}
+
+void GetITSResiduals() {
+
+
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ // cout << "1 " << ipart << " " << icharge << endl;
+
+// gSystem->ProcessEvents();
+// gSystem->Sleep(1000);
+ TF1* f = (TF1*) hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0);
+ TH1F * hres1, *hres2;
+ AliBWTools::GetResiduals(hSpectra[kITS][ipart][icharge], f, &hres1, &hres2);
+
+ TCanvas * c1 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge],TString(partFlag[ipart])+"_"+chargeFlag[icharge]);
+ c1->SetLogy();
+ hSpectra[kCombTOFTPC][ipart][icharge]->Draw();
+ hSpectra[kITS][ipart][icharge]->SetMarkerStyle(24);
+ hSpectra[kCombTOFTPC][ipart][icharge]->SetMarkerStyle(20);
+ hSpectra[kITS][ipart][icharge]->Draw("same");
+ hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
+ TLegend* l = new TLegend( 0.182886, 0.192308, 0.505034,0.384615, TString(partLabel[ipart][icharge])+" "+chargeFlag[icharge]);
+ l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge], "TOF+TPC");
+ l->AddEntry(hSpectra[kITS][ipart][icharge], "ITS");
+ l->AddEntry(f, "Fit to TOF+TPC");
+ l->Draw();
+
+ TCanvas * c2 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res",
+ TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res");
+ c2->SetGridy();
+ hres2->SetMinimum(-1);
+ hres2->SetMaximum(1);
+ hres2->Draw();
+ hres2->GetYaxis()->SetTitleOffset(1.2);
+ Float_t x = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
+ TLine * line = new TLine(x,-1,x,1);
+ line->SetLineStyle(kDashed);
+ line->Draw("same");
+
+ if (doPrint) {
+ c1->Update();
+ c2->Update();
+ gSystem->ProcessEvents();
+ c1->Print(TString(c1->GetName()) + ".png");
+ c2->Print(TString(c2->GetName()) + ".png");
+ }
+ }
+ }
+}
+
+void DrawWithModels() {
+
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+
+ // Draw with models
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ // Multipad canvas
+ TCanvas * c1 = new TCanvas(TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge],TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge] );
+ TPad *p1 = new TPad(TString("p1")+partFlag[ipart]+chargeFlag[icharge], "p1", 0.0, 0.3, 1.0, 0.95, 0, 0, 0);
+ p1->SetBottomMargin(0);
+ p1->Draw();
+
+ TPad *p2 = new TPad(TString("p2")+partFlag[ipart]+chargeFlag[icharge], "p2", 0.0, 0.05, 1.0, 0.3, 0, 0, 0);
+ p2->SetTopMargin(0);
+ p2->SetBottomMargin(0.4);
+ p2->Draw();
+
+ Float_t scaleFonts = (0.95-0.3)/(0.3-0.05);
+
+ // Draw spectra
+ p1->cd();
+ p1->SetLogy();
+ TH2F * hempty = new TH2F(TString("hempty")+long(ipart),"hempty",100,0.,4, 100, 0.0015,5);
+ hempty->SetXTitle("p_{t} (GeV/c)");
+ hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
+ hempty->Draw();
+ c1->SetLogy();
+
+
+ TLegend * l =new TLegend( 0.543624, 0.431818, 0.892617,0.629371);
+ l->SetFillColor(kWhite);
+ hSpectra[kCombTOFTPC][ipart][icharge]->Draw("same");
+ l->AddEntry(hSpectra[kTOF][ipart][icharge],TString ("Data: ")+partLabel[ipart][icharge]);
+ for(Int_t itune = 0; itune < kNTunes; itune++){
+ l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
+ hSpectraMC[itune][ipart][icharge]->SetLineWidth(2);
+ hSpectraMC[itune][ipart][icharge]->Draw("same,chist");
+ }
+ l->Draw("same");
+
+ // Draw ratio
+ p2->cd();
+ TH2F * hemptyr = new TH2F(TString("hemptyratio")+long(ipart),"hempty",100,0.,4, 100, 0.01,2.99);
+ hemptyr->SetXTitle("p_{t} (GeV/c)");
+ hemptyr->SetYTitle("Data/MC");
+ hemptyr->GetXaxis()->SetLabelSize(0.04*scaleFonts);
+ hemptyr->GetYaxis()->SetLabelSize(0.04*scaleFonts);
+ hemptyr->GetYaxis()->SetTitleSize(0.05*scaleFonts);
+ hemptyr->GetYaxis()->SetTitleOffset(1.4/scaleFonts);
+ hemptyr->GetXaxis()->SetTitleSize(0.05*scaleFonts);
+ hemptyr->SetTickLength(0.03*scaleFonts, "X");
+ hemptyr->SetTickLength(0.02*scaleFonts, "Y");
+ // hemptyr->GetXaxis()->SetTitleOffset(1.4/scaleFonts);
+ hemptyr->GetYaxis()->SetNdivisions(5);
+ hemptyr->Draw("");
+
+ AliBWFunc fm;
+ for(Int_t itune = 0; itune < kNTunes; itune++){
+ TF1* f = fm.GetHistoFunc(hSpectraMC[itune][ipart][icharge], TString("f")+mcTuneName[itune]);
+
+ // l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
+ TH1F* hRatio = AliBWTools::DivideHistoByFunc(hSpectra[kCombTOFTPC][ipart][icharge],f);
+ hRatio->SetLineStyle(hSpectraMC[itune][ipart][icharge]->GetLineStyle());
+ hRatio->SetLineColor(hSpectraMC[itune][ipart][icharge]->GetLineColor());
+ hRatio->SetLineWidth(hSpectraMC[itune][ipart][icharge]->GetLineWidth());
+ hRatio->Draw("lhist,same");
+ }
+
+
+ // print
+ if(doPrint) {
+ c1->Update();
+ gSystem->ProcessEvents();
+ c1->Print(TString(c1->GetName())+".eps");
+ ALICEWorkInProgress(c1,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
+ c1->Print(TString(c1->GetName())+".png");
+
+ }
+ }
+ }
+
+
+
+}
+
+void DrawAllAndKaons() {
+
+
+ // gROOT->LoadMacro("ALICEWorkInProgress.C");
+
+ gStyle->SetOptFit(0);
+
+ TH1F * hKaonsAllTPCTOF = (TH1F*) hSpectra[kCombTOFTPC][kKaon][kPos]->Clone();
+ hKaonsAllTPCTOF->Add(hSpectra[kCombTOFTPC][kKaon][kNeg]);
+
+ TH1F * hK0Scaled = (TH1F*) hSpectra[kK0][kKaon][kPos]->Clone();
+ hK0Scaled->Add(hSpectra[kK0][kKaon][kPos]);
+
+ hSpectra[kKinks][kKaon][kPos]->SetMarkerStyle(25);
+ hSpectra[kKinks][kKaon][kPos]->SetStats(0);
+ TH1F * hKinksAll = (TH1F*) hSpectra[kKinks][kKaon][kPos]->Clone();
+ hKinksAll->Add(hSpectra[kKinks][kKaon][kNeg]);
+
+ TCanvas * c1 = new TCanvas("cKaons","cKaons");
+ c1->SetLogy();
+ TH2F * hempty = new TH2F("hempty_allkaons","hempty",100,0.,3, 100, 1e-4,10);
+ hempty->SetXTitle("p_{t} (GeV/c)");
+ hempty->SetYTitle("dN / dp_{t} (A.U.)");
+ hempty->Draw();
+ hKaonsAllTPCTOF->Draw("same");
+ hK0Scaled->Draw("same");
+ hKinksAll->Draw("same");
+
+ TLegend * leg = new TLegend(0.2013423,0.2255245,0.5503356,0.4335664,NULL,"brNDC");
+// leg->SetBorderSize(0);
+// leg->SetLineColor(1);
+// leg->SetLineStyle(1);
+// leg->SetLineWidth(1);
+// leg->SetFillColor(19);
+// leg->SetFillStyle(1001);
+ TLegendEntry *entry=leg->AddEntry("hkaPos_h_tpcKaonPos","K^{+} + K^{-}, TPC+TOF ","lpf");
+ entry=leg->AddEntry("h1PtSpectraOff_inv","K^{0} #times 2","lpf");
+ entry=leg->AddEntry("fptallK","K^{+} + K ^{-}, Kinks","lpf");
+ leg->Draw();
+
+ ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Performance}{Not fully corrected}");
+ TLatex * tex = new TLatex(0.2120805,0.01288336,"Statistical error only");
+ tex->SetTextColor(2);
+ tex->SetTextFont(42);
+ tex->SetTextSize(0.03496503);
+ tex->Draw();
+
+ c1->Update();
+ if(doPrint) c1->Print(TString(c1->GetName())+".png");
+
+ // Draw all "stable" hadrons
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ TCanvas * c1 = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge]);
+ c1->SetLogy();
+ TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,4, 100, 1e-4,10);
+ hempty->SetXTitle("p_{t} (GeV/c)");
+ hempty->SetYTitle("dN / dp_{t} (A.U.)");
+ hempty->Draw();
+ leg = new TLegend( 0.645973, 0.325175, 0.892617,0.636364, NULL,"brNDC");
+ for(Int_t ipart = 0; ipart < kNPart; ipart++) {
+ for(Int_t idet = 0; idet <= kITSTPC; idet++){
+ // if (idet == kITS) continue;
+ // if (idet == kITSTPC) hSpectra[idet][ipart][icharge]->SetMarkerColor(kGreen);
+ hSpectra[idet][ipart][icharge]->SetMarkerStyle(marker[idet]);
+ hSpectra[idet][ipart][icharge]->Draw("same");
+ leg->AddEntry(hSpectra[idet][ipart][icharge],TString(partLabel[ipart][icharge])+" (" + detFlag[idet] + ")","lpf");
+ }
+ // leg->AddLine();
+ }
+ leg->Draw();
+ ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Performance}{Not Fully Corrected}");
+ c1->Update();
+ if(doPrint) c1->Print(TString(c1->GetName())+".png");
+ }
+
+
+ // Draw ratios (tmp)
+
+// new TCanvas;
+// hSpectra[kTPC][kKaon][kNeg]->Divide(hSpectra[kTPC][kKaon][kPos]);
+// hSpectra[kTPC][kKaon][kNeg]->Draw();
+// new TCanvas;
+// hSpectra[kITS][kKaon][kNeg]->Divide(hSpectra[kITS][kKaon][kPos]);
+// hSpectra[kITS][kKaon][kNeg]->Draw("");
+ // hSpectra[kTOF][kProton][kPos]->Draw("same");
+
+// for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+// TCanvas * c1 = new TCanvas(TString("cAllRatio_")+chargeFlag[icharge],TString("cAllRatio_")+chargeFlag[icharge]);
+// c1->SetGridy();
+// TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.2,1, 100, 0.0,2.0);
+// hempty->SetXTitle("p_{t} (GeV/c)");
+// hempty->SetYTitle("ITSsa / TPC");
+// hempty->Draw();
+
+
+// for(Int_t ipart = 0; ipart < kNPart; ipart++) {
+// hSpectra[kITS][ipart][icharge]->Divide(hSpectra[kTPC][ipart][icharge]);
+// hSpectra[kITS][ipart][icharge]->Draw("same");
+// }
+// if(doPrint) c1->Print(TString(c1->GetName())+".png");
+// }
+
+// end of tmp
+
+}
+
+void DrawWithJacek() {
+
+ //1. Convert spectra to dNdeta and sum
+ TH1F * hsum = (TH1F*) htemplate->Clone();
+ hsum->Reset();
+ Int_t idet= kCombTOFTPC;
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ TH1 * h = hSpectra[idet][ipart][icharge];
+ Int_t nbin = h->GetNbinsX();
+ for(Int_t ibin = 1; ibin <= nbin; ibin++){
+ Double_t pt = h->GetBinCenter(ibin);
+ Double_t mt = TMath::Sqrt(pt*pt + mass[ipart]*mass[ipart]);
+ Double_t jacobian = pt/mt;
+ h->SetBinContent(ibin,h->GetBinContent(ibin)*jacobian);
+ h->SetBinError (ibin,h->GetBinError (ibin)*jacobian);
+ Int_t ibinSum = hsum->FindBin(pt);
+ Double_t epsilon = 0.0001;
+ if ( h->GetBinContent(ibin) > 0 &&
+ (TMath::Abs(h->GetBinLowEdge(ibin) - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
+ TMath::Abs(h->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
+ ) {
+ cout << "DISCREPANCY IN BIN RANGES" << endl;
+ cout << pt << " " << ibinSum << " " << ibin << "; " << h->GetBinContent(ibin) << endl
+ << h->GetBinLowEdge(ibin) << "-" << h->GetBinLowEdge(ibin+1) << endl
+ << hsum->GetBinLowEdge(ibinSum) << "-" << hsum->GetBinLowEdge(ibinSum+1) << endl;
+ cout << "" << endl;
+ }
+
+
+ hsum->SetBinContent(ibinSum,hsum->GetBinContent(ibinSum)+h->GetBinContent(ibin)); // EROOR FIXME
+ hsum->SetBinError (ibinSum,0);
+ }
+ // hsum->Add(h);
+ }
+ }
+
+
+ // Load Jacek and Draw both:
+ new TFile ("./Files/dNdPt_Data_Points_ALICE_900GeV.root");
+ TGraphErrors * gJacek = (TGraphErrors*) gDirectory->Get("inel");
+ gJacek->Draw("AP");
+ hsum->Draw("same");
+
+ TGraphErrors * gRatio = AliBWTools::DivideGraphByHisto(gJacek,hsum);
+
+
+ new TCanvas();
+ gRatio->Draw("AP");
+
+
+
+}
+
+
+void DrawRatioToStar() {
+
+ // Star data
+ // gROOT->LoadMacro("StarPPSpectra.C");
+ TGraphErrors ** gStar = StarPPSpectra();
+
+ // ALICE, INEL -> NSD
+ Double_t scaleYield = 3.58/3.02; // from paper 2
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ hSpectra[kCombTOFTPC][ipart][icharge]->Scale(scaleYield);
+ }
+ }
+
+
+ TCanvas * c1 = new TCanvas("cRatioToStarNeg","cRatioToStarNeg");
+ TH2F * hempty = new TH2F(TString("hemptyNeg"),"hemptyNeg",100,0.,1.5, 100, 0.001,1.8);
+ hempty->SetXTitle("p_{t} (GeV/c)");
+ hempty->SetYTitle("ALICE/STAR (NSD)");
+
+ TLegend *leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Negative","brNDC");
+ leg->SetBorderSize(0);
+ leg->SetLineColor(1);
+ leg->SetLineStyle(1);
+ leg->SetLineWidth(1);
+ leg->SetFillColor(0);
+ leg->SetFillStyle(1001);
+
+
+ hempty->Draw();
+ TGraphErrors * g ;
+ g = AliBWTools::DivideGraphByHisto(gStar[0],hSpectra[kCombTOFTPC][kPion][kNeg],1);
+ g->SetMarkerStyle(kFullCircle);
+ g->SetMarkerColor(kBlack);
+ g->Draw("p");
+ leg->AddEntry(g,"#pi^{-}","lp");
+ g = AliBWTools::DivideGraphByHisto(gStar[2],hSpectra[kCombTOFTPC][kKaon][kNeg],1);
+ g->SetMarkerStyle(kOpenCircle);
+ g->SetMarkerColor(kRed);
+ g->Draw("p");
+ leg->AddEntry(g,"K^{-}","lp");
+ g = AliBWTools::DivideGraphByHisto(gStar[4],hSpectra[kCombTOFTPC][kProton][kNeg],1);
+ g->SetMarkerStyle(kOpenSquare);
+ g->SetMarkerColor(kBlue);
+ g->Draw("p");
+ leg->AddEntry(g,"#bar{p}","lp");
+ leg->Draw();
+
+
+
+
+ TCanvas * c2 = new TCanvas("cRatioToStarPos","cRatioToStarPos");
+ hempty->Draw();
+ leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Positive","brNDC");
+ leg->SetBorderSize(0);
+ leg->SetLineColor(1);
+ leg->SetLineStyle(1);
+ leg->SetLineWidth(1);
+ leg->SetFillColor(0);
+ leg->SetFillStyle(1001);
+ // TGraphErrors * g ;
+ g = AliBWTools::DivideGraphByHisto(gStar[1],hSpectra[kCombTOFTPC][kPion][kPos],1);
+ g->SetMarkerStyle(kFullCircle);
+ g->SetMarkerColor(kBlack);
+ g->Draw("p");
+ leg->AddEntry(g,"#pi^{+}","lp");
+ g = AliBWTools::DivideGraphByHisto(gStar[3],hSpectra[kCombTOFTPC][kKaon][kPos],1);
+ g->SetMarkerStyle(kOpenCircle);
+ g->SetMarkerColor(kRed);
+ g->Draw("p");
+ leg->AddEntry(g,"K^{+}","lp");
+ g = AliBWTools::DivideGraphByHisto(gStar[5],hSpectra[kCombTOFTPC][kProton][kPos],1);
+ g->SetMarkerStyle(kOpenSquare);
+ g->SetMarkerColor(kBlue);
+ g->Draw("p");
+ leg->AddEntry(g,"p","lp");
+
+
+ leg->Draw();
+
+ c1->Update();
+ c2->Update();
+ gSystem->ProcessEvents();
+ c1->Print(TString(c1->GetName()) + ".eps");
+ c2->Print(TString(c2->GetName()) + ".eps");
+ ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
+ ALICEWorkInProgress(c2,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
+ c1->Update();
+ c2->Update();
+ c1->Print(TString(c1->GetName()) + ".png");
+ c2->Print(TString(c2->GetName()) + ".png");
+
+
+
+
+}
+
+
+
+void DrawRatios() {
+
+ // Draws ratios of combined spectra
+
+ // Compute ratios
+ TH1F * hPosNegRatio[kNPart];
+
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ hPosNegRatio[ipart] = (TH1F*) hSpectra[kCombTOFTPC][ipart][kPos]->Clone();
+ hPosNegRatio[ipart]->Divide(hSpectra[kCombTOFTPC][ipart][kNeg]);
+ hPosNegRatio[ipart]->SetYTitle(TString(partLabel[ipart][kPos])+"/"+partLabel[ipart][kNeg]);
+ hPosNegRatio[ipart]->SetMinimum(0.5);
+ hPosNegRatio[ipart]->SetMaximum(1.5);
+ }
+
+ TH1F * hKPiRatio = (TH1F*) hSpectra[kCombTOFTPC][kKaon][kPos]->Clone();
+ hKPiRatio->Add(hSpectra[kCombTOFTPC][kKaon][kNeg]);
+ TH1F * htmp = (TH1F*) hSpectra[kCombTOFTPC][kPion][kPos]->Clone();
+ htmp->Add(hSpectra[kCombTOFTPC][kPion][kNeg]);
+ hKPiRatio->Divide(htmp);
+ hKPiRatio->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
+
+ TH1F * hKPiRatioMC[kNTunes];
+ for(Int_t itune = 0; itune < kNTunes; itune++){
+ hKPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kKaon][kPos]->Clone();
+ hKPiRatioMC[itune]->Add(hSpectraMC[itune][kKaon][kNeg]);
+ TH1F * htmp = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
+ htmp->Add(hSpectraMC[itune][kPion][kNeg]);
+ hKPiRatioMC[itune]->Divide(htmp);
+ hKPiRatioMC[itune]->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
+ }
+
+
+
+ TH1F * hPPiRatio = (TH1F*) hSpectra[kCombTOFTPC][kProton][kPos]->Clone();
+ hPPiRatio->Add(hSpectra[kCombTOFTPC][kProton][kNeg]);
+ htmp = (TH1F*) hSpectra[kCombTOFTPC][kPion][kPos]->Clone();
+ htmp->Add(hSpectra[kCombTOFTPC][kPion][kNeg]);
+ hPPiRatio->Divide(htmp);
+ hPPiRatio->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
+
+ TH1F * hPPiRatioMC[kNTunes];
+ for(Int_t itune = 0; itune < kNTunes; itune++){
+ hPPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kProton][kPos]->Clone();
+ hPPiRatioMC[itune]->Add(hSpectraMC[itune][kProton][kNeg]);
+ TH1F * htmp = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
+ htmp->Add(hSpectraMC[itune][kPion][kNeg]);
+ hPPiRatioMC[itune]->Divide(htmp);
+ hPPiRatioMC[itune]->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
+ }
+
+ // Draw
+// TH2F * hempty = new TH2F(TString("hempty"),"hempty",100,0.,1.5, 100, 0.001,1.8);
+// hempty->SetXTitle("p_{t} (GeV/c)");
+ // hempty->SetYTitle("");
+
+ // tmp: overlay levi fits
+ AliBWFunc * fm2 = new AliBWFunc;
+ fm2->SetVarType(AliBWFunc::kdNdpt);
+ TF1 * fLevi[kNPart] [kNCharge];
+ fLevi[kPion] [kPos] = fm2->GetLevi (mass[0], 0.1243, 7.614785, 1.524167, "fLeviPiPlus");
+ fLevi[kKaon] [kPos] = fm2->GetLevi (mass[1], 0.1625, 5.869318, 0.186361, "fLeviKPlus");
+ fLevi[kProton][kPos] = fm2->GetLevi (mass[2], 0.1773, 6.918065, 0.086389, "fLeviPPlus");
+ fLevi[kPion] [kNeg] = fm2->GetLevi (mass[0], 0.1267, 7.979582, 1.515908, "fLeviPiNeg");
+ fLevi[kKaon] [kNeg] = fm2->GetLevi (mass[1], 0.1721, 6.927956, 0.191140, "fLeviKNeg");
+ fLevi[kProton][kNeg] = fm2->GetLevi (mass[2], 0.1782, 8.160362, 0.082091, "fLeviPNeg");
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ fLevi[ipart][icharge]->SetRange(0,4);
+ }
+ }
+
+
+ for(Int_t ipart = 0; ipart < kNPart; ipart++){
+ TCanvas * c1 = new TCanvas(TString("cRatio_")+partFlag[ipart], TString("cRatio_")+partFlag[ipart]);
+ hPosNegRatio[ipart]->Draw();
+ TF1 * fRatio = new TF1 (TString("fRatio")+partFlag[ipart], TString(fLevi[ipart][kPos]->GetName())+"/"+fLevi[ipart][kNeg]->GetName());
+ // fRatio->Draw("same");
+ fRatio->SetRange(0,5);
+ if (doPrint) {
+ c1->Update();
+ gSystem->ProcessEvents();
+ c1->Print(TString(c1->GetName()) + ".eps");
+ }
+
+ }
+
+
+ TCanvas * c2 = new TCanvas(TString("cRatio_KPi"),TString("cRatio_KPi"));
+ hKPiRatio->Draw();
+ TLegend * lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
+ lMC->SetFillColor(kWhite);
+
+ // gROOT->LoadMacro("GetE735Ratios.C");
+ GetE735Ratios(0,0)->Draw("EX0,same");
+ GetE735Ratios(0,1)->Draw("EX0,same");
+ GetE735Ratios(0,2)->Draw("EX0,same");
+ GetE735Ratios(0,3)->Draw("EX0,same");
+ hKPiRatio->SetMarkerStyle(20);
+ hKPiRatio->Draw("same");
+
+ for(Int_t itune = 0; itune < kNTunes; itune++){
+ lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
+ hKPiRatioMC[itune]->SetLineWidth(2);
+ hKPiRatioMC[itune]->Draw("same,chist");
+ }
+
+ lMC->Draw();
+
+
+ TLegend * l = new TLegend( 0.1879, 0.68, 0.54,0.92);
+ l->SetFillColor(kWhite);
+ l->AddEntry(hKPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
+ l->AddEntry(GetE735Ratios(0,0), "E735, #sqrt{s} = 300 GeV","lpf");
+ l->AddEntry(GetE735Ratios(0,1), "E735, #sqrt{s} = 540 GeV","lpf");
+ l->AddEntry(GetE735Ratios(0,2), "E735, #sqrt{s} = 1000 GeV","lpf");
+ l->AddEntry(GetE735Ratios(0,3), "E735, #sqrt{s} = 1800 GeV","lpf");
+ l->Draw();
+
+
+ TCanvas * c3 = new TCanvas(TString("cRatio_PPi"),TString("cRatio_PPi"));
+ hPPiRatio->Draw();
+ hPPiRatio->SetMaximum(0.39);
+ lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
+ lMC->SetFillColor(kWhite);
+ for(Int_t itune = 0; itune < kNTunes; itune++){
+ lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
+ hKPiRatioMC[itune]->SetLineWidth(2);
+ hKPiRatioMC[itune]->Draw("same,chist");
+ }
+
+ lMC->Draw();
+
+ if (doPrint) {
+ c2->Update();
+ gSystem->ProcessEvents();
+ c2->Print(TString(c2->GetName()) + ".eps");
+ c3->Update();
+ gSystem->ProcessEvents();
+ c3->Print(TString(c3->GetName()) + ".eps");
+ }
+
+
+}
+
--- /dev/null
+TH1F* GetE735Ratios(Int_t ratio, Int_t energy) {
+
+ Int_t nbins = 10;
+
+ const Float_t binsKpi [] = {0.25,0.30, 0.35, 0.40,0.50,0.60,0.70,0.95,1.15,1.35,1.55};
+ const Float_t binsPpi [] = {0.30, 0.35, 0.40,0.50,0.60,0.70,0.95,1.15,1.35,1.55,1.75};
+ enum {kEn300,kEn540,kEn1000,kEn1800,kNenergy};
+
+
+ TH1F * hKpi[kNenergy] ;
+ TH1F * hPpi[kNenergy] ;
+
+ const Int_t markers[] = {24,25,27,28};
+
+ for(Int_t iener = 0; iener < kNenergy; iener++){
+ hKpi[iener] = new TH1F(TString("hKpi")+long(iener),TString("hKpi")+long(iener),nbins,binsKpi);
+ hPpi[iener] = new TH1F(TString("hPpi")+long(iener),TString("hPpi")+long(iener),nbins,binsPpi);
+ hKpi[iener]->SetMarkerStyle(markers[iener]);
+ hPpi[iener]->SetMarkerStyle(markers[iener]);
+ }
+
+
+
+ // 300
+ hKpi[kEn300] ->SetBinContent(0 , 8.1 /100);
+ hKpi[kEn300] ->SetBinContent(1 , 7.5 /100);
+ hKpi[kEn300] ->SetBinContent(2 , 9.5 /100);
+ hKpi[kEn300] ->SetBinContent(3 , 16.1 /100);
+ hKpi[kEn300] ->SetBinContent(4 , 16.4 /100);
+ hKpi[kEn300] ->SetBinContent(5 , 17.9 /100);
+ hKpi[kEn300] ->SetBinContent(6 , 19.4 /100);
+ hKpi[kEn300] ->SetBinContent(7 , 23.8 /100);
+ hKpi[kEn300] ->SetBinContent(8 , 25.7 /100);
+ hKpi[kEn300] ->SetBinContent(9 , 0 /100);
+ hKpi[kEn300] ->SetBinError (0 , 1.7/100);
+ hKpi[kEn300] ->SetBinError (1 , 1.4/100);
+ hKpi[kEn300] ->SetBinError (2 , 1.6/100);
+ hKpi[kEn300] ->SetBinError (3 , 1.5/100);
+ hKpi[kEn300] ->SetBinError (4 , 2.5/100);
+ hKpi[kEn300] ->SetBinError (5 , 3.5/100);
+ hKpi[kEn300] ->SetBinError (6 , 2.4/100);
+ hKpi[kEn300] ->SetBinError (7 , 4.9/100);
+ hKpi[kEn300] ->SetBinError (8 , 7.5/100);
+ hKpi[kEn300] ->SetBinError (9 , 0 /100);
+
+ hKpi[kEn540] ->SetBinContent(0 , 3.3 /100);
+ hKpi[kEn540] ->SetBinContent(1 , 12.7/100);
+ hKpi[kEn540] ->SetBinContent(2 , 10.6/100);
+ hKpi[kEn540] ->SetBinContent(3 , 13.1/100);
+ hKpi[kEn540] ->SetBinContent(4 , 16.9/100);
+ hKpi[kEn540] ->SetBinContent(5 , 17.3/100);
+ hKpi[kEn540] ->SetBinContent(6 , 20.2/100);
+ hKpi[kEn540] ->SetBinContent(7 , 21.2/100);
+ hKpi[kEn540] ->SetBinContent(8 , 25.0/100);
+ hKpi[kEn540] ->SetBinContent(9 , 31.2/100);
+ hKpi[kEn540] ->SetBinError (0 , 0.9/100);
+ hKpi[kEn540] ->SetBinError (1 , 1.1/100);
+ hKpi[kEn540] ->SetBinError (2 , 0.9/100);
+ hKpi[kEn540] ->SetBinError (3 , 0.7/100);
+ hKpi[kEn540] ->SetBinError (4 , 1.6/100);
+ hKpi[kEn540] ->SetBinError (5 , 1.7/100);
+ hKpi[kEn540] ->SetBinError (6 , 1.4/100);
+ hKpi[kEn540] ->SetBinError (7 , 2.1/100);
+ hKpi[kEn540] ->SetBinError (8 , 3.2/100);
+ hKpi[kEn540] ->SetBinError (9 , 4.9/100);
+
+ hKpi[kEn1000] ->SetBinContent(0 , 9.2 /100);
+ hKpi[kEn1000] ->SetBinContent(1 , 10.4/100);
+ hKpi[kEn1000] ->SetBinContent(2 , 9.3 /100);
+ hKpi[kEn1000] ->SetBinContent(3 , 11.8/100);
+ hKpi[kEn1000] ->SetBinContent(4 , 15.3/100);
+ hKpi[kEn1000] ->SetBinContent(5 , 19.8/100);
+ hKpi[kEn1000] ->SetBinContent(6 , 19.8/100);
+ hKpi[kEn1000] ->SetBinContent(7 , 25.9/100);
+ hKpi[kEn1000] ->SetBinContent(8 , 27.8/100);
+ hKpi[kEn1000] ->SetBinContent(9 , 32.8/100);
+ hKpi[kEn1000] ->SetBinError (0 , 0.7/100);
+ hKpi[kEn1000] ->SetBinError (1 , 0.7/100);
+ hKpi[kEn1000] ->SetBinError (2 , 0.7/100);
+ hKpi[kEn1000] ->SetBinError (3 , 0.7/100);
+ hKpi[kEn1000] ->SetBinError (4 , 1.2/100);
+ hKpi[kEn1000] ->SetBinError (5 , 1.5/100);
+ hKpi[kEn1000] ->SetBinError (6 , 1.2/100);
+ hKpi[kEn1000] ->SetBinError (7 , 1.9/100);
+ hKpi[kEn1000] ->SetBinError (8 , 2.6/100);
+ hKpi[kEn1000] ->SetBinError (9 , 3.9/100);
+
+ hKpi[kEn1800] ->SetBinContent(0 , 7.8 /100);
+ hKpi[kEn1800] ->SetBinContent(1 , 9.1 /100);
+ hKpi[kEn1800] ->SetBinContent(2 , 9.7 /100);
+ hKpi[kEn1800] ->SetBinContent(3 , 13.9/100);
+ hKpi[kEn1800] ->SetBinContent(4 , 16.2/100);
+ hKpi[kEn1800] ->SetBinContent(5 , 18.1/100);
+ hKpi[kEn1800] ->SetBinContent(6 , 21.6/100);
+ hKpi[kEn1800] ->SetBinContent(7 , 24.6/100);
+ hKpi[kEn1800] ->SetBinContent(8 , 25.8/100);
+ hKpi[kEn1800] ->SetBinContent(9 , 25.2/100);
+ hKpi[kEn1800] ->SetBinError (0 , 0.3/100);
+ hKpi[kEn1800] ->SetBinError (1 , 0.3/100);
+ hKpi[kEn1800] ->SetBinError (2 , 0.3/100);
+ hKpi[kEn1800] ->SetBinError (3 , 0.3/100);
+ hKpi[kEn1800] ->SetBinError (4 , 0.6/100);
+ hKpi[kEn1800] ->SetBinError (5 , 0.7/100);
+ hKpi[kEn1800] ->SetBinError (6 , 0.7/100);
+ hKpi[kEn1800] ->SetBinError (7 , 0.8/100);
+ hKpi[kEn1800] ->SetBinError (8 , 0.9/100);
+ hKpi[kEn1800] ->SetBinError (9 , 1.0/100);
+
+ if (ratio == 0) return hKpi[energy];
+ else return hPpi[energy];
+
+}
--- /dev/null
+TGraphErrors ** StarPPSpectra() {
+
+ TGraphErrors ** gStar = new TGraphErrors*[6];
+ const Int_t kNMaxPoints = 19;
+ Double_t pt[kNMaxPoints];
+ Double_t errpt[kNMaxPoints] = {0};
+ Double_t pim [kNMaxPoints];
+ Double_t pime[kNMaxPoints];
+ Double_t pip [kNMaxPoints];
+ Double_t pipe[kNMaxPoints];
+ Double_t km [kNMaxPoints];
+ Double_t kme[kNMaxPoints];
+ Double_t kp [kNMaxPoints];
+ Double_t kpe[kNMaxPoints];
+ Double_t pm [kNMaxPoints];
+ Double_t pme[kNMaxPoints];
+ Double_t pp [kNMaxPoints];
+ Double_t ppe[kNMaxPoints];
+
+ Int_t ipoint = 0;
+ pt[ipoint++] = 0.225;
+ pt[ipoint++] = 0.275;
+ pt[ipoint++] = 0.325;
+ pt[ipoint++] = 0.375;
+ pt[ipoint++] = 0.425;
+ pt[ipoint++] = 0.475;
+ pt[ipoint++] = 0.525;
+ pt[ipoint++] = 0.575;
+ pt[ipoint++] = 0.625;
+ pt[ipoint++] = 0.675;
+ pt[ipoint++] = 0.725;
+ pt[ipoint++] = 0.775;
+ pt[ipoint++] = 0.825;
+ pt[ipoint++] = 0.875;
+ pt[ipoint++] = 0.925;
+ pt[ipoint++] = 0.975;
+ pt[ipoint++] = 1.025;
+ pt[ipoint++] = 1.075;
+ pt[ipoint++] = 1.125;
+
+
+ // pi-
+ ipoint = 0;
+ pim[ipoint++] = 2.02;
+ pim[ipoint++] = 1.52;
+ pim[ipoint++] = 1.13;
+ pim[ipoint++] = 8.44*1e-1;
+ pim[ipoint++] = 6.35*1e-1;
+ pim[ipoint++] = 4.69*1e-1;
+ pim[ipoint++] = 3.54*1e-1;
+ pim[ipoint++] = 2.67*1e-1;
+ pim[ipoint++] = 2.02*1e-1;
+ pim[ipoint++] = 1.53*1e-1;
+ pim[ipoint++] = 1.16*1e-1;
+ pim[ipoint++] = 8.65*1e-2;
+
+ ipoint = 0;
+ pime[ipoint++] = 0.06;
+ pime[ipoint++] = 0.03;
+ pime[ipoint++] = 0.02;
+ pime[ipoint++] = 0.09*1e-1;
+ pime[ipoint++] = 0.07*1e-1;
+ pime[ipoint++] = 0.05*1e-1;
+ pime[ipoint++] = 0.04*1e-1;
+ pime[ipoint++] = 0.03*1e-1;
+ pime[ipoint++] = 0.04*1e-1;
+ pime[ipoint++] = 0.03*1e-1;
+ pime[ipoint++] = 0.03*1e-1;
+ pime[ipoint++] = 0.30*1e-2;
+
+ gStar[0] = new TGraphErrors(ipoint,pt,pim,errpt,pime);
+
+ // pi+
+ ipoint = 0;
+ pip[ipoint++] = 2.07;
+ pip[ipoint++] = 1.54;
+ pip[ipoint++] = 1.14;
+ pip[ipoint++] = 8.57*1e-1;
+ pip[ipoint++] = 6.38*1e-1;
+ pip[ipoint++] = 4.76*1e-1;
+ pip[ipoint++] = 3.59*1e-1;
+ pip[ipoint++] = 2.73*1e-1;
+ pip[ipoint++] = 2.07*1e-1;
+ pip[ipoint++] = 1.55*1e-1;
+ pip[ipoint++] = 1.16*1e-1;
+ pip[ipoint++] = 8.95*1e-2;
+
+ ipoint = 0;
+ pipe[ipoint++] = 0.06;
+ pipe[ipoint++] = 0.03;
+ pipe[ipoint++] = 0.02;
+ pipe[ipoint++] = 0.09*1e-1;
+ pipe[ipoint++] = 0.07*1e-1;
+ pipe[ipoint++] = 0.05*1e-1;
+ pipe[ipoint++] = 0.04*1e-1;
+ pipe[ipoint++] = 0.03*1e-1;
+ pipe[ipoint++] = 0.04*1e-1;
+ pipe[ipoint++] = 0.03*1e-1;
+ pipe[ipoint++] = 0.03*1e-1;
+ pipe[ipoint++] = 0.31*1e-2;
+
+ gStar[1] = new TGraphErrors(ipoint,pt,pip,errpt,pipe);
+
+ // k-
+ ipoint = 0;
+ km[ipoint++] = 1.43*1e-1;
+ km[ipoint++] = 1.26*1e-1;
+ km[ipoint++] = 1.08*1e-1;
+ km[ipoint++] = 8.77*1e-2;
+ km[ipoint++] = 7.34*1e-2;
+ km[ipoint++] = 6.17*1e-2;
+ km[ipoint++] = 4.87*1e-2;
+ km[ipoint++] = 4.11*1e-2;
+ km[ipoint++] = 3.70*1e-2;
+ km[ipoint++] = 2.86*1e-2;
+ km[ipoint++] = 2.42*1e-2;
+
+ ipoint = 0;
+ kme[ipoint++] = 0.11*1e-1;
+ kme[ipoint++] = 0.05*1e-1;
+ kme[ipoint++] = 0.02*1e-1;
+ kme[ipoint++] = 0.24*1e-2;
+ kme[ipoint++] = 0.26*1e-2;
+ kme[ipoint++] = 0.63*1e-2;
+ kme[ipoint++] = 0.50*1e-2;
+ kme[ipoint++] = 0.42*1e-2;
+ kme[ipoint++] = 0.28*1e-2;
+ kme[ipoint++] = 0.22*1e-2;
+ kme[ipoint++] = 0.27*1e-2;
+
+ gStar[2] = new TGraphErrors(ipoint,pt,km,errpt,kme);
+
+
+ // k+
+ ipoint = 0;
+ kp[ipoint++] = 1.52*1e-1;
+ kp[ipoint++] = 1.30*1e-1;
+ kp[ipoint++] = 1.08*1e-1;
+ kp[ipoint++] = 9.16*1e-2;
+ kp[ipoint++] = 7.47*1e-2;
+ kp[ipoint++] = 6.26*1e-2;
+ kp[ipoint++] = 5.26*1e-2;
+ kp[ipoint++] = 4.41*1e-2;
+ kp[ipoint++] = 3.81*1e-2;
+ kp[ipoint++] = 3.06*1e-2;
+ kp[ipoint++] = 2.49*1e-2;
+
+ ipoint = 0;
+ kpe[ipoint++] = 0.11*1e-1;
+ kpe[ipoint++] = 0.05*1e-1;
+ kpe[ipoint++] = 0.02*1e-1;
+ kpe[ipoint++] = 0.25*1e-2;
+ kpe[ipoint++] = 0.27*1e-2;
+ kpe[ipoint++] = 0.64*1e-2;
+ kpe[ipoint++] = 0.54*1e-2;
+ kpe[ipoint++] = 0.45*1e-2;
+ kpe[ipoint++] = 0.28*1e-2;
+ kpe[ipoint++] = 0.24*1e-2;
+ kpe[ipoint++] = 0.28*1e-2;
+
+ gStar[3] = new TGraphErrors(ipoint,pt,kp,errpt,kpe);
+
+
+ // pbar
+ ipoint = 0;
+ pm[ipoint++] = 0;
+ pm[ipoint++] = 0;
+ pm[ipoint++] = 0;
+ pm[ipoint++] = 5.54*1e-2;
+ pm[ipoint++] = 4.81*1e-2;
+ pm[ipoint++] = 4.25*1e-2;
+ pm[ipoint++] = 3.77*1e-2;
+ pm[ipoint++] = 3.38*1e-2;
+ pm[ipoint++] = 2.78*1e-2;
+ pm[ipoint++] = 2.49*1e-2;
+ pm[ipoint++] = 2.03*1e-2;
+ pm[ipoint++] = 1.75*1e-2;
+ pm[ipoint++] = 1.52*1e-2;
+ pm[ipoint++] = 1.27*1e-2;
+ pm[ipoint++] = 1.05*1e-2;
+ pm[ipoint++] = 8.95*1e-3;
+ pm[ipoint++] = 7.36*1e-3;
+ pm[ipoint++] = 6.59*1e-3;
+ pm[ipoint++] = 5.21*1e-3;
+
+ ipoint = 0;
+ pme[ipoint++] = 0;
+ pme[ipoint++] = 0;
+ pme[ipoint++] = 0;
+ pme[ipoint++] = 0.13*1e-2;
+ pme[ipoint++] = 0.11*1e-2;
+ pme[ipoint++] = 0.10*1e-2;
+ pme[ipoint++] = 0.09*1e-2;
+ pme[ipoint++] = 0.08*1e-2;
+ pme[ipoint++] = 0.07*1e-2;
+ pme[ipoint++] = 0.06*1e-2;
+ pme[ipoint++] = 0.06*1e-2;
+ pme[ipoint++] = 0.06*1e-2;
+ pme[ipoint++] = 0.05*1e-2;
+ pme[ipoint++] = 0.04*1e-2;
+ pme[ipoint++] = 0.04*1e-2;
+ pme[ipoint++] = 0.34*1e-3;
+ pme[ipoint++] = 0.34*1e-3;
+ pme[ipoint++] = 0.32*1e-3;
+ pme[ipoint++] = 0.29*1e-3;
+
+ gStar[4] = new TGraphErrors(ipoint,pt,pm,errpt,pme);
+
+
+ // p
+ ipoint = 0;
+ pp[ipoint++] = 0;
+ pp[ipoint++] = 0;
+ pp[ipoint++] = 0;
+ pp[ipoint++] = 0;
+ pp[ipoint++] = 0;
+ pp[ipoint++] = 5.07*1e-2;
+ pp[ipoint++] = 4.69*1e-2;
+ pp[ipoint++] = 4.08*1e-2;
+ pp[ipoint++] = 3.42*1e-2;
+ pp[ipoint++] = 2.87*1e-2;
+ pp[ipoint++] = 2.41*1e-2;
+ pp[ipoint++] = 2.13*1e-2;
+ pp[ipoint++] = 1.82*1e-2;
+ pp[ipoint++] = 1.54*1e-2;
+ pp[ipoint++] = 1.31*1e-2;
+ pp[ipoint++] = 1.11*1e-2;
+ pp[ipoint++] = 9.78*1e-3;
+ pp[ipoint++] = 8.56*1e-3;
+ pp[ipoint++] = 7.38*1e-3;
+
+ ipoint = 0;
+ ppe[ipoint++] = 0;
+ ppe[ipoint++] = 0;
+ ppe[ipoint++] = 0;
+ ppe[ipoint++] = 0;
+ ppe[ipoint++] = 0;
+ ppe[ipoint++] = 0.25*1e-2;
+ ppe[ipoint++] = 0.19*1e-2;
+ ppe[ipoint++] = 0.14*1e-2;
+ ppe[ipoint++] = 0.10*1e-2;
+ ppe[ipoint++] = 0.07*1e-2;
+ ppe[ipoint++] = 0.07*1e-2;
+ ppe[ipoint++] = 0.06*1e-2;
+ ppe[ipoint++] = 0.05*1e-2;
+ ppe[ipoint++] = 0.05*1e-2;
+ ppe[ipoint++] = 0.04*1e-2;
+ ppe[ipoint++] = 0.04*1e-2;
+ ppe[ipoint++] = 0.40*1e-3;
+ ppe[ipoint++] = 0.37*1e-3;
+ ppe[ipoint++] = 0.38*1e-3;
+
+
+
+ gStar[5] = new TGraphErrors(ipoint,pt,pp,errpt,ppe);
+
+ gStar[0]->SetMarkerStyle(kOpenTriangleUp);
+ gStar[1]->SetMarkerStyle(kFullTriangleUp);
+ gStar[2]->SetMarkerStyle(kOpenCircle);
+ gStar[3]->SetMarkerStyle(kFullCircle);
+ gStar[4]->SetMarkerStyle(kOpenSquare);
+ gStar[5]->SetMarkerStyle(kFullSquare);
+
+ gStar[0]->SetMarkerColor(kBlack);
+ gStar[1]->SetMarkerColor(kBlack);
+ gStar[2]->SetMarkerColor(kRed);
+ gStar[3]->SetMarkerColor(kRed);
+ gStar[4]->SetMarkerColor(kBlue);
+ gStar[5]->SetMarkerColor(kBlue);
+
+ gStar[0]->SetLineColor (kBlack);
+ gStar[1]->SetLineColor (kBlack);
+ gStar[2]->SetLineColor (kRed);
+ gStar[3]->SetLineColor (kRed);
+ gStar[4]->SetLineColor (kBlue);
+ gStar[5]->SetLineColor (kBlue);
+
+ // Scale all for 2*pi and pt in order to go to dNdpt
+ for(Int_t istar = 0; istar < 6; istar++){
+ Int_t npoint = gStar[istar]->GetN();
+ for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+ Double_t ptbin = gStar[istar]->GetX()[ipoint];
+ gStar[istar]->SetPoint(ipoint,
+ gStar[istar]->GetX()[ipoint],
+ gStar[istar]->GetY()[ipoint]*2*TMath::Pi()*ptbin);
+ gStar[istar]->SetPointError(ipoint,
+ gStar[istar]->GetEX()[ipoint],
+ gStar[istar]->GetEY()[ipoint]*2*TMath::Pi()*ptbin);
+ }
+
+ }
+
+
+
+// gStar[0]->Draw("AP");
+// gStar[1]->Draw("P");
+// gStar[2]->Draw("P");
+// gStar[3]->Draw("P");
+// gStar[4]->Draw("P");
+// gStar[5]->Draw("P");
+
+
+ return gStar;
+
+
+
+ // table from PRC79 79 (2009) 34909
+// pi- pi+ k- K+ pbar p
+// 0.225 2.02+-0.06 2.07+-0.06 (1.43+-0.11)*10-1 (1.52+-0.11)*10-1
+// 0.275 1.52+-0.03 1.54+-0.03 (1.26+-0.05)*10-1 (1.30+-0.05)*10-1
+// 0.325 1.13+-0.02 1.14+-0.02 (1.08+-0.02)*10-1 (1.08+-0.02)*10-1
+// 0.375 (8.44+-0.09)*10-1 (8.57+-0.09)*10-1 (8.77+-0.24)*10-2 (9.16+-0.25)*10-2 (5.54+-0.13)*10-2
+// 0.425 (6.35+-0.07)*10-1 (6.38+-0.07)*10-1 (7.34+-0.26)*10-2 (7.47+-0.27)*10-2 (4.81+-0.11)*10-2
+// 0.475 (4.69+-0.05)*10-1 (4.76+-0.05)*10-1 (6.17+-0.63)*10-2 (6.26+-0.64)*10-2 (4.25+-0.10)*10-2 (5.07+-0.25)*10-2
+// 0.525 (3.54+-0.04)*10-1 (3.59+-0.04)*10-1 (4.87+-0.50)*10-2 (5.26+-0.54)*10-2 (3.77+-0.09)*10-2 (4.69+-0.19)*10-2
+// 0.575 (2.67+-0.03)*10-1 (2.73+-0.03)*10-1 (4.11+-0.42)*10-2 (4.41+-0.45)*10-2 (3.38+-0.08)*10-2 (4.08+-0.14)*10-2
+// 0.625 (2.02+-0.04)*10-1 (2.07+-0.04)*10-1 (3.70+-0.28)*10-2 (3.81+-0.28)*10-2 (2.78+-0.07)*10-2 (3.42+-0.10)*10-2
+// 0.675 (1.53+-0.03)*10-1 (1.55+-0.03)*10-1 (2.86+-0.22)*10-2 (3.06+-0.24)*10-2 (2.49+-0.06)*10-2 (2.87+-0.07)*10-2
+// 0.725 (1.16+-0.03)*10-1 (1.16+-0.03)*10-1 (2.42+-0.27)*10-2 (2.49+-0.28)*10-2 (2.03+-0.06)*10-2 (2.41+-0.07)*10-2
+// 0.775 (8.65+-0.30)*10-2 (8.95+-0.31)*10-2 (1.75+-0.06)*10-2 (2.13+-0.06)*10-2
+// 0.825 (1.52+-0.05)*10-2 (1.82+-0.05)*10-2
+// 0.875 (1.27+-0.04)*10-2 (1.54+-0.05)*10-2
+// 0.925 (1.05+-0.04)*10-2 (1.31+-0.04)*10-2
+// 0.975 (8.95+-0.34)*10-3 (1.11+-0.04)*10-2
+// 1.025 (7.36+-0.34)*10-3 (9.78+-0.40)*10-3
+// 1.075 (6.59+-0.32)*10-3 (8.56+-0.37)*10-3
+// 1.125 (5.21+-0.29)*10-3 (7.38+-0.38)*10-3
+
+
+ }
SPECTRA/AliAnalysisCentralCutEvtMC.cxx \
SPECTRA/AliAnalysisCentralCutEvtESD.cxx \
SPECTRA/AliAnalysisCentralExtrapolate.cxx \
- SPECTRA/AliAnalysisTaskCentral.cxx
+ SPECTRA/AliAnalysisTaskCentral.cxx \
+ SPECTRA/Fit/AliBWTools.cxx \
+ SPECTRA/Fit/AliBWFunc.cxx \
+ SPECTRA/Fit/AliLatexTable.cxx
HDRS= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff