PWG2/SPECTRA/Fit/AliBWFunc.cxx

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