AliPhysics  a5cd6b6 (a5cd6b6)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliPWGFunc.cxx
Go to the documentation of this file.
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 
30 
31 AliPWGFunc::AliPWGFunc () : fLastFunc(0), fLineWidth(1), fVarType(kdNdpt) {
32 
33  // ctor
34  fLineWidth = 1;
35 }
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 
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 
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;
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;
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 
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 
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
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
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 
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 
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 
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 
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 // }
TF1 * GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t temp, Double_t n, Double_t norm, const char *name="fBGBW")
Definition: AliPWGFunc.cxx:534
static Double_t IntegrandTsallis(const double *x, const double *p)
Definition: AliPWGFunc.cxx:552
TF1 * GetMTExp(Double_t mass, Double_t T, Double_t norm, const char *name="fMtExp")
Definition: AliPWGFunc.cxx:144
TF1 * GetBoseEinsteindNdpt(Double_t mass, Double_t T, Double_t norm, const char *name="fBoseEinstein")
Definition: AliPWGFunc.cxx:979
const Double_t ymax
Definition: AddTaskCFDStar.C:7
double Double_t
Definition: External.C:58
static Double_t StaticUA1FuncOneOverPt(const double *x, const double *p)
Definition: AliPWGFunc.cxx:422
TF1 * GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char *name="fLevi")
Definition: AliPWGFunc.cxx:804
TF1 * GetMTExpdNdmt(Double_t mass, Double_t temp, Double_t norm, const char *name, VarType_t var)
Definition: AliPWGFunc.cxx:954
TF1 * GetBoseEinstein(Double_t mass, Double_t T, Double_t norm, const char *name="fBoseEinstein")
Definition: AliPWGFunc.cxx:173
long long Long64_t
Definition: External.C:43
TF1 * GetPTExpdNdptTimesPt(Double_t T, Double_t norm, const char *name="fPtExpTimesPt")
Definition: AliPWGFunc.cxx:737
static Double_t StaticUA1Func(const double *x, const double *p)
Definition: AliPWGFunc.cxx:386
Double_t mass
VarType_t fVarType
Definition: AliPWGFunc.h:201
TF1 * GetTsallisBW(Double_t mass, Double_t beta, Double_t T, Double_t q, Double_t norm, Double_t ymax=0.5, const char *name="fTsallisBW")
Definition: AliPWGFunc.cxx:121
TF1 * GetFermiDiracdNdpt(Double_t mass, Double_t T, Double_t norm, const char *name="fFermiDirac")
Definition: AliPWGFunc.cxx:991
TF1 * GetBoseEinsteindNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char *name="fBoseEinstein")
Definition: AliPWGFunc.cxx:705
TF1 * GetGraphFunc(TGraph *h, const char *name="fHisto")
Definition: AliPWGFunc.cxx:60
TF1 * GetLevidNdptTimesPt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char *name="fLevi")
Definition: AliPWGFunc.cxx:786
static Double_t StaticTsallisdNdPtTimesPt(const double *x, const double *p)
Definition: AliPWGFunc.cxx:617
TF1 * GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T, Double_t q, Double_t norm, Double_t ymax=0.5, const char *name="fTsallisBWTimesPt")
Definition: AliPWGFunc.cxx:674
TF1 * GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T, Double_t n, Double_t norm, const char *name="fBGBWTimesPt")
Definition: AliPWGFunc.cxx:643
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
static Double_t StaticBGdNdMtTimesMt(const double *x, const double *p)
Definition: AliPWGFunc.cxx:526
TF1 * GetHistoFunc(TH1 *h, const char *name="fHisto")
Definition: AliPWGFunc.cxx:44
static Double_t IntegrandTest(const double *x, const double *p)
Definition: AliPWGFunc.cxx:883
TF1 * GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char *name="fPowerLaw")
Definition: AliPWGFunc.cxx:268
TF1 * GetBGBW(Double_t mass, Double_t beta, Double_t T, Double_t n, Double_t norm, const char *name="fBGBW")
Definition: AliPWGFunc.cxx:77
TF1 * GetLevidNdmt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char *name="fLeviMt", VarType_t var=kOneOverMtdNdmt)
Definition: AliPWGFunc.cxx:855
static Double_t StaticBGdNdPt(const double *x, const double *p)
Definition: AliPWGFunc.cxx:497
Width_t fLineWidth
Definition: AliPWGFunc.h:200
TF1 * GetBGBWdNdptTimesMt(Double_t mass, Double_t beta, Double_t T, Double_t n, Double_t norm, const char *name="fBGBWTimesMt")
Definition: AliPWGFunc.cxx:659
TF1 * GetMTExpdNdpt(Double_t mass, Double_t T, Double_t norm, const char *name="fExp")
Definition: AliPWGFunc.cxx:940
TF1 * GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char *name="fUA1")
Definition: AliPWGFunc.cxx:292
TF1 * GetFermiDiracdNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char *name="fFermiDirac")
Definition: AliPWGFunc.cxx:720
TF1 * GetPTExp(Double_t T, Double_t norm, const char *name="fPtExp")
Definition: AliPWGFunc.cxx:218
TF1 * GetMTExpdNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char *name="fMtExpTimesPt")
Definition: AliPWGFunc.cxx:690
static Double_t StaticTsallisdNdPt(const double *x, const double *p)
Definition: AliPWGFunc.cxx:586
TF1 * fLastFunc
Definition: AliPWGFunc.h:199
TF1 * GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t T, Double_t q, Double_t norm, Double_t ymax=0.5, const char *name="fTsallisBW")
Definition: AliPWGFunc.cxx:624
TF1 * GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char *name="fLevi")
Definition: AliPWGFunc.cxx:241
ClassImp(AliPWGFunc) AliPWGFunc
Definition: AliPWGFunc.cxx:29
TF1 * GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char *name="fLevi")
Definition: AliPWGFunc.cxx:821
TF1 * GetFermiDirac(Double_t mass, Double_t T, Double_t norm, const char *name="fFermiDirac")
Definition: AliPWGFunc.cxx:195
TF1 * GetTestFunc(Double_t mass, Double_t T, Double_t norm, Double_t ymax, const char *name="fTest")
Definition: AliPWGFunc.cxx:924
TF1 * GetBoltzmanndNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char *name="fBoltzmannTimesPt")
Definition: AliPWGFunc.cxx:753
static Double_t StaticHistoFunc(const double *x, const double *p)
Definition: AliPWGFunc.cxx:351
static Double_t StaticBGdNdPtTimesPt(const double *x, const double *p)
Definition: AliPWGFunc.cxx:521
static Double_t IntegrandBG(const double *x, const double *p)
Definition: AliPWGFunc.cxx:460
TF1 * GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char *name="fBoltzmann")
Definition: AliPWGFunc.cxx:101
TF1 * GetLevidNdpt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char *name="fLevi")
Definition: AliPWGFunc.cxx:838
static Double_t StaticTest(const double *x, const double *p)
Definition: AliPWGFunc.cxx:899
Definition: External.C:196