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