Fixing coding conventions.
[u/mrichter/AliRoot.git] / 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
47   fLastFunc = new TF1 (name, StaticHistoFunc, 0.0, 10, 2);
48   fLastFunc->SetParameter(0,1);
49   fLastFunc->FixParameter(1,Double_t(Long64_t(h)));
50   fLastFunc->SetParNames("norm", "pointer to histo");
51   fLastFunc->SetLineWidth(fLineWidth);
52   return fLastFunc;
53   
54
55 }
56
57
58 TF1 * AliBWFunc::GetBGBW(Double_t mass, Double_t beta, Double_t T,
59               Double_t norm, const char * name){
60
61   // Boltzmann-Gibbs blast wave
62
63   switch (fVarType) {
64   case kdNdpt:
65     return GetBGBWdNdptTimesPt(mass,beta,T,norm,name);
66     break;
67   case kOneOverPtdNdpt:
68     return GetBGBWdNdpt(mass,beta,T,norm,name);
69     break;
70   case kOneOverMtdNdmt:
71     AliFatal("Not implemented");
72     break;
73   default:
74     AliFatal("Not implemented");
75   }
76   
77   return 0;
78
79 }
80   
81
82 TF1 * AliBWFunc::GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char * name){
83   // Boltzmann
84   switch (fVarType) {
85   case kdNdpt:
86     return GetBoltzmanndNdptTimesPt(mass, T, norm, name);
87   case kOneOverPtdNdpt:
88     AliFatal("Not implemented");
89     break;
90   case kOneOverMtdNdmt:
91     AliFatal("Not implemented");
92     break;
93   default:
94     AliFatal("Not implemented");
95   }
96   
97   return 0;
98
99 }
100
101
102 TF1 * AliBWFunc::GetTsallisBW(Double_t mass, Double_t beta, Double_t T, Double_t q,
103                               Double_t norm, Double_t ymax, const char * name){
104
105   // Tsallis blast wave
106   switch (fVarType) {
107   case kdNdpt:
108     return GetTsallisBWdNdptTimesPt(mass,beta,T,q,norm,ymax,name);
109     break;
110   case kOneOverPtdNdpt:
111     return GetTsallisBWdNdpt(mass,beta,T,q,norm,ymax,name);
112     break;
113   case kOneOverMtdNdmt:
114     AliFatal("Not implemented");
115     break;
116   default:
117     AliFatal("Not implemented");
118   }
119   
120   return 0;
121   
122 }
123
124
125 TF1 * AliBWFunc::GetMTExp(Double_t mass, Double_t T, Double_t norm, const char * name){
126
127   // Simple exponential in 1/mt*MT
128   switch (fVarType) {
129   case kdNdpt:
130     return GetMTExpdNdptTimesPt(mass,T,norm,name);
131     break;
132   case kOneOverPtdNdpt:
133     return GetMTExpdNdpt(mass,T,norm,name);
134     break;
135   case kOneOverMtdNdmt:
136     AliFatal("Not implemented");
137     break;
138   default:
139     AliFatal("Not implemented");
140   }
141   
142   return 0;
143
144
145 }
146
147 TF1 * AliBWFunc::GetPTExp(Double_t T, Double_t norm, const char * name){
148
149   // Simple exponential in 1/mt*MT
150   switch (fVarType) {
151   case kdNdpt:
152     return GetPTExpdNdptTimesPt(T,norm,name);
153     break;
154   case kOneOverPtdNdpt:
155     AliFatal("Not implemented");
156     break;
157   case kOneOverMtdNdmt:
158     AliFatal("Not implemented");
159     break;
160   default:
161     AliFatal("Not implemented");
162   }
163   
164   return 0;
165
166
167 }
168
169
170 TF1 * AliBWFunc::GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
171   // Levi function (aka Tsallis)
172   switch (fVarType) {
173   case kdNdpt:
174     return GetLevidNdptTimesPt(mass,T,n,norm,name);
175     break;
176   case kOneOverPtdNdpt:
177     return GetLevidNdpt(mass,T,n,norm,name);
178     break;
179   case kOneOverMtdNdmt:
180     return GetLevidNdmt(mass,T,n,norm,name);
181     break;
182   default:
183     AliFatal("Not implemented");
184   }
185   
186   return 0;
187
188
189 }
190
191 TF1 * AliBWFunc::GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name){
192   // power law Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
193   // This is sometimes also called Hagedorn or modified Hagedorn
194
195   switch (fVarType) {
196   case kdNdpt:
197     return GetPowerLawdNdptTimesPt(pt0,n,norm,name);
198     break;
199   case kOneOverPtdNdpt:
200     return GetPowerLawdNdpt(pt0,n,norm,name);
201     break;
202   case kOneOverMtdNdmt:
203     AliFatal("Not Implemented");
204     //    return GetUA1dNdmt(mass,T,n,norm,name);
205     break;
206   default:
207     AliFatal("Not implemented");
208   }
209   
210   return 0;
211
212
213 }
214
215 TF1 * AliBWFunc::GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char * name) {
216   // UA1 parametrization Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
217
218   switch (fVarType) {
219   case kdNdpt:
220
221     fLastFunc = new TF1 (name, StaticUA1Func, 0.0, 10, 6);
222     fLastFunc->FixParameter(0,mass);
223     fLastFunc->SetParameter(1,p0star);
224     fLastFunc->SetParameter(2,pt0);
225     fLastFunc->SetParameter(3,n);
226     fLastFunc->SetParameter(4,T);
227     fLastFunc->SetParameter(5,norm);
228     fLastFunc->SetParLimits(1,0.01,1);
229     fLastFunc->SetParLimits(2,0.01,100);
230     fLastFunc->SetParLimits(3,0.01,100);
231     fLastFunc->SetParLimits(4,0.01,100);
232     fLastFunc->SetParNames("mass","p0star","pt0","n","T","norm");
233     fLastFunc->SetNpx(5000);
234     fLastFunc->SetLineWidth(fLineWidth);
235     return fLastFunc;
236
237     break;
238   case kOneOverPtdNdpt:
239     AliFatal("Not Implemented");
240     break;
241   case kOneOverMtdNdmt:
242     AliFatal("Not Implemented");
243     //    return GetUA1dNdmt(mass,T,n,norm,name);
244     break;
245   default:
246     AliFatal("Not implemented");
247   }
248
249   return 0;
250 }
251
252
253
254
255 // ________________________________________________________________________
256
257 // Backend (private functions and support functions for numerical integration)
258
259 Double_t AliBWFunc::StaticHistoFunc(const double * x, const double* p){
260
261   // provides a function interpolating a histo with a spline; 
262   // using double to store a pointer... This is a bad hack. To be replaced
263
264   double norm = p[0];
265   
266   TH1 * h     = (TH1*) Long64_t(p[1]);
267
268 //    Int_t bin = h->FindBin(x[0]);
269 //    double value = h->GetBinContent(bin);
270
271   if (h->FindBin(x[0]) > h->GetNbinsX()) return 0;
272
273   static TH1 * oldptr = 0;
274   TSpline3 * spl = 0;
275   if (h!=oldptr) {
276     spl  = new TSpline3(h);
277   }
278   double value = spl->Eval(x[0]);
279   delete spl;
280
281   return value*norm;
282   
283 }
284
285 Double_t AliBWFunc::StaticUA1Func(const double * x, const double* p) {
286   
287
288   // "mass","p0star","pt0","n","T","norm"
289   Double_t mass   = p[0];
290   Double_t p0star = p[1];
291   Double_t pt0    = p[2];
292   Double_t n      = p[3];
293   Double_t temp   = p[4];
294   Double_t norm   = p[5];
295   
296   Double_t xx = x[0];
297
298   static AliBWFunc * self = new AliBWFunc;
299   static TF1 * fPLaw   = self->GetPowerLawdNdptTimesPt(pt0, n, norm, "fLocalPLawUA1");
300   static TF1 * fPMTExp = self->GetMTExpdNdptTimesPt   (mass, temp, norm, "fLocalMTexpUA1");
301
302   fPLaw->SetParameters(norm,pt0,n);
303   fPMTExp->SetParameters(1,temp);
304   
305
306   Double_t normMT =fPMTExp->Eval(p0star) >0 ? fPLaw->Eval(p0star) / fPMTExp->Eval(p0star) *  fPMTExp->GetParameter(0) : 1;
307   fPMTExp->SetParameter(0,normMT);
308   
309   
310   if (TMath::Abs(fPMTExp->Eval(p0star) - fPLaw->Eval(p0star)) > 0.0001 ) {
311     Printf("AliBWFunc::StaticUA1Func - Wrong norm") ; 
312     Printf(" p0* %f  NMT: %f  N: %f  PL: %f  MT: %f", p0star, normMT, norm, fPLaw->Eval(p0star), fPMTExp->Eval(p0star));
313   }
314
315   if (xx > p0star)  return fPLaw->Eval(xx);
316   return fPMTExp->Eval(xx);    
317   
318   
319
320 }
321
322
323 Double_t AliBWFunc::IntegrandBG(const double * x, const double* p){
324   // integrand for boltzman-gibbs blast wave
325
326   double x0 = x[0]; 
327   
328   double mass = p[0];
329   double pT   = p[1];
330   double beta = p[2];
331   double temp    = p[3];
332   
333   double mT      = TMath::Sqrt(mass*mass+pT*pT);
334
335   double rho0   = TMath::ATanH(beta*x0);  
336   double arg00 = pT*TMath::SinH(rho0)/temp;
337   double arg01 = mT*TMath::CosH(rho0)/temp;
338   double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01);
339
340   return f0;
341 }
342
343
344
345 Double_t AliBWFunc::StaticBGdNdPt(const double * x, const double* p) {
346
347   // implementation of BGBW (1/pt dNdpt)
348
349   double pT = x[0];;
350   
351
352   double mass = p[0];
353   double beta = p[1];
354   double temp    = p[2];
355
356   static TF1 * fIntBG = 0;
357   if(!fIntBG)
358     fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 4);
359
360   fIntBG->SetParameters(mass, pT, beta, temp);
361   double result = fIntBG->Integral(0,1);
362   return result*p[3];//*1e30;;
363
364 }
365
366 Double_t AliBWFunc::StaticBGdNdPtTimesPt(const double * x, const double* p) {
367   // BGBW dNdpt implementation
368   return x[0]*StaticBGdNdPt(x,p);
369 }
370
371
372 TF1 * AliBWFunc::GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t temp,
373                             Double_t norm, const char * name){
374   
375   // BGBW 1/pt dNdpt
376
377   fLastFunc = new TF1 (name, StaticBGdNdPt, 0.0, 10, 4);
378   fLastFunc->SetParameters(mass,beta,temp,norm);    
379   fLastFunc->SetParNames("mass", "#beta", "T", "norm");
380   fLastFunc->SetLineWidth(fLineWidth);
381   return fLastFunc;
382   
383 }
384
385
386 //_____________________________________________________________________
387 // Tsallis
388
389 Double_t AliBWFunc::IntegrandTsallis(const double * x, const double* p){
390
391   // integrand for numerical integration (tsallis)
392
393   Double_t r   = x[0]; 
394   Double_t phi = x[1];
395   Double_t y   = x[2];
396
397   Double_t mass = p[0];
398   Double_t pt   = p[1];
399   Double_t beta = p[2];
400   Double_t temp    = p[3];
401   Double_t q    = p[4];
402   
403   Double_t mt      = TMath::Sqrt(mass*mass+pt*pt);
404
405   Double_t rho    = TMath::ATanH(beta*r); // TODO: implement different velocity profiles  
406
407   Double_t res = mt*
408     r*TMath::CosH(y) *TMath::Power( (
409                                      1+(q-1)/temp * (
410                                                   mt*TMath::CosH(y)*TMath::CosH(rho) -
411                                                   pt*TMath::SinH(rho)*TMath::Cos(phi)
412                                                   )
413                                      ),
414                                        -1/(q-1)
415                                     );                  
416
417
418   return res;
419 }
420
421
422
423 Double_t AliBWFunc::StaticTsallisdNdPt(const double * x, const double* p) {
424
425   // tsallis BW implementation 1/pt dNdpt
426
427   double pT = x[0];;
428   
429
430   double mass = p[0];
431   double beta = p[1];
432   double temp    = p[2];
433   double q    = p[3];
434
435   Double_t ymax = p[5];
436
437
438   static TF3 * fInt = 0;
439   if(!fInt){
440     fInt = new TF3 ("fIntTsa", IntegrandTsallis, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
441 //     fInt->SetNpx(10000);
442 //     fInt->SetNpy(10000);
443 //     fInt->SetNpz(10000);
444   }
445   
446   fInt->SetParameters(mass, pT, beta, temp, q);
447   double result = fInt->Integral(0,1, -TMath::Pi(), TMath::Pi(), -ymax, ymax);
448   //  double result = fInt->Integral(0,1, -2, 2, -ymax, ymax);
449   
450   return result*p[4];//*1e30;;
451
452 }
453
454 Double_t AliBWFunc::StaticTsallisdNdPtTimesPt(const double * x, const double* p) {
455
456   // tsallis BW , implementatio of dNdpt
457   return x[0]*StaticTsallisdNdPt(x,p);
458
459 }
460
461 TF1 * AliBWFunc::GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t temp, Double_t q,
462                                    Double_t norm, Double_t ymax,const char * name){
463   
464
465   // tsallis BW, 1/pt dNdpt
466
467   fLastFunc = new TF1 (name, StaticTsallisdNdPt, 0.0, 10, 6);
468   fLastFunc->SetParameters(mass,beta,temp,q,norm,ymax);
469   fLastFunc->SetParLimits(1,0.0,0.99);
470   fLastFunc->SetParLimits(2,0.01,0.99);
471   fLastFunc->SetParLimits(3,1.0001,1.9);
472   fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax");
473   fLastFunc->SetLineWidth(fLineWidth);
474   return fLastFunc;
475   
476 }
477
478 // Times Pt funcs
479 // Boltzmann-Gibbs Blast Wave
480 TF1 * AliBWFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp,
481                                      Double_t norm, const char * name){
482
483   // BGBW, dNdpt
484
485   fLastFunc = new TF1 (name, StaticBGdNdPtTimesPt, 0.0, 10, 4);
486   fLastFunc->SetParameters(mass,beta,temp,norm);    
487   fLastFunc->SetParNames("mass", "#beta", "temp", "norm");
488   fLastFunc->SetLineWidth(fLineWidth);
489   return fLastFunc;
490
491
492 }
493
494
495
496 TF1 * AliBWFunc::GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t q,
497                                           Double_t norm, Double_t ymax, const char * name){
498
499 // Tsallis blast wave, dNdpt
500
501   fLastFunc = new TF1 (name, StaticTsallisdNdPtTimesPt, 0.0, 10, 6);
502   fLastFunc->SetParameters(mass,beta,temp,q,norm,ymax);    
503   fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax");
504   fLastFunc->SetLineWidth(fLineWidth);
505   return fLastFunc;
506  
507
508 }
509
510
511
512 TF1 * AliBWFunc::GetMTExpdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
513
514   // Simple exponential in 1/mt*MT, as a function of dNdpt
515   char formula[500];
516   sprintf(formula,"[0]*x*exp(-sqrt(x**2+%f**2)/[1])", mass);
517   fLastFunc=new TF1(name,formula,0,10);
518   fLastFunc->SetParameters(norm, temp);
519   fLastFunc->SetParLimits(1, 0.01, 10);
520   fLastFunc->SetParNames("norm", "T");
521   fLastFunc->SetLineWidth(fLineWidth);
522   return fLastFunc;
523
524
525 }
526
527
528 TF1 * AliBWFunc::GetPTExpdNdptTimesPt(Double_t temp, Double_t norm, const char * name){
529
530   // Simple exponential in 1/pt*dNdpT, as a function of dNdpt
531   char formula[500];
532   sprintf(formula,"[0]*x*exp(-x/[1])");
533   fLastFunc=new TF1(name,formula,0,10);
534   fLastFunc->SetParameters(norm, temp);
535   fLastFunc->SetParLimits(1, 0.01, 10);
536   fLastFunc->SetParNames("norm", "T");
537   fLastFunc->SetLineWidth(fLineWidth);
538   return fLastFunc;
539
540
541 }
542
543
544 TF1 * AliBWFunc::GetBoltzmanndNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
545   // Boltzmann (exp in 1/mt*dNdmT times mt) as a function of dNdpt
546  char formula[500];
547  sprintf(formula,"[0]*x*sqrt(x**2+%f**2)*exp(-sqrt(x**2+%f**2)/[1])", mass,mass);
548   fLastFunc=new TF1(name,formula,0,10);
549   fLastFunc->SetParameters(norm, temp);
550   fLastFunc->SetParLimits(1, 0.01, 10);
551   fLastFunc->SetParNames("norm", "T");
552   fLastFunc->SetLineWidth(fLineWidth);
553   return fLastFunc;
554
555
556 }
557
558
559 // Tsallis (no BW, a la CMS)
560 // TF1 * AliBWFunc::GetTsallisdNdptTimesPt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name){
561
562 //   char formula[500];
563 //   //  sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2)-%f)),(-1/([2]-1)))", mass, mass); //CMS
564 //   sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass);  // STAR
565 //   //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
566 //   fLastFunc=new TF1(name,formula,0,10);
567 //   fLastFunc->SetParameters(norm, T, q);
568 //   fLastFunc->SetParLimits(1, 0.001, 10);
569 //   fLastFunc->SetParNames("norm", "T", "q");
570 //   fLastFunc->SetLineWidth(fLineWidth);
571 //   return fLastFunc;
572
573
574 // }
575
576
577 TF1 * AliBWFunc::GetLevidNdptTimesPt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
578
579   // Levi function, dNdpt
580   char formula[500];
581
582   sprintf(formula,"( x*[0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2])  )^(-[1])");
583   //  sprintf(formula,"( x*[0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x))/([1]*[2])  )^(-[1])");
584   fLastFunc=new TF1(name,formula,0,10);
585   fLastFunc->SetParameters(norm, n, temp,mass);
586   fLastFunc->SetParLimits(2, 0.01, 10);
587   fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
588   fLastFunc->FixParameter(3,mass);
589   fLastFunc->SetLineWidth(fLineWidth);
590   return fLastFunc;
591
592
593 }
594
595 TF1 * AliBWFunc::GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name){
596
597   // PowerLaw function, dNdpt
598   char formula[500];
599
600   sprintf(formula,"x*[0]*( 1 + x/[1] )^(-[2])");
601   fLastFunc=new TF1(name,formula,0,10);
602   fLastFunc->SetParameters(norm, pt0, n);
603   fLastFunc->SetParLimits(1, 0.01, 10);
604   //fLastFunc->SetParLimits(2, 0.01, 50);
605   fLastFunc->SetParNames("norm", "pt0", "n");
606   fLastFunc->SetLineWidth(fLineWidth);
607   return fLastFunc;
608
609
610 }
611
612 TF1 * AliBWFunc::GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name){
613
614   // PowerLaw function, 1/pt dNdpt
615   char formula[500];
616
617   sprintf(formula," [0]*( 1 + x/[1] )^(-[2])");
618   fLastFunc=new TF1(name,formula,0,10);
619   fLastFunc->SetParameters(norm, pt0, n);
620   //  fLastFunc->SetParLimits(2, 0.01, 10);
621   fLastFunc->SetParNames("norm", "pt0", "n");
622   fLastFunc->SetLineWidth(fLineWidth);
623   return fLastFunc;
624
625
626 }
627
628
629 TF1 * AliBWFunc::GetLevidNdpt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
630
631   // Levi function, dNdpt
632   char formula[500];
633
634   sprintf(formula,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2])  )^(-[1])");
635   fLastFunc=new TF1(name,formula,0,10);
636   fLastFunc->SetParameters(norm, n, temp,mass);
637   fLastFunc->SetParLimits(2, 0.01, 10);
638   fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
639   fLastFunc->FixParameter(3,mass);
640   fLastFunc->SetLineWidth(fLineWidth);
641   return fLastFunc;
642
643
644 }
645
646 TF1 * AliBWFunc::GetLevidNdmt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
647
648   // Levi function, dNdmt
649   char formula[500];
650
651   //  sprintf(formula,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (x -[3])/([1]*[2])  )^(-[1])");
652   sprintf(formula,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + x/([1]*[2])  )^(-[1])");
653   //  sprintf(formula,"[0] * ( 1 + x/([1]*[2])  )^(-[1])");
654   fLastFunc=new TF1(name,formula,0,10);
655   fLastFunc->SetParameters(norm, n, temp,mass);
656   fLastFunc->SetParLimits(2, 0.01, 10);
657   fLastFunc->SetParNames("norm", "n", "T", "mass");
658   fLastFunc->FixParameter(3,mass);
659   fLastFunc->SetLineWidth(fLineWidth);
660   return fLastFunc;
661
662
663 }
664
665
666
667 // Test Function
668 Double_t AliBWFunc::IntegrandTest(const double * x, const double* p){
669
670   // test function
671
672   Double_t y = x[0];
673
674   Double_t mass = p[0];
675   Double_t pt   = p[1];
676   Double_t temp    = p[2];
677
678   Double_t mt      = TMath::Sqrt(mass*mass+pt*pt);    
679   
680   return mt*TMath::CosH(y)*TMath::Exp(-mt*TMath::CosH(y)/temp);
681
682 }
683
684 Double_t AliBWFunc::StaticTest(const double * x, const double* p) {
685
686   // test function
687
688   double pT = x[0];;
689   
690
691   double mass = p[0];
692   double temp    = p[1];
693   Double_t ymax = p[3];
694
695
696   static TF3 * fIntTest = 0;
697   if(!fIntTest){
698     fIntTest = new TF3 ("fIntTest", IntegrandTest, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
699     //    fInt->SetNpx(10000);
700   }
701   
702   fIntTest->SetParameters(mass, pT, temp);
703   double result = fIntTest->Integral(-ymax, ymax);
704   
705   return result*p[2];//*1e30;;
706
707 }
708
709 TF1 * AliBWFunc::GetTestFunc(Double_t mass, Double_t temp, Double_t norm, Double_t ymax, const char * name){
710   
711   // test function
712   
713   fLastFunc = new TF1 (name, StaticTest, 0.0, 10, 4);
714   fLastFunc->SetParameters(mass,temp,norm,ymax);    
715   fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
716   fLastFunc->SetLineWidth(fLineWidth);
717   return fLastFunc;
718   
719 }
720
721
722 //___________________________________________________________
723
724
725 TF1 * AliBWFunc::GetMTExpdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
726   // Simple exp in 1/mt dNdmt, as a function of dNdpt
727   // mt scaling
728   char formula[500];
729   sprintf(formula,"[0]*exp(-sqrt(x**2+%f**2)/[1])", 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 // // Simple tsallis (a la CMS)
740 // TF1 * AliBWFunc::GetTsallisdNdpt(Double_t mass, Double_t temp, Double_t q, Double_t norm, const char * name){
741   
742 //   char formula[500];
743 //   sprintf(formula,"[0]*sqrt(x**2+%f**2)*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass,mass); 
744 //   fLastFunc=new TF1(name,formula,0,10);
745 //   fLastFunc->SetParameters(norm, temp, q);
746 //   fLastFunc->SetParLimits(1, 0.01, 10);
747 //   fLastFunc->SetParNames("norm", "T", "q");
748 //   fLastFunc->SetLineWidth(fLineWidth);
749 //   return fLastFunc;
750 // }