Migrating PWG2/SPECTRA/Fit to new PWG structure
[u/mrichter/AliRoot.git] / PWG / Tools / AliPWGFunc.cxx
CommitLineData
4c0d7bc7 1//----------------------------------------------------------------------
fef6f2b3 2// AliPWGFunc
4c0d7bc7 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
fef6f2b3 11// variable you want to use calling AliPWGFunc::SetVarType with one of
4c0d7bc7 12// the elements of the VarType_t enum.
13//
fef6f2b3 14// TODO: code cleaup, make the naming convention of function more transparent
15//
8ffc8c1c 16// Warning: not all variables are implemented for all the functions.
4c0d7bc7 17//
18// Author: M. Floris, CERN
19//----------------------------------------------------------------------
20
fef6f2b3 21#include "AliPWGFunc.h"
4c0d7bc7 22#include "TMath.h"
23#include "TF1.h"
24#include "TF3.h"
25#include "TH1.h"
26#include "TSpline.h"
8ffc8c1c 27#include "AliLog.h"
4c0d7bc7 28
fef6f2b3 29ClassImp(AliPWGFunc)
4c0d7bc7 30
fef6f2b3 31AliPWGFunc::AliPWGFunc () : fLastFunc(0), fLineWidth(1), fVarType(kdNdpt) {
4c0d7bc7 32
33 // ctor
34 fLineWidth = 1;
35}
fef6f2b3 36AliPWGFunc::~AliPWGFunc(){
4c0d7bc7 37
38 // dtor
39 if (fLastFunc) delete fLastFunc;
40
41}
42
43
fef6f2b3 44TF1 * AliPWGFunc::GetHistoFunc(TH1 * h, const char * name) {
4c0d7bc7 45
46 // Regardless of the variable type, this returns a function made
47 // from the histo * a multiplicative normalization.
a4097895 48 // This uses a bad hack...
4c0d7bc7 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
a4097895 58
59}
fef6f2b3 60TF1 * AliPWGFunc::GetGraphFunc(TGraph * g, const char * name) {
a4097895 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
4c0d7bc7 74}
75
76
fef6f2b3 77TF1 * AliPWGFunc::GetBGBW(Double_t mass, Double_t beta, Double_t T,
a4097895 78 Double_t n, Double_t norm, const char * name){
4c0d7bc7 79
80 // Boltzmann-Gibbs blast wave
81
82 switch (fVarType) {
83 case kdNdpt:
a4097895 84 return GetBGBWdNdptTimesPt(mass,beta,T,n,norm,name);
4c0d7bc7 85 break;
86 case kOneOverPtdNdpt:
a4097895 87 return GetBGBWdNdpt(mass,beta,T,n,norm,name);
4c0d7bc7 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
fef6f2b3 101TF1 * AliPWGFunc::GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char * name){
4c0d7bc7 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
fef6f2b3 121TF1 * AliPWGFunc::GetTsallisBW(Double_t mass, Double_t beta, Double_t T, Double_t q,
4c0d7bc7 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
fef6f2b3 144TF1 * AliPWGFunc::GetMTExp(Double_t mass, Double_t T, Double_t norm, const char * name){
4c0d7bc7 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
a4097895 166
fef6f2b3 167TF1 * AliPWGFunc::GetBoseEinstein(Double_t mass, Double_t T, Double_t norm, const char * name){
a4097895 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
fef6f2b3 189TF1 * AliPWGFunc::GetFermiDirac(Double_t mass, Double_t T, Double_t norm, const char * name){
a4097895 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
fef6f2b3 212TF1 * AliPWGFunc::GetPTExp(Double_t T, Double_t norm, const char * name){
4c0d7bc7 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
fef6f2b3 235TF1 * AliPWGFunc::GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
4c0d7bc7 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:
0e019ec6 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);
4c0d7bc7 252 break;
253 default:
254 AliFatal("Not implemented");
255 }
256
257 return 0;
258
259
260}
261
fef6f2b3 262TF1 * AliPWGFunc::GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name){
4c0d7bc7 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
fef6f2b3 286TF1 * AliPWGFunc::GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char * name) {
4c0d7bc7 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:
fef6f2b3 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
4c0d7bc7 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
fef6f2b3 345Double_t AliPWGFunc::StaticHistoFunc(const double * x, const double* p){
4c0d7bc7 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
a4097895 352 TObject * h = (TObject*) Long64_t(p[1]);
4c0d7bc7 353
354// Int_t bin = h->FindBin(x[0]);
355// double value = h->GetBinContent(bin);
356
4c0d7bc7 357
58eeda84 358 // static TH1 * oldptr = 0;
359 // static TSpline3 * spl = 0;
360 // if (h!=oldptr) {
361 // FIXME: recheck static pointers
a4097895 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 {
fef6f2b3 369 Printf("AliPWGFunc::StaticHistoFunc: Unsupported type");
a4097895 370 return 0;
371 }
58eeda84 372 // }
4c0d7bc7 373 double value = spl->Eval(x[0]);
374 delete spl;
375
376 return value*norm;
377
378}
379
fef6f2b3 380Double_t AliPWGFunc::StaticUA1Func(const double * x, const double* p) {
4c0d7bc7 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];
8ffc8c1c 388 Double_t temp = p[4];
4c0d7bc7 389 Double_t norm = p[5];
390
391 Double_t xx = x[0];
392
fef6f2b3 393 static AliPWGFunc * self = new AliPWGFunc;
4c0d7bc7 394 static TF1 * fPLaw = self->GetPowerLawdNdptTimesPt(pt0, n, norm, "fLocalPLawUA1");
8ffc8c1c 395 static TF1 * fPMTExp = self->GetMTExpdNdptTimesPt (mass, temp, norm, "fLocalMTexpUA1");
4c0d7bc7 396
397 fPLaw->SetParameters(norm,pt0,n);
8ffc8c1c 398 fPMTExp->SetParameters(1,temp);
4c0d7bc7 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 ) {
fef6f2b3 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}
416Double_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") ;
4c0d7bc7 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
fef6f2b3 454Double_t AliPWGFunc::IntegrandBG(const double * x, const double* p){
4c0d7bc7 455 // integrand for boltzman-gibbs blast wave
a4097895 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
4c0d7bc7 463
464 double x0 = x[0];
465
a4097895 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
4c0d7bc7 476 double mT = TMath::Sqrt(mass*mass+pT*pT);
477
a4097895 478 double rho0 = TMath::ATanH(beta);
8ffc8c1c 479 double arg00 = pT*TMath::SinH(rho0)/temp;
a4097895 480 if (arg00 > 700.) arg00 = 700.; // avoid FPE
8ffc8c1c 481 double arg01 = mT*TMath::CosH(rho0)/temp;
482 double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01);
4c0d7bc7 483
a4097895 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
4c0d7bc7 486 return f0;
487}
488
489
490
fef6f2b3 491Double_t AliPWGFunc::StaticBGdNdPt(const double * x, const double* p) {
4c0d7bc7 492
493 // implementation of BGBW (1/pt dNdpt)
494
495 double pT = x[0];;
496
497
a4097895 498 double mass = p[0];
499 double beta = p[1];
8ffc8c1c 500 double temp = p[2];
a4097895 501 double n = p[3];
502 double norm = p[4];
4c0d7bc7 503
504 static TF1 * fIntBG = 0;
505 if(!fIntBG)
a4097895 506 fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 5);
4c0d7bc7 507
a4097895 508 fIntBG->SetParameters(mass, pT, beta, temp,n);
4c0d7bc7 509 double result = fIntBG->Integral(0,1);
a4097895 510 // printf ("[%4.4f], Int :%f\n", pT, result);
511 return result*norm;//*1e30;;
4c0d7bc7 512
513}
514
fef6f2b3 515Double_t AliPWGFunc::StaticBGdNdPtTimesPt(const double * x, const double* p) {
4c0d7bc7 516 // BGBW dNdpt implementation
517 return x[0]*StaticBGdNdPt(x,p);
518}
519
520
fef6f2b3 521TF1 * AliPWGFunc::GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t temp,
a4097895 522 Double_t n, Double_t norm, const char * name){
4c0d7bc7 523
524 // BGBW 1/pt dNdpt
525
a4097895 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");
4c0d7bc7 530 fLastFunc->SetLineWidth(fLineWidth);
531 return fLastFunc;
532
533}
534
535
536//_____________________________________________________________________
537// Tsallis
538
fef6f2b3 539Double_t AliPWGFunc::IntegrandTsallis(const double * x, const double* p){
4c0d7bc7 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];
8ffc8c1c 550 Double_t temp = p[3];
4c0d7bc7 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( (
8ffc8c1c 559 1+(q-1)/temp * (
4c0d7bc7 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
fef6f2b3 573Double_t AliPWGFunc::StaticTsallisdNdPt(const double * x, const double* p) {
4c0d7bc7 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];
8ffc8c1c 582 double temp = p[2];
4c0d7bc7 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
8ffc8c1c 596 fInt->SetParameters(mass, pT, beta, temp, q);
4c0d7bc7 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
fef6f2b3 604Double_t AliPWGFunc::StaticTsallisdNdPtTimesPt(const double * x, const double* p) {
4c0d7bc7 605
606 // tsallis BW , implementatio of dNdpt
607 return x[0]*StaticTsallisdNdPt(x,p);
608
609}
610
fef6f2b3 611TF1 * AliPWGFunc::GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t temp, Double_t q,
4c0d7bc7 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);
8ffc8c1c 618 fLastFunc->SetParameters(mass,beta,temp,q,norm,ymax);
4c0d7bc7 619 fLastFunc->SetParLimits(1,0.0,0.99);
620 fLastFunc->SetParLimits(2,0.01,0.99);
621 fLastFunc->SetParLimits(3,1.0001,1.9);
8ffc8c1c 622 fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax");
4c0d7bc7 623 fLastFunc->SetLineWidth(fLineWidth);
624 return fLastFunc;
625
626}
627
628// Times Pt funcs
629// Boltzmann-Gibbs Blast Wave
fef6f2b3 630TF1 * AliPWGFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t n,
4c0d7bc7 631 Double_t norm, const char * name){
632
633 // BGBW, dNdpt
634
a4097895 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");
4c0d7bc7 639 fLastFunc->SetLineWidth(fLineWidth);
640 return fLastFunc;
641
642
643}
644
645
646
fef6f2b3 647TF1 * AliPWGFunc::GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t q,
4c0d7bc7 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);
8ffc8c1c 653 fLastFunc->SetParameters(mass,beta,temp,q,norm,ymax);
654 fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax");
4c0d7bc7 655 fLastFunc->SetLineWidth(fLineWidth);
656 return fLastFunc;
657
658
659}
660
661
662
fef6f2b3 663TF1 * AliPWGFunc::GetMTExpdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
4c0d7bc7 664
665 // Simple exponential in 1/mt*MT, as a function of dNdpt
666 char formula[500];
29b53d49 667 snprintf(formula,500,"[0]*x*exp(-sqrt(x**2+%f**2)/[1])", mass);
4c0d7bc7 668 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 669 fLastFunc->SetParameters(norm, temp);
4c0d7bc7 670 fLastFunc->SetParLimits(1, 0.01, 10);
671 fLastFunc->SetParNames("norm", "T");
672 fLastFunc->SetLineWidth(fLineWidth);
673 return fLastFunc;
674
675
676}
677
fef6f2b3 678TF1 * AliPWGFunc::GetBoseEinsteindNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
a4097895 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
fef6f2b3 693TF1 * AliPWGFunc::GetFermiDiracdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
a4097895 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
4c0d7bc7 709
fef6f2b3 710TF1 * AliPWGFunc::GetPTExpdNdptTimesPt(Double_t temp, Double_t norm, const char * name){
4c0d7bc7 711
712 // Simple exponential in 1/pt*dNdpT, as a function of dNdpt
713 char formula[500];
29b53d49 714 snprintf(formula,500,"[0]*x*exp(-x/[1])");
4c0d7bc7 715 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 716 fLastFunc->SetParameters(norm, temp);
4c0d7bc7 717 fLastFunc->SetParLimits(1, 0.01, 10);
718 fLastFunc->SetParNames("norm", "T");
719 fLastFunc->SetLineWidth(fLineWidth);
720 return fLastFunc;
721
722
723}
724
725
fef6f2b3 726TF1 * AliPWGFunc::GetBoltzmanndNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
4c0d7bc7 727 // Boltzmann (exp in 1/mt*dNdmT times mt) as a function of dNdpt
728 char formula[500];
29b53d49 729 snprintf(formula,500,"[0]*x*sqrt(x**2+%f**2)*exp(-sqrt(x**2+%f**2)/[1])", mass,mass);
4c0d7bc7 730 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 731 fLastFunc->SetParameters(norm, temp);
4c0d7bc7 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)
fef6f2b3 742// TF1 * AliPWGFunc::GetTsallisdNdptTimesPt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name){
4c0d7bc7 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
fef6f2b3 759TF1 * AliPWGFunc::GetLevidNdptTimesPt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
4c0d7bc7 760
761 // Levi function, dNdpt
762 char formula[500];
763
29b53d49 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])");
4c0d7bc7 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);
8ffc8c1c 767 fLastFunc->SetParameters(norm, n, temp,mass);
4c0d7bc7 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
fef6f2b3 777TF1 * AliPWGFunc::GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name){
4c0d7bc7 778
779 // PowerLaw function, dNdpt
780 char formula[500];
781
29b53d49 782 snprintf(formula,500,"x*[0]*( 1 + x/[1] )^(-[2])");
4c0d7bc7 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
fef6f2b3 794TF1 * AliPWGFunc::GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name){
4c0d7bc7 795
796 // PowerLaw function, 1/pt dNdpt
797 char formula[500];
798
29b53d49 799 snprintf(formula,500," [0]*( 1 + x/[1] )^(-[2])");
4c0d7bc7 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
fef6f2b3 811TF1 * AliPWGFunc::GetLevidNdpt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
4c0d7bc7 812
813 // Levi function, dNdpt
814 char formula[500];
815
29b53d49 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])");
4c0d7bc7 817 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 818 fLastFunc->SetParameters(norm, n, temp,mass);
4c0d7bc7 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
fef6f2b3 828TF1 * AliPWGFunc::GetLevidNdmt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name, VarType_t var){
4c0d7bc7 829
0e019ec6 830 // Levi function, 1/mt dNdmt
4c0d7bc7 831 char formula[500];
0e019ec6 832 if (var == kOneOverMtdNdmt)
29b53d49 833 snprintf(formula,500,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x -[3])/([1]*[2]) )^(-[1])");
0e019ec6 834 else if (var == kdNdmt)
29b53d49 835 snprintf(formula,500,"( x*[0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x-[3])/([1]*[2]) )^(-[1])");
0e019ec6 836 if (var == kOneOverMtdNdmtMinusM)
29b53d49 837 snprintf(formula,500,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x)/([1]*[2]) )^(-[1])");
0e019ec6 838
839 //sprintf(formula,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + x/([1]*[2]) )^(-[1])");
4c0d7bc7 840 // sprintf(formula,"[0] * ( 1 + x/([1]*[2]) )^(-[1])");
841 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 842 fLastFunc->SetParameters(norm, n, temp,mass);
4c0d7bc7 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
0e019ec6 854
4c0d7bc7 855// Test Function
fef6f2b3 856Double_t AliPWGFunc::IntegrandTest(const double * x, const double* p){
4c0d7bc7 857
858 // test function
859
860 Double_t y = x[0];
861
862 Double_t mass = p[0];
863 Double_t pt = p[1];
8ffc8c1c 864 Double_t temp = p[2];
4c0d7bc7 865
866 Double_t mt = TMath::Sqrt(mass*mass+pt*pt);
867
8ffc8c1c 868 return mt*TMath::CosH(y)*TMath::Exp(-mt*TMath::CosH(y)/temp);
4c0d7bc7 869
870}
871
fef6f2b3 872Double_t AliPWGFunc::StaticTest(const double * x, const double* p) {
4c0d7bc7 873
874 // test function
875
876 double pT = x[0];;
877
878
879 double mass = p[0];
8ffc8c1c 880 double temp = p[1];
4c0d7bc7 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
8ffc8c1c 890 fIntTest->SetParameters(mass, pT, temp);
4c0d7bc7 891 double result = fIntTest->Integral(-ymax, ymax);
892
893 return result*p[2];//*1e30;;
894
895}
896
fef6f2b3 897TF1 * AliPWGFunc::GetTestFunc(Double_t mass, Double_t temp, Double_t norm, Double_t ymax, const char * name){
4c0d7bc7 898
899 // test function
900
901 fLastFunc = new TF1 (name, StaticTest, 0.0, 10, 4);
8ffc8c1c 902 fLastFunc->SetParameters(mass,temp,norm,ymax);
4c0d7bc7 903 fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
904 fLastFunc->SetLineWidth(fLineWidth);
905 return fLastFunc;
906
907}
908
909
910//___________________________________________________________
911
912
fef6f2b3 913TF1 * AliPWGFunc::GetMTExpdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
4c0d7bc7 914 // Simple exp in 1/mt dNdmt, as a function of dNdpt
915 // mt scaling
916 char formula[500];
29b53d49 917 snprintf(formula,500,"[0]*exp(-sqrt(x**2+%f**2)/[1])", mass);
4c0d7bc7 918 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 919 fLastFunc->SetParameters(norm, temp);
4c0d7bc7 920 fLastFunc->SetParLimits(1, 0.01, 10);
921 fLastFunc->SetParNames("norm", "T");
922 fLastFunc->SetLineWidth(fLineWidth);
923 return fLastFunc;
924}
925
fef6f2b3 926TF1 * AliPWGFunc::GetBoseEinsteindNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
a4097895 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
fef6f2b3 938TF1 * AliPWGFunc::GetFermiDiracdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
a4097895 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
4c0d7bc7 950
951// // Simple tsallis (a la CMS)
fef6f2b3 952// TF1 * AliPWGFunc::GetTsallisdNdpt(Double_t mass, Double_t temp, Double_t q, Double_t norm, const char * name){
4c0d7bc7 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);
8ffc8c1c 957// fLastFunc->SetParameters(norm, temp, q);
4c0d7bc7 958// fLastFunc->SetParLimits(1, 0.01, 10);
959// fLastFunc->SetParNames("norm", "T", "q");
960// fLastFunc->SetLineWidth(fLineWidth);
961// return fLastFunc;
962// }