Removing try/catch-es
[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   // 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 // }