PWG/Tools/AliPWGFunc.cxx

   1 //----------------------------------------------------------------------
   2 //                              AliPWGFunc
   3 //
   4 // This class implements several function useful to fit pt spectra,
   5 // including but not limited to blast wave models.
   6 //
   7 // It can return the same functional for as a function of different
   8 // variables: dNdpt vs pt, 1/pt dNdpt vs pt, 1/mt dNdmt vs mt.
   9 //
  10 // Before getting the function you need, you have to chose the
  11 // variable you want to use calling AliPWGFunc::SetVarType with one of
  12 // the elements of the VarType_t enum.
  13 //
  14 // TODO: code cleaup, make the naming convention of function more transparent
  15 //
  16 // Warning: not all variables are implemented for all the functions.
  17 //
  18 // Author: M. Floris, CERN
  19 //----------------------------------------------------------------------
  20
  21 #include "AliPWGFunc.h"
  22 #include "TMath.h"
  23 #include "TF1.h"
  24 #include "TF3.h"
  25 #include "TH1.h"
  26 #include "TSpline.h"
  27 #include "AliLog.h"
  28
  29 ClassImp(AliPWGFunc)
  30
  31 AliPWGFunc::AliPWGFunc () : fLastFunc(0),  fLineWidth(1), fVarType(kdNdpt) {
  32
  33   // ctor
  34   fLineWidth = 1;
  35 }
  36 AliPWGFunc::~AliPWGFunc(){
  37
  38   // dtor
  39   if (fLastFunc) delete fLastFunc;
  40
  41 }
  42
  43
  44 TF1 * AliPWGFunc::GetHistoFunc(TH1 * h, const char * name) {
  45
  46   // Regardless of the variable type, this returns a function made
  47   // from the histo * a multiplicative normalization.
  48   // This uses a bad hack...
  49
  50   fLastFunc = new TF1 (name, StaticHistoFunc, 0.0, 10, 2);
  51   fLastFunc->SetParameter(0,1);
  52   fLastFunc->FixParameter(1,Double_t(Long64_t(h)));
  53   fLastFunc->SetParNames("norm", "pointer to histo");
  54   fLastFunc->SetLineWidth(fLineWidth);
  55   return fLastFunc;
  56
  57
  58
  59 }
  60 TF1 * AliPWGFunc::GetGraphFunc(TGraph * g, const char * name) {
  61
  62   // Regardless of the variable type, this returns a function made
  63   // from the graph * a multiplicative normalization.
  64   // This uses a bad hack...
  65
  66   fLastFunc = new TF1 (name, StaticHistoFunc, 0.0, 10, 2);
  67   fLastFunc->SetParameter(0,1);
  68   fLastFunc->FixParameter(1,Double_t(Long64_t(g)));
  69   fLastFunc->SetParNames("norm", "pointer to histo");
  70   fLastFunc->SetLineWidth(fLineWidth);
  71   return fLastFunc;
  72
  73
  74 }
  75
  76
  77 TF1 * AliPWGFunc::GetBGBW(Double_t mass, Double_t beta, Double_t T,
  78                          Double_t n, Double_t norm, const char * name){
  79
  80   // Boltzmann-Gibbs blast wave
  81
  82   switch (fVarType) {
  83   case kdNdpt:
  84     return GetBGBWdNdptTimesPt(mass,beta,T,n,norm,name);
  85     break;
  86   case kOneOverPtdNdpt:
  87     return GetBGBWdNdpt(mass,beta,T,n,norm,name);
  88     break;
  89   case kOneOverMtdNdmt:
  90     AliFatal("Not implemented");
  91     break;
  92   default:
  93     AliFatal("Not implemented");
  94   }
  95
  96   return 0;
  97
  98 }
  99
 100
 101 TF1 * AliPWGFunc::GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char * name){
 102   // Boltzmann
 103   switch (fVarType) {
 104   case kdNdpt:
 105     return GetBoltzmanndNdptTimesPt(mass, T, norm, name);
 106   case kOneOverPtdNdpt:
 107     AliFatal("Not implemented");
 108     break;
 109   case kOneOverMtdNdmt:
 110     AliFatal("Not implemented");
 111     break;
 112   default:
 113     AliFatal("Not implemented");
 114   }
 115
 116   return 0;
 117
 118 }
 119
 120
 121 TF1 * AliPWGFunc::GetTsallisBW(Double_t mass, Double_t beta, Double_t T, Double_t q,
 122                               Double_t norm, Double_t ymax, const char * name){
 123
 124   // Tsallis blast wave
 125   switch (fVarType) {
 126   case kdNdpt:
 127     return GetTsallisBWdNdptTimesPt(mass,beta,T,q,norm,ymax,name);
 128     break;
 129   case kOneOverPtdNdpt:
 130     return GetTsallisBWdNdpt(mass,beta,T,q,norm,ymax,name);
 131     break;
 132   case kOneOverMtdNdmt:
 133     AliFatal("Not implemented");
 134     break;
 135   default:
 136     AliFatal("Not implemented");
 137   }
 138
 139   return 0;
 140
 141 }
 142
 143
 144 TF1 * AliPWGFunc::GetMTExp(Double_t mass, Double_t T, Double_t norm, const char * name){
 145
 146   // Simple exponential in 1/mt*MT
 147   switch (fVarType) {
 148   case kdNdpt:
 149     return GetMTExpdNdptTimesPt(mass,T,norm,name);
 150     break;
 151   case kOneOverPtdNdpt:
 152     return GetMTExpdNdpt(mass,T,norm,name);
 153     break;
 154   case kOneOverMtdNdmt:
 155     return GetMTExpdNdmt(mass,T,norm,name,kOneOverMtdNdmt);
 156     break;
 157   case kdNdmt:
 158     return GetMTExpdNdmt(mass,T,norm,name,kdNdmt);
 159     break;
 160   case kOneOverMtdNdmtMinusM:
 161     return GetMTExpdNdmt(mass,T,norm,name,kOneOverMtdNdmtMinusM);
 162     break;
 163   default:
 164     AliFatal("Not implemented");
 165   }
 166
 167   return 0;
 168
 169
 170 }
 171
 172
 173 TF1 * AliPWGFunc::GetBoseEinstein(Double_t mass, Double_t T, Double_t norm, const char * name){
 174
 175   // Bose einstein
 176   switch (fVarType) {
 177   case kdNdpt:
 178     return GetBoseEinsteindNdptTimesPt(mass,T,norm,name);
 179     break;
 180   case kOneOverPtdNdpt:
 181     return GetBoseEinsteindNdpt(mass,T,norm,name);
 182     break;
 183   case kOneOverMtdNdmt:
 184     AliFatal("Not implemented");
 185     break;
 186   default:
 187     AliFatal("Not implemented");
 188   }
 189
 190   return 0;
 191
 192
 193 }
 194
 195 TF1 * AliPWGFunc::GetFermiDirac(Double_t mass, Double_t T, Double_t norm, const char * name){
 196
 197   // Simple exponential in 1/mt*MT
 198   switch (fVarType) {
 199   case kdNdpt:
 200     return GetFermiDiracdNdptTimesPt(mass,T,norm,name);
 201     break;
 202   case kOneOverPtdNdpt:
 203     return GetFermiDiracdNdpt(mass,T,norm,name);
 204     break;
 205   case kOneOverMtdNdmt:
 206     AliFatal("Not implemented");
 207     break;
 208   default:
 209     AliFatal("Not implemented");
 210   }
 211
 212   return 0;
 213
 214
 215 }
 216
 217
 218 TF1 * AliPWGFunc::GetPTExp(Double_t T, Double_t norm, const char * name){
 219
 220   // Simple exponential in 1/mt*MT
 221   switch (fVarType) {
 222   case kdNdpt:
 223     return GetPTExpdNdptTimesPt(T,norm,name);
 224     break;
 225   case kOneOverPtdNdpt:
 226     AliFatal("Not implemented");
 227     break;
 228   case kOneOverMtdNdmt:
 229     AliFatal("Not implemented");
 230     break;
 231   default:
 232     AliFatal("Not implemented");
 233   }
 234
 235   return 0;
 236
 237
 238 }
 239
 240
 241 TF1 * AliPWGFunc::GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
 242   // Levi function (aka Tsallis)
 243   switch (fVarType) {
 244   case kdNdpt:
 245     return GetLevidNdptTimesPt(mass,T,n,norm,name);
 246     break;
 247   case kOneOverPtdNdpt:
 248     return GetLevidNdpt(mass,T,n,norm,name);
 249     break;
 250   case kOneOverMtdNdmt:
 251     return GetLevidNdmt(mass,T,n,norm,name,kOneOverMtdNdmt);
 252     break;
 253   case kdNdmt:
 254     return GetLevidNdmt(mass,T,n,norm,name,kdNdmt);
 255     break;
 256   case kOneOverMtdNdmtMinusM:
 257     return GetLevidNdmt(mass,T,n,norm,name,kOneOverMtdNdmtMinusM);
 258     break;
 259   default:
 260     AliFatal("Not implemented");
 261   }
 262
 263   return 0;
 264
 265
 266 }
 267
 268 TF1 * AliPWGFunc::GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name){
 269   // power law Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
 270   // This is sometimes also called Hagedorn or modified Hagedorn
 271
 272   switch (fVarType) {
 273   case kdNdpt:
 274     return GetPowerLawdNdptTimesPt(pt0,n,norm,name);
 275     break;
 276   case kOneOverPtdNdpt:
 277     return GetPowerLawdNdpt(pt0,n,norm,name);
 278     break;
 279   case kOneOverMtdNdmt:
 280     AliFatal("Not Implemented");
 281     //    return GetUA1dNdmt(mass,T,n,norm,name);
 282     break;
 283   default:
 284     AliFatal("Not implemented");
 285   }
 286
 287   return 0;
 288
 289
 290 }
 291
 292 TF1 * AliPWGFunc::GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char * name) {
 293   // UA1 parametrization Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
 294
 295   switch (fVarType) {
 296   case kdNdpt:
 297
 298     fLastFunc = new TF1 (name, StaticUA1Func, 0.0, 10, 6);
 299     fLastFunc->FixParameter(0,mass);
 300     fLastFunc->SetParameter(1,p0star);
 301     fLastFunc->SetParameter(2,pt0);
 302     fLastFunc->SetParameter(3,n);
 303     fLastFunc->SetParameter(4,T);
 304     fLastFunc->SetParameter(5,norm);
 305     fLastFunc->SetParLimits(1,0.01,1);
 306     fLastFunc->SetParLimits(2,0.01,100);
 307     fLastFunc->SetParLimits(3,0.01,100);
 308     fLastFunc->SetParLimits(4,0.01,100);
 309     fLastFunc->SetParNames("mass","p0star","pt0","n","T","norm");
 310     fLastFunc->SetNpx(5000);
 311     fLastFunc->SetLineWidth(fLineWidth);
 312     return fLastFunc;
 313
 314     break;
 315   case kOneOverPtdNdpt:
 316     fLastFunc = new TF1 (name, StaticUA1FuncOneOverPt, 0.0, 10, 6);
 317     fLastFunc->FixParameter(0,mass);
 318     fLastFunc->SetParameter(1,p0star);
 319     fLastFunc->SetParameter(2,pt0);
 320     fLastFunc->SetParameter(3,n);
 321     fLastFunc->SetParameter(4,T);
 322     fLastFunc->SetParameter(5,norm);
 323     fLastFunc->SetParLimits(1,0.01,1);
 324     fLastFunc->SetParLimits(2,0.01,100);
 325     fLastFunc->SetParLimits(3,0.01,100);
 326     fLastFunc->SetParLimits(4,0.01,100);
 327     fLastFunc->SetParNames("mass","p0star","pt0","n","T","norm");
 328     fLastFunc->SetNpx(5000);
 329     fLastFunc->SetLineWidth(fLineWidth);
 330     return fLastFunc;
 331
 332     break;
 333   case kOneOverMtdNdmt:
 334     AliFatal("Not Implemented");
 335     //    return GetUA1dNdmt(mass,T,n,norm,name);
 336     break;
 337   default:
 338     AliFatal("Not implemented");
 339   }
 340
 341   return 0;
 342 }
 343
 344
 345
 346
 347 // ________________________________________________________________________
 348
 349 // Backend (private functions and support functions for numerical integration)
 350
 351 Double_t AliPWGFunc::StaticHistoFunc(const double * x, const double* p){
 352
 353   // provides a function interpolating a histo with a spline;
 354   // using double to store a pointer... This is a bad hack. To be replaced
 355
 356   double norm = p[0];
 357
 358   TObject * h     = (TObject*) Long64_t(p[1]);
 359
 360 //    Int_t bin = h->FindBin(x[0]);
 361 //    double value = h->GetBinContent(bin);
 362
 363
 364   // static TH1 * oldptr = 0;
 365   // static TSpline3 * spl = 0;
 366   // if (h!=oldptr) {
 367   // FIXME: recheck static pointers
 368   TSpline3 * spl  = 0;
 369   if(h->InheritsFrom("TH1")) {
 370     if ( ((TH1*)h)->FindBin(x[0]) > ((TH1*)h)->GetNbinsX()) return 0;
 371     spl= new TSpline3((TH1*)h);
 372   }
 373   else if(h->InheritsFrom("TGraph")) spl= new TSpline3("fGraph",(TGraph*)h);
 374   else {
 375     Printf("AliPWGFunc::StaticHistoFunc: Unsupported type");
 376     return 0;
 377   }
 378     //  }
 379   double value = spl->Eval(x[0]);
 380   delete spl;
 381
 382   return value*norm;
 383
 384 }
 385
 386 Double_t AliPWGFunc::StaticUA1Func(const double * x, const double* p) {
 387
 388
 389   // "mass","p0star","pt0","n","T","norm"
 390   Double_t mass   = p[0];
 391   Double_t p0star = p[1];
 392   Double_t pt0    = p[2];
 393   Double_t n      = p[3];
 394   Double_t temp   = p[4];
 395   Double_t norm   = p[5];
 396
 397   Double_t xx = x[0];
 398
 399   static AliPWGFunc * self = new AliPWGFunc;
 400   static TF1 * fPLaw   = self->GetPowerLawdNdptTimesPt(pt0, n, norm, "fLocalPLawUA1");
 401   static TF1 * fPMTExp = self->GetMTExpdNdptTimesPt   (mass, temp, norm, "fLocalMTexpUA1");
 402
 403   fPLaw->SetParameters(norm,pt0,n);
 404   fPMTExp->SetParameters(1,temp);
 405
 406
 407   Double_t normMT =fPMTExp->Eval(p0star) >0 ? fPLaw->Eval(p0star) / fPMTExp->Eval(p0star) *  fPMTExp->GetParameter(0) : 1;
 408   fPMTExp->SetParameter(0,normMT);
 409
 410
 411   if (TMath::Abs(fPMTExp->Eval(p0star) - fPLaw->Eval(p0star)) > 0.0001 ) {
 412     Printf("AliPWGFunc::StaticUA1Func - Wrong norm") ;
 413     Printf(" p0* %f  NMT: %f  N: %f  PL: %f  MT: %f", p0star, normMT, norm, fPLaw->Eval(p0star), fPMTExp->Eval(p0star));
 414   }
 415
 416   if (xx > p0star)  return fPLaw->Eval(xx);
 417   return fPMTExp->Eval(xx);
 418
 419
 420
 421 }
 422 Double_t AliPWGFunc::StaticUA1FuncOneOverPt(const double * x, const double* p) {
 423   // UA1 func
 424
 425   // "mass","p0star","pt0","n","T","norm"
 426   Double_t mass   = p[0];
 427   Double_t p0star = p[1];
 428   Double_t pt0    = p[2];
 429   Double_t n      = p[3];
 430   Double_t temp   = p[4];
 431   Double_t norm   = p[5];
 432
 433   Double_t xx = x[0];
 434
 435   static AliPWGFunc * self = new AliPWGFunc;
 436   static TF1 * fPLaw   = self->GetPowerLawdNdpt(pt0, n, norm, "fLocalPLawUA1");
 437   static TF1 * fPMTExp = self->GetMTExpdNdpt   (mass, temp, norm, "fLocalMTexpUA1");
 438
 439   fPLaw->SetParameters(norm,pt0,n);
 440   fPMTExp->SetParameters(1,temp);
 441
 442
 443   Double_t normMT =fPMTExp->Eval(p0star) >0 ? fPLaw->Eval(p0star) / fPMTExp->Eval(p0star) *  fPMTExp->GetParameter(0) : 1;
 444   fPMTExp->SetParameter(0,normMT);
 445
 446
 447   if (TMath::Abs(fPMTExp->Eval(p0star) - fPLaw->Eval(p0star)) > 0.0001 ) {
 448     Printf("AliPWGFunc::StaticUA1Func - Wrong norm") ;
 449     Printf(" p0* %f  NMT: %f  N: %f  PL: %f  MT: %f", p0star, normMT, norm, fPLaw->Eval(p0star), fPMTExp->Eval(p0star));
 450   }
 451
 452   if (xx > p0star)  return fPLaw->Eval(xx);
 453   return fPMTExp->Eval(xx);
 454
 455
 456
 457 }
 458
 459
 460 Double_t AliPWGFunc::IntegrandBG(const double * x, const double* p){
 461   // integrand for boltzman-gibbs blast wave
 462      // x[0] -> r (radius)
 463      // p[0] -> mass
 464      // p[1] -> pT (transverse momentum)
 465      // p[2] -> beta_max (surface velocity)
 466      // p[3] -> T (freezout temperature)
 467      // p[4] -> n (velocity profile)
 468
 469
 470   double x0 = x[0];
 471
 472   double mass     = p[0];
 473   double pT       = p[1];
 474   double beta_max = p[2];
 475   double temp     = p[3];
 476   Double_t n      = p[4];
 477
 478   // Keep beta within reasonable limits
 479   Double_t beta = beta_max * TMath::Power(x0, n);
 480   if (beta > 0.9999999999999999) beta = 0.9999999999999999;
 481
 482   double mT      = TMath::Sqrt(mass*mass+pT*pT);
 483
 484   double rho0   = TMath::ATanH(beta);
 485   double arg00 = pT*TMath::SinH(rho0)/temp;
 486   if (arg00 > 700.) arg00 = 700.; // avoid FPE
 487   double arg01 = mT*TMath::CosH(rho0)/temp;
 488   double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01);
 489
 490   //  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);
 491
 492   return f0;
 493 }
 494
 495
 496
 497 Double_t AliPWGFunc::StaticBGdNdPt(const double * x, const double* p) {
 498
 499   // implementation of BGBW (1/pt dNdpt)
 500
 501   double pT = x[0];;
 502
 503
 504   double mass    = p[0];
 505   double beta    = p[1];
 506   double temp    = p[2];
 507   double n       = p[3];
 508   double norm    = p[4];
 509
 510   static TF1 * fIntBG = 0;
 511   if(!fIntBG)
 512     fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 5);
 513
 514   fIntBG->SetParameters(mass, pT, beta, temp,n);
 515   double result = fIntBG->Integral(0,1);
 516   //  printf ("[%4.4f], Int :%f\n", pT, result);
 517   return result*norm;//*1e30;;
 518
 519 }
 520
 521 Double_t AliPWGFunc::StaticBGdNdPtTimesPt(const double * x, const double* p) {
 522   // BGBW dNdpt implementation
 523   return x[0]*StaticBGdNdPt(x,p);
 524 }
 525
 526 Double_t AliPWGFunc::StaticBGdNdMtTimesMt(const double * x, const double* p) {
 527   // BGBW dNdpt implementation
 528   // X0 is mt here
 529   Double_t pt = TMath::Sqrt(x[0]*x[0]-p[0]*p[0]);
 530   return pt*StaticBGdNdPt(&pt,p);
 531 }
 532
 533
 534 TF1 * AliPWGFunc::GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t temp,
 535                               Double_t n, Double_t norm, const char * name){
 536
 537   // BGBW 1/pt dNdpt
 538
 539   fLastFunc = new TF1 (name, StaticBGdNdPt, 0.0, 10, 5);
 540   fLastFunc->SetParameters(mass,beta,temp,n,norm);
 541   fLastFunc->FixParameter(0,mass);
 542   fLastFunc->SetParNames("mass", "#beta", "T", "n", "norm");
 543   fLastFunc->SetLineWidth(fLineWidth);
 544   return fLastFunc;
 545
 546 }
 547
 548
 549 //_____________________________________________________________________
 550 // Tsallis
 551
 552 Double_t AliPWGFunc::IntegrandTsallis(const double * x, const double* p){
 553
 554   // integrand for numerical integration (tsallis)
 555
 556   Double_t r   = x[0];
 557   Double_t phi = x[1];
 558   Double_t y   = x[2];
 559
 560   Double_t mass = p[0];
 561   Double_t pt   = p[1];
 562   Double_t beta = p[2];
 563   Double_t temp    = p[3];
 564   Double_t q    = p[4];
 565
 566   Double_t mt      = TMath::Sqrt(mass*mass+pt*pt);
 567
 568   Double_t rho    = TMath::ATanH(beta*r); // TODO: implement different velocity profiles
 569
 570   Double_t res = mt*
 571     r*TMath::CosH(y) *TMath::Power( (
 572                                      1+(q-1)/temp * (
 573                                                   mt*TMath::CosH(y)*TMath::CosH(rho) -
 574                                                   pt*TMath::SinH(rho)*TMath::Cos(phi)
 575                                                   )
 576                                      ),
 577                                        -1/(q-1)
 578                                     );
 579
 580
 581   return res;
 582 }
 583
 584
 585
 586 Double_t AliPWGFunc::StaticTsallisdNdPt(const double * x, const double* p) {
 587
 588   // tsallis BW implementation 1/pt dNdpt
 589
 590   double pT = x[0];;
 591
 592
 593   double mass = p[0];
 594   double beta = p[1];
 595   double temp    = p[2];
 596   double q    = p[3];
 597
 598   Double_t ymax = p[5];
 599
 600
 601   static TF3 * fInt = 0;
 602   if(!fInt){
 603     fInt = new TF3 ("fIntTsa", IntegrandTsallis, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
 604 //     fInt->SetNpx(10000);
 605 //     fInt->SetNpy(10000);
 606 //     fInt->SetNpz(10000);
 607   }
 608
 609   fInt->SetParameters(mass, pT, beta, temp, q);
 610   double result = fInt->Integral(0,1, -TMath::Pi(), TMath::Pi(), -ymax, ymax);
 611   //  double result = fInt->Integral(0,1, -2, 2, -ymax, ymax);
 612
 613   return result*p[4];//*1e30;;
 614
 615 }
 616
 617 Double_t AliPWGFunc::StaticTsallisdNdPtTimesPt(const double * x, const double* p) {
 618
 619   // tsallis BW , implementatio of dNdpt
 620   return x[0]*StaticTsallisdNdPt(x,p);
 621
 622 }
 623
 624 TF1 * AliPWGFunc::GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t temp, Double_t q,
 625                                    Double_t norm, Double_t ymax,const char * name){
 626
 627
 628   // tsallis BW, 1/pt dNdpt
 629
 630   fLastFunc = new TF1 (name, StaticTsallisdNdPt, 0.0, 10, 6);
 631   fLastFunc->SetParameters(mass,beta,temp,q,norm,ymax);
 632   fLastFunc->SetParLimits(1,0.0,0.99);
 633   fLastFunc->SetParLimits(2,0.01,0.99);
 634   fLastFunc->SetParLimits(3,1.0001,1.9);
 635   fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax");
 636   fLastFunc->SetLineWidth(fLineWidth);
 637   return fLastFunc;
 638
 639 }
 640
 641 // Times Pt funcs
 642 // Boltzmann-Gibbs Blast Wave
 643 TF1 * AliPWGFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t n,
 644                                      Double_t norm, const char * name){
 645
 646   // BGBW, dNdpt
 647
 648   fLastFunc = new TF1 (name, StaticBGdNdPtTimesPt, 0.0, 10, 5);
 649   fLastFunc->SetParameters(mass,beta,temp,n,norm);
 650   fLastFunc->FixParameter(0,mass);
 651   fLastFunc->SetParNames("mass", "#beta", "temp", "n", "norm");
 652   fLastFunc->SetLineWidth(fLineWidth);
 653   return fLastFunc;
 654
 655
 656 }
 657
 658 // Boltzmann-Gibbs Blast Wave
 659 TF1 * AliPWGFunc::GetBGBWdNdptTimesMt(Double_t mass, Double_t beta, Double_t temp, Double_t n,
 660                                       Double_t norm, const char * name){
 661
 662   // BGBW, dNdpt
 663   // 1/Mt dN/dmt
 664   fLastFunc = new TF1 (name, StaticBGdNdMtTimesMt, 0.0, 10, 5);
 665   fLastFunc->SetParameters(mass,beta,temp,n,norm);
 666   fLastFunc->FixParameter(0,mass);
 667   fLastFunc->SetParNames("mass", "#beta", "temp", "n", "norm");
 668   fLastFunc->SetLineWidth(fLineWidth);
 669   return fLastFunc;
 670
 671
 672 }
 673
 674 TF1 * AliPWGFunc::GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t q,
 675                                           Double_t norm, Double_t ymax, const char * name){
 676
 677 // Tsallis blast wave, dNdpt
 678
 679   fLastFunc = new TF1 (name, StaticTsallisdNdPtTimesPt, 0.0, 10, 6);
 680   fLastFunc->SetParameters(mass,beta,temp,q,norm,ymax);
 681   fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax");
 682   fLastFunc->SetLineWidth(fLineWidth);
 683   return fLastFunc;
 684
 685
 686 }
 687
 688
 689
 690 TF1 * AliPWGFunc::GetMTExpdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
 691
 692   // Simple exponential in 1/mt*MT, as a function of dNdpt
 693   char formula[500];
 694   snprintf(formula,500,"[0]*x*exp(-sqrt(x**2+%f**2)/[1])", mass);
 695   fLastFunc=new TF1(name,formula,0,10);
 696   fLastFunc->SetParameters(norm, temp);
 697   fLastFunc->SetParLimits(1, 0.01, 10);
 698   fLastFunc->SetParNames("norm", "T");
 699   fLastFunc->SetLineWidth(fLineWidth);
 700   return fLastFunc;
 701
 702
 703 }
 704
 705 TF1 * AliPWGFunc::GetBoseEinsteindNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
 706
 707   // Bose einstein distribution as a function of dNdpt
 708   char formula[500];
 709   snprintf(formula,500,"[0]*x*1./(exp(sqrt(x**2+%f**2)/[1])-1)", mass);
 710   fLastFunc=new TF1(name,formula,0,10);
 711   fLastFunc->SetParameters(norm, temp);
 712   fLastFunc->SetParLimits(1, 0.01, 10);
 713   fLastFunc->SetParNames("norm", "T");
 714   fLastFunc->SetLineWidth(fLineWidth);
 715   return fLastFunc;
 716
 717
 718 }
 719
 720 TF1 * AliPWGFunc::GetFermiDiracdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
 721
 722   // Bose einstein distribution as a function of dNdpt
 723   char formula[500];
 724   snprintf(formula,500,"[0]*x*1./(exp(sqrt(x**2+%f**2)/[1])+1)", mass);
 725   fLastFunc=new TF1(name,formula,0,10);
 726   fLastFunc->SetParameters(norm, temp);
 727   fLastFunc->SetParLimits(1, 0.01, 10);
 728   fLastFunc->SetParNames("norm", "T");
 729   fLastFunc->SetLineWidth(fLineWidth);
 730   return fLastFunc;
 731
 732
 733 }
 734
 735
 736
 737 TF1 * AliPWGFunc::GetPTExpdNdptTimesPt(Double_t temp, Double_t norm, const char * name){
 738
 739   // Simple exponential in 1/pt*dNdpT, as a function of dNdpt
 740   char formula[500];
 741   snprintf(formula,500,"[0]*x*exp(-x/[1])");
 742   fLastFunc=new TF1(name,formula,0,10);
 743   fLastFunc->SetParameters(norm, temp);
 744   fLastFunc->SetParLimits(1, 0.01, 10);
 745   fLastFunc->SetParNames("norm", "T");
 746   fLastFunc->SetLineWidth(fLineWidth);
 747   return fLastFunc;
 748
 749
 750 }
 751
 752
 753 TF1 * AliPWGFunc::GetBoltzmanndNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
 754   // Boltzmann (exp in 1/mt*dNdmT times mt) as a function of dNdpt
 755  char formula[500];
 756  snprintf(formula,500,"[0]*x*sqrt(x**2+%f**2)*exp(-sqrt(x**2+%f**2)/[1])", mass,mass);
 757   fLastFunc=new TF1(name,formula,0,10);
 758   fLastFunc->SetParameters(norm, temp);
 759   fLastFunc->SetParLimits(1, 0.01, 10);
 760   fLastFunc->SetParNames("norm", "T");
 761   fLastFunc->SetLineWidth(fLineWidth);
 762   return fLastFunc;
 763
 764
 765 }
 766
 767
 768 // Tsallis (no BW, a la CMS)
 769 // TF1 * AliPWGFunc::GetTsallisdNdptTimesPt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name){
 770
 771 //   char formula[500];
 772 //   //  sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2)-%f)),(-1/([2]-1)))", mass, mass); //CMS
 773 //   sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass);  // STAR
 774 //   //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
 775 //   fLastFunc=new TF1(name,formula,0,10);
 776 //   fLastFunc->SetParameters(norm, T, q);
 777 //   fLastFunc->SetParLimits(1, 0.001, 10);
 778 //   fLastFunc->SetParNames("norm", "T", "q");
 779 //   fLastFunc->SetLineWidth(fLineWidth);
 780 //   return fLastFunc;
 781
 782
 783 // }
 784
 785
 786 TF1 * AliPWGFunc::GetLevidNdptTimesPt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
 787
 788   // Levi function, dNdpt
 789   char formula[500];
 790
 791   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])");
 792   //  sprintf(formula,"( x*[0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x))/([1]*[2])  )^(-[1])");
 793   fLastFunc=new TF1(name,formula,0,10);
 794   fLastFunc->SetParameters(norm, n, temp,mass);
 795   fLastFunc->SetParLimits(2, 0.01, 10);
 796   fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
 797   fLastFunc->FixParameter(3,mass);
 798   fLastFunc->SetLineWidth(fLineWidth);
 799   return fLastFunc;
 800
 801
 802 }
 803
 804 TF1 * AliPWGFunc::GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name){
 805
 806   // PowerLaw function, dNdpt
 807   char formula[500];
 808
 809   snprintf(formula,500,"x*[0]*( 1 + x/[1] )^(-[2])");
 810   fLastFunc=new TF1(name,formula,0,10);
 811   fLastFunc->SetParameters(norm, pt0, n);
 812   fLastFunc->SetParLimits(1, 0.01, 10);
 813   //fLastFunc->SetParLimits(2, 0.01, 50);
 814   fLastFunc->SetParNames("norm", "pt0", "n");
 815   fLastFunc->SetLineWidth(fLineWidth);
 816   return fLastFunc;
 817
 818
 819 }
 820
 821 TF1 * AliPWGFunc::GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name){
 822
 823   // PowerLaw function, 1/pt dNdpt
 824   char formula[500];
 825
 826   snprintf(formula,500," [0]*( 1 + x/[1] )^(-[2])");
 827   fLastFunc=new TF1(name,formula,0,10);
 828   fLastFunc->SetParameters(norm, pt0, n);
 829   //  fLastFunc->SetParLimits(2, 0.01, 10);
 830   fLastFunc->SetParNames("norm", "pt0", "n");
 831   fLastFunc->SetLineWidth(fLineWidth);
 832   return fLastFunc;
 833
 834
 835 }
 836
 837
 838 TF1 * AliPWGFunc::GetLevidNdpt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
 839
 840   // Levi function, dNdpt
 841   char formula[500];
 842
 843   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])");
 844   fLastFunc=new TF1(name,formula,0,10);
 845   fLastFunc->SetParameters(norm, n, temp,mass);
 846   fLastFunc->SetParLimits(2, 0.01, 10);
 847   fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
 848   fLastFunc->FixParameter(3,mass);
 849   fLastFunc->SetLineWidth(fLineWidth);
 850   return fLastFunc;
 851
 852
 853 }
 854
 855 TF1 * AliPWGFunc::GetLevidNdmt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name, VarType_t var){
 856
 857   // Levi function, 1/mt dNdmt
 858   char formula[500];
 859   if (var == kOneOverMtdNdmt)
 860     snprintf(formula,500,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (x -[3])/([1]*[2])  )^(-[1])");
 861   else if (var == kdNdmt)
 862     snprintf(formula,500,"( x*[0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (x-[3])/([1]*[2])  )^(-[1])");
 863   if (var == kOneOverMtdNdmtMinusM)
 864     snprintf(formula,500,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (x)/([1]*[2])  )^(-[1])");
 865
 866   //sprintf(formula,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + x/([1]*[2])  )^(-[1])");
 867   //  sprintf(formula,"[0] * ( 1 + x/([1]*[2])  )^(-[1])");
 868   fLastFunc=new TF1(name,formula,0,10);
 869   fLastFunc->SetParameters(norm, n, temp,mass);
 870   fLastFunc->SetParLimits(2, 0.01, 10);
 871   fLastFunc->SetParNames("norm", "n", "T", "mass");
 872   fLastFunc->FixParameter(3,mass);
 873   fLastFunc->SetLineWidth(fLineWidth);
 874   return fLastFunc;
 875
 876
 877 }
 878
 879
 880
 881
 882 // Test Function
 883 Double_t AliPWGFunc::IntegrandTest(const double * x, const double* p){
 884
 885   // test function
 886
 887   Double_t y = x[0];
 888
 889   Double_t mass = p[0];
 890   Double_t pt   = p[1];
 891   Double_t temp    = p[2];
 892
 893   Double_t mt      = TMath::Sqrt(mass*mass+pt*pt);
 894
 895   return mt*TMath::CosH(y)*TMath::Exp(-mt*TMath::CosH(y)/temp);
 896
 897 }
 898
 899 Double_t AliPWGFunc::StaticTest(const double * x, const double* p) {
 900
 901   // test function
 902
 903   double pT = x[0];;
 904
 905
 906   double mass = p[0];
 907   double temp    = p[1];
 908   Double_t ymax = p[3];
 909
 910
 911   static TF3 * fIntTest = 0;
 912   if(!fIntTest){
 913     fIntTest = new TF3 ("fIntTest", IntegrandTest, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
 914     //    fInt->SetNpx(10000);
 915   }
 916
 917   fIntTest->SetParameters(mass, pT, temp);
 918   double result = fIntTest->Integral(-ymax, ymax);
 919
 920   return result*p[2];//*1e30;;
 921
 922 }
 923
 924 TF1 * AliPWGFunc::GetTestFunc(Double_t mass, Double_t temp, Double_t norm, Double_t ymax, const char * name){
 925
 926   // test function
 927
 928   fLastFunc = new TF1 (name, StaticTest, 0.0, 10, 4);
 929   fLastFunc->SetParameters(mass,temp,norm,ymax);
 930   fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
 931   fLastFunc->SetLineWidth(fLineWidth);
 932   return fLastFunc;
 933
 934 }
 935
 936
 937 //___________________________________________________________
 938
 939
 940 TF1 * AliPWGFunc::GetMTExpdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
 941   // Simple exp in 1/mt dNdmt, as a function of dNdpt
 942   // mt scaling
 943   char formula[500];
 944   snprintf(formula,500,"[0]*exp(-sqrt(x**2+%f**2)/[1])", mass);
 945   fLastFunc=new TF1(name,formula,0,10);
 946   fLastFunc->SetParameters(norm, temp);
 947   fLastFunc->SetParLimits(1, 0.01, 10);
 948   fLastFunc->SetParNames("norm", "T");
 949   fLastFunc->SetLineWidth(fLineWidth);
 950   return fLastFunc;
 951 }
 952
 953
 954 TF1 * AliPWGFunc::GetMTExpdNdmt(Double_t mass, Double_t temp, Double_t norm, const char * name, VarType_t var){
 955
 956   // Levi function, 1/mt dNdmt1/mt dNdmt,
 957   char formula[500];
 958   if (var == kOneOverMtdNdmt)
 959     snprintf(formula,500,"[0] * exp (-x/[1]) + %f ", mass);
 960   else if (var == kdNdmt)
 961     snprintf(formula,500,"[0] * x * exp (-x/[1]) + %f ", mass);
 962   if (var == kOneOverMtdNdmtMinusM)
 963     snprintf(formula,500,"[0] * exp (-x/[1])");
 964
 965   //sprintf(formula,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + x/([1]*[2])  )^(-[1])");
 966   //  sprintf(formula,"[0] * ( 1 + x/([1]*[2])  )^(-[1])");
 967   fLastFunc=new TF1(name,formula,0,10);
 968   fLastFunc->SetParameters(norm, temp);
 969   fLastFunc->SetParLimits(1, 0.01, 10);
 970   fLastFunc->SetParNames("norm", "T");
 971   fLastFunc->SetLineWidth(fLineWidth);
 972   return fLastFunc;
 973
 974
 975 }
 976
 977
 978
 979 TF1 * AliPWGFunc::GetBoseEinsteindNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
 980   // bose einstein
 981   char formula[500];
 982   snprintf(formula,500,"[0]*1./(exp(sqrt(x**2+%f**2)/[1])-1)", mass);
 983   fLastFunc=new TF1(name,formula,0,10);
 984   fLastFunc->SetParameters(norm, temp);
 985   fLastFunc->SetParLimits(1, 0.01, 10);
 986   fLastFunc->SetParNames("norm", "T");
 987   fLastFunc->SetLineWidth(fLineWidth);
 988   return fLastFunc;
 989 }
 990
 991 TF1 * AliPWGFunc::GetFermiDiracdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
 992   // bose einstein
 993   char formula[500];
 994   snprintf(formula,500,"[0]*1./(exp(sqrt(x**2+%f**2)/[1])+1)", mass);
 995   fLastFunc=new TF1(name,formula,0,10);
 996   fLastFunc->SetParameters(norm, temp);
 997   fLastFunc->SetParLimits(1, 0.01, 10);
 998   fLastFunc->SetParNames("norm", "T");
 999   fLastFunc->SetLineWidth(fLineWidth);
1000   return fLastFunc;
1001 }
1002
1003
1004 // // Simple tsallis (a la CMS)
1005 // TF1 * AliPWGFunc::GetTsallisdNdpt(Double_t mass, Double_t temp, Double_t q, Double_t norm, const char * name){
1006
1007 //   char formula[500];
1008 //   sprintf(formula,"[0]*sqrt(x**2+%f**2)*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass,mass);
1009 //   fLastFunc=new TF1(name,formula,0,10);
1010 //   fLastFunc->SetParameters(norm, temp, q);
1011 //   fLastFunc->SetParLimits(1, 0.01, 10);
1012 //   fLastFunc->SetParNames("norm", "T", "q");
1013 //   fLastFunc->SetLineWidth(fLineWidth);
1014 //   return fLastFunc;
1015 // }