//---------------------------------------------------------------------- // AliPWGFunc // // 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 AliPWGFunc::SetVarType with one of // the elements of the VarType_t enum. // // TODO: code cleaup, make the naming convention of function more transparent // // Warning: not all variables are implemented for all the functions. // // Author: M. Floris, CERN //---------------------------------------------------------------------- #include "AliPWGFunc.h" #include "TMath.h" #include "TF1.h" #include "TF3.h" #include "TH1.h" #include "TSpline.h" #include "AliLog.h" ClassImp(AliPWGFunc) AliPWGFunc::AliPWGFunc () : fLastFunc(0), fLineWidth(1), fVarType(kdNdpt) { // ctor fLineWidth = 1; } AliPWGFunc::~AliPWGFunc(){ // dtor if (fLastFunc) delete fLastFunc; } TF1 * AliPWGFunc::GetHistoFunc(TH1 * h, const char * name) { // Regardless of the variable type, this returns a function made // from the histo * a multiplicative normalization. // This uses a bad hack... 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 * AliPWGFunc::GetGraphFunc(TGraph * g, const char * name) { // Regardless of the variable type, this returns a function made // from the graph * a multiplicative normalization. // This uses a bad hack... fLastFunc = new TF1 (name, StaticHistoFunc, 0.0, 10, 2); fLastFunc->SetParameter(0,1); fLastFunc->FixParameter(1,Double_t(Long64_t(g))); fLastFunc->SetParNames("norm", "pointer to histo"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetBGBW(Double_t mass, Double_t beta, Double_t T, Double_t n, Double_t norm, const char * name){ // Boltzmann-Gibbs blast wave switch (fVarType) { case kdNdpt: return GetBGBWdNdptTimesPt(mass,beta,T,n,norm,name); break; case kOneOverPtdNdpt: return GetBGBWdNdpt(mass,beta,T,n,norm,name); break; case kOneOverMtdNdmt: AliFatal("Not implemented"); break; default: AliFatal("Not implemented"); } return 0; } TF1 * AliPWGFunc::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 * AliPWGFunc::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 * AliPWGFunc::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: return GetMTExpdNdmt(mass,T,norm,name,kOneOverMtdNdmt); break; case kdNdmt: return GetMTExpdNdmt(mass,T,norm,name,kdNdmt); break; case kOneOverMtdNdmtMinusM: return GetMTExpdNdmt(mass,T,norm,name,kOneOverMtdNdmtMinusM); break; default: AliFatal("Not implemented"); } return 0; } TF1 * AliPWGFunc::GetBoseEinstein(Double_t mass, Double_t T, Double_t norm, const char * name){ // Bose einstein switch (fVarType) { case kdNdpt: return GetBoseEinsteindNdptTimesPt(mass,T,norm,name); break; case kOneOverPtdNdpt: return GetBoseEinsteindNdpt(mass,T,norm,name); break; case kOneOverMtdNdmt: AliFatal("Not implemented"); break; default: AliFatal("Not implemented"); } return 0; } TF1 * AliPWGFunc::GetFermiDirac(Double_t mass, Double_t T, Double_t norm, const char * name){ // Simple exponential in 1/mt*MT switch (fVarType) { case kdNdpt: return GetFermiDiracdNdptTimesPt(mass,T,norm,name); break; case kOneOverPtdNdpt: return GetFermiDiracdNdpt(mass,T,norm,name); break; case kOneOverMtdNdmt: AliFatal("Not implemented"); break; default: AliFatal("Not implemented"); } return 0; } TF1 * AliPWGFunc::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 * AliPWGFunc::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,kOneOverMtdNdmt); break; case kdNdmt: return GetLevidNdmt(mass,T,n,norm,name,kdNdmt); break; case kOneOverMtdNdmtMinusM: return GetLevidNdmt(mass,T,n,norm,name,kOneOverMtdNdmtMinusM); break; default: AliFatal("Not implemented"); } return 0; } TF1 * AliPWGFunc::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 * AliPWGFunc::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: fLastFunc = new TF1 (name, StaticUA1FuncOneOverPt, 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 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 AliPWGFunc::StaticHistoFunc(const double * x, const 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]; TObject * h = (TObject*) Long64_t(p[1]); // Int_t bin = h->FindBin(x[0]); // double value = h->GetBinContent(bin); // static TH1 * oldptr = 0; // static TSpline3 * spl = 0; // if (h!=oldptr) { // FIXME: recheck static pointers TSpline3 * spl = 0; if(h->InheritsFrom("TH1")) { if ( ((TH1*)h)->FindBin(x[0]) > ((TH1*)h)->GetNbinsX()) return 0; spl= new TSpline3((TH1*)h); } else if(h->InheritsFrom("TGraph")) spl= new TSpline3("fGraph",(TGraph*)h); else { Printf("AliPWGFunc::StaticHistoFunc: Unsupported type"); return 0; } // } double value = spl->Eval(x[0]); delete spl; return value*norm; } Double_t AliPWGFunc::StaticUA1Func(const double * x, const 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 temp = p[4]; Double_t norm = p[5]; Double_t xx = x[0]; static AliPWGFunc * self = new AliPWGFunc; static TF1 * fPLaw = self->GetPowerLawdNdptTimesPt(pt0, n, norm, "fLocalPLawUA1"); static TF1 * fPMTExp = self->GetMTExpdNdptTimesPt (mass, temp, norm, "fLocalMTexpUA1"); fPLaw->SetParameters(norm,pt0,n); fPMTExp->SetParameters(1,temp); 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("AliPWGFunc::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 AliPWGFunc::StaticUA1FuncOneOverPt(const double * x, const double* p) { // UA1 func // "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 temp = p[4]; Double_t norm = p[5]; Double_t xx = x[0]; static AliPWGFunc * self = new AliPWGFunc; static TF1 * fPLaw = self->GetPowerLawdNdpt(pt0, n, norm, "fLocalPLawUA1"); static TF1 * fPMTExp = self->GetMTExpdNdpt (mass, temp, norm, "fLocalMTexpUA1"); fPLaw->SetParameters(norm,pt0,n); fPMTExp->SetParameters(1,temp); 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("AliPWGFunc::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 AliPWGFunc::IntegrandBG(const double * x, const double* p){ // integrand for boltzman-gibbs blast wave // x[0] -> r (radius) // p[0] -> mass // p[1] -> pT (transverse momentum) // p[2] -> beta_max (surface velocity) // p[3] -> T (freezout temperature) // p[4] -> n (velocity profile) double x0 = x[0]; double mass = p[0]; double pT = p[1]; double beta_max = p[2]; double temp = p[3]; Double_t n = p[4]; // Keep beta within reasonable limits Double_t beta = beta_max * TMath::Power(x0, n); if (beta > 0.9999999999999999) beta = 0.9999999999999999; double mT = TMath::Sqrt(mass*mass+pT*pT); double rho0 = TMath::ATanH(beta); double arg00 = pT*TMath::SinH(rho0)/temp; if (arg00 > 700.) arg00 = 700.; // avoid FPE double arg01 = mT*TMath::CosH(rho0)/temp; double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01); // printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", x0, pT, beta_max, temp, n, mT, beta, rho0, arg00, arg01); return f0; } Double_t AliPWGFunc::StaticBGdNdPt(const double * x, const double* p) { // implementation of BGBW (1/pt dNdpt) double pT = x[0];; double mass = p[0]; double beta = p[1]; double temp = p[2]; double n = p[3]; double norm = p[4]; static TF1 * fIntBG = 0; if(!fIntBG) fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 5); fIntBG->SetParameters(mass, pT, beta, temp,n); double result = fIntBG->Integral(0,1); // printf ("[%4.4f], Int :%f\n", pT, result); return result*norm;//*1e30;; } Double_t AliPWGFunc::StaticBGdNdPtTimesPt(const double * x, const double* p) { // BGBW dNdpt implementation return x[0]*StaticBGdNdPt(x,p); } Double_t AliPWGFunc::StaticBGdNdMtTimesMt(const double * x, const double* p) { // BGBW dNdpt implementation // X0 is mt here Double_t pt = TMath::Sqrt(x[0]*x[0]-p[0]*p[0]); return pt*StaticBGdNdPt(&pt,p); } TF1 * AliPWGFunc::GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t temp, Double_t n, Double_t norm, const char * name){ // BGBW 1/pt dNdpt fLastFunc = new TF1 (name, StaticBGdNdPt, 0.0, 10, 5); fLastFunc->SetParameters(mass,beta,temp,n,norm); fLastFunc->FixParameter(0,mass); fLastFunc->SetParNames("mass", "#beta", "T", "n", "norm"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } //_____________________________________________________________________ // Tsallis Double_t AliPWGFunc::IntegrandTsallis(const double * x, const 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 temp = 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)/temp * ( mt*TMath::CosH(y)*TMath::CosH(rho) - pt*TMath::SinH(rho)*TMath::Cos(phi) ) ), -1/(q-1) ); return res; } Double_t AliPWGFunc::StaticTsallisdNdPt(const double * x, const double* p) { // tsallis BW implementation 1/pt dNdpt double pT = x[0];; double mass = p[0]; double beta = p[1]; double temp = 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, temp, 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 AliPWGFunc::StaticTsallisdNdPtTimesPt(const double * x, const double* p) { // tsallis BW , implementatio of dNdpt return x[0]*StaticTsallisdNdPt(x,p); } TF1 * AliPWGFunc::GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t temp, 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,temp,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", "temp", "q", "norm", "ymax"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } // Times Pt funcs // Boltzmann-Gibbs Blast Wave TF1 * AliPWGFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t n, Double_t norm, const char * name){ // BGBW, dNdpt fLastFunc = new TF1 (name, StaticBGdNdPtTimesPt, 0.0, 10, 5); fLastFunc->SetParameters(mass,beta,temp,n,norm); fLastFunc->FixParameter(0,mass); fLastFunc->SetParNames("mass", "#beta", "temp", "n", "norm"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } // Boltzmann-Gibbs Blast Wave TF1 * AliPWGFunc::GetBGBWdNdptTimesMt(Double_t mass, Double_t beta, Double_t temp, Double_t n, Double_t norm, const char * name){ // BGBW, dNdpt // 1/Mt dN/dmt fLastFunc = new TF1 (name, StaticBGdNdMtTimesMt, 0.0, 10, 5); fLastFunc->SetParameters(mass,beta,temp,n,norm); fLastFunc->FixParameter(0,mass); fLastFunc->SetParNames("mass", "#beta", "temp", "n", "norm"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, 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,temp,q,norm,ymax); fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetMTExpdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){ // Simple exponential in 1/mt*MT, as a function of dNdpt char formula[500]; snprintf(formula,500,"[0]*x*exp(-sqrt(x**2+%f**2)/[1])", mass); fLastFunc=new TF1(name,formula,0,10); fLastFunc->SetParameters(norm, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetBoseEinsteindNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){ // Bose einstein distribution as a function of dNdpt char formula[500]; snprintf(formula,500,"[0]*x*1./(exp(sqrt(x**2+%f**2)/[1])-1)", mass); fLastFunc=new TF1(name,formula,0,10); fLastFunc->SetParameters(norm, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetFermiDiracdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){ // Bose einstein distribution as a function of dNdpt char formula[500]; snprintf(formula,500,"[0]*x*1./(exp(sqrt(x**2+%f**2)/[1])+1)", mass); fLastFunc=new TF1(name,formula,0,10); fLastFunc->SetParameters(norm, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetPTExpdNdptTimesPt(Double_t temp, Double_t norm, const char * name){ // Simple exponential in 1/pt*dNdpT, as a function of dNdpt char formula[500]; snprintf(formula,500,"[0]*x*exp(-x/[1])"); fLastFunc=new TF1(name,formula,0,10); fLastFunc->SetParameters(norm, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetBoltzmanndNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){ // Boltzmann (exp in 1/mt*dNdmT times mt) as a function of dNdpt char formula[500]; snprintf(formula,500,"[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, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } // Tsallis (no BW, a la CMS) // TF1 * AliPWGFunc::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 * AliPWGFunc::GetLevidNdptTimesPt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){ // Levi function, dNdpt char formula[500]; snprintf(formula,500,"( 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, temp,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 * AliPWGFunc::GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name){ // PowerLaw function, dNdpt char formula[500]; snprintf(formula,500,"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 * AliPWGFunc::GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name){ // PowerLaw function, 1/pt dNdpt char formula[500]; snprintf(formula,500," [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 * AliPWGFunc::GetLevidNdpt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){ // Levi function, dNdpt char formula[500]; snprintf(formula,500,"( [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, temp,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 * AliPWGFunc::GetLevidNdmt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name, VarType_t var){ // Levi function, 1/mt dNdmt char formula[500]; if (var == kOneOverMtdNdmt) snprintf(formula,500,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x -[3])/([1]*[2]) )^(-[1])"); else if (var == kdNdmt) snprintf(formula,500,"( x*[0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x-[3])/([1]*[2]) )^(-[1])"); if (var == kOneOverMtdNdmtMinusM) snprintf(formula,500,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x)/([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, temp,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 AliPWGFunc::IntegrandTest(const double * x, const double* p){ // test function Double_t y = x[0]; Double_t mass = p[0]; Double_t pt = p[1]; Double_t temp = p[2]; Double_t mt = TMath::Sqrt(mass*mass+pt*pt); return mt*TMath::CosH(y)*TMath::Exp(-mt*TMath::CosH(y)/temp); } Double_t AliPWGFunc::StaticTest(const double * x, const double* p) { // test function double pT = x[0];; double mass = p[0]; double temp = 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, temp); double result = fIntTest->Integral(-ymax, ymax); return result*p[2];//*1e30;; } TF1 * AliPWGFunc::GetTestFunc(Double_t mass, Double_t temp, Double_t norm, Double_t ymax, const char * name){ // test function fLastFunc = new TF1 (name, StaticTest, 0.0, 10, 4); fLastFunc->SetParameters(mass,temp,norm,ymax); fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } //___________________________________________________________ TF1 * AliPWGFunc::GetMTExpdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){ // Simple exp in 1/mt dNdmt, as a function of dNdpt // mt scaling char formula[500]; snprintf(formula,500,"[0]*exp(-sqrt(x**2+%f**2)/[1])", mass); fLastFunc=new TF1(name,formula,0,10); fLastFunc->SetParameters(norm, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetMTExpdNdmt(Double_t mass, Double_t temp, Double_t norm, const char * name, VarType_t var){ // Levi function, 1/mt dNdmt1/mt dNdmt, char formula[500]; if (var == kOneOverMtdNdmt) snprintf(formula,500,"[0] * exp (-x/[1]) + %f ", mass); else if (var == kdNdmt) snprintf(formula,500,"[0] * x * exp (-x/[1]) + %f ", mass); if (var == kOneOverMtdNdmtMinusM) snprintf(formula,500,"[0] * exp (-x/[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, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetBoseEinsteindNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){ // bose einstein char formula[500]; snprintf(formula,500,"[0]*1./(exp(sqrt(x**2+%f**2)/[1])-1)", mass); fLastFunc=new TF1(name,formula,0,10); fLastFunc->SetParameters(norm, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } TF1 * AliPWGFunc::GetFermiDiracdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){ // bose einstein char formula[500]; snprintf(formula,500,"[0]*1./(exp(sqrt(x**2+%f**2)/[1])+1)", mass); fLastFunc=new TF1(name,formula,0,10); fLastFunc->SetParameters(norm, temp); fLastFunc->SetParLimits(1, 0.01, 10); fLastFunc->SetParNames("norm", "T"); fLastFunc->SetLineWidth(fLineWidth); return fLastFunc; } // // Simple tsallis (a la CMS) // TF1 * AliPWGFunc::GetTsallisdNdpt(Double_t mass, Double_t temp, 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, temp, q); // fLastFunc->SetParLimits(1, 0.01, 10); // fLastFunc->SetParNames("norm", "T", "q"); // fLastFunc->SetLineWidth(fLineWidth); // return fLastFunc; // }