Migrating PWG2/SPECTRA/Fit to new PWG structure
[u/mrichter/AliRoot.git] / 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 // }