]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/Tools/AliPWGFunc.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[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:
4fd89c20 155 return GetMTExpdNdmt(mass,T,norm,name,kOneOverMtdNdmt);
156 break;
157 case kdNdmt:
158 return GetMTExpdNdmt(mass,T,norm,name,kdNdmt);
159 break;
160 case kOneOverMtdNdmtMinusM:
161 return GetMTExpdNdmt(mass,T,norm,name,kOneOverMtdNdmtMinusM);
4c0d7bc7 162 break;
163 default:
164 AliFatal("Not implemented");
165 }
166
167 return 0;
168
169
170}
171
a4097895 172
fef6f2b3 173TF1 * AliPWGFunc::GetBoseEinstein(Double_t mass, Double_t T, Double_t norm, const char * name){
a4097895 174
175 // Bose einstein
176 switch (fVarType) {
177 case kdNdpt:
178 return GetBoseEinsteindNdptTimesPt(mass,T,norm,name);
179 break;
180 case kOneOverPtdNdpt:
181 return GetBoseEinsteindNdpt(mass,T,norm,name);
182 break;
183 case kOneOverMtdNdmt:
184 AliFatal("Not implemented");
185 break;
186 default:
187 AliFatal("Not implemented");
188 }
189
190 return 0;
191
192
193}
194
fef6f2b3 195TF1 * AliPWGFunc::GetFermiDirac(Double_t mass, Double_t T, Double_t norm, const char * name){
a4097895 196
197 // Simple exponential in 1/mt*MT
198 switch (fVarType) {
199 case kdNdpt:
200 return GetFermiDiracdNdptTimesPt(mass,T,norm,name);
201 break;
202 case kOneOverPtdNdpt:
203 return GetFermiDiracdNdpt(mass,T,norm,name);
204 break;
205 case kOneOverMtdNdmt:
206 AliFatal("Not implemented");
207 break;
208 default:
209 AliFatal("Not implemented");
210 }
211
212 return 0;
213
214
215}
216
217
fef6f2b3 218TF1 * AliPWGFunc::GetPTExp(Double_t T, Double_t norm, const char * name){
4c0d7bc7 219
220 // Simple exponential in 1/mt*MT
221 switch (fVarType) {
222 case kdNdpt:
223 return GetPTExpdNdptTimesPt(T,norm,name);
224 break;
225 case kOneOverPtdNdpt:
226 AliFatal("Not implemented");
227 break;
228 case kOneOverMtdNdmt:
229 AliFatal("Not implemented");
230 break;
231 default:
232 AliFatal("Not implemented");
233 }
234
235 return 0;
236
237
238}
239
240
fef6f2b3 241TF1 * AliPWGFunc::GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
4c0d7bc7 242 // Levi function (aka Tsallis)
243 switch (fVarType) {
244 case kdNdpt:
245 return GetLevidNdptTimesPt(mass,T,n,norm,name);
246 break;
247 case kOneOverPtdNdpt:
248 return GetLevidNdpt(mass,T,n,norm,name);
249 break;
250 case kOneOverMtdNdmt:
0e019ec6 251 return GetLevidNdmt(mass,T,n,norm,name,kOneOverMtdNdmt);
252 break;
253 case kdNdmt:
254 return GetLevidNdmt(mass,T,n,norm,name,kdNdmt);
255 break;
256 case kOneOverMtdNdmtMinusM:
257 return GetLevidNdmt(mass,T,n,norm,name,kOneOverMtdNdmtMinusM);
4c0d7bc7 258 break;
259 default:
260 AliFatal("Not implemented");
261 }
262
263 return 0;
264
265
266}
267
fef6f2b3 268TF1 * AliPWGFunc::GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name){
4c0d7bc7 269 // power law Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
270 // This is sometimes also called Hagedorn or modified Hagedorn
271
272 switch (fVarType) {
273 case kdNdpt:
274 return GetPowerLawdNdptTimesPt(pt0,n,norm,name);
275 break;
276 case kOneOverPtdNdpt:
277 return GetPowerLawdNdpt(pt0,n,norm,name);
278 break;
279 case kOneOverMtdNdmt:
280 AliFatal("Not Implemented");
281 // return GetUA1dNdmt(mass,T,n,norm,name);
282 break;
283 default:
284 AliFatal("Not implemented");
285 }
286
287 return 0;
288
289
290}
291
fef6f2b3 292TF1 * AliPWGFunc::GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char * name) {
4c0d7bc7 293 // UA1 parametrization Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
294
295 switch (fVarType) {
296 case kdNdpt:
297
298 fLastFunc = new TF1 (name, StaticUA1Func, 0.0, 10, 6);
299 fLastFunc->FixParameter(0,mass);
300 fLastFunc->SetParameter(1,p0star);
301 fLastFunc->SetParameter(2,pt0);
302 fLastFunc->SetParameter(3,n);
303 fLastFunc->SetParameter(4,T);
304 fLastFunc->SetParameter(5,norm);
305 fLastFunc->SetParLimits(1,0.01,1);
306 fLastFunc->SetParLimits(2,0.01,100);
307 fLastFunc->SetParLimits(3,0.01,100);
308 fLastFunc->SetParLimits(4,0.01,100);
309 fLastFunc->SetParNames("mass","p0star","pt0","n","T","norm");
310 fLastFunc->SetNpx(5000);
311 fLastFunc->SetLineWidth(fLineWidth);
312 return fLastFunc;
313
314 break;
315 case kOneOverPtdNdpt:
fef6f2b3 316 fLastFunc = new TF1 (name, StaticUA1FuncOneOverPt, 0.0, 10, 6);
317 fLastFunc->FixParameter(0,mass);
318 fLastFunc->SetParameter(1,p0star);
319 fLastFunc->SetParameter(2,pt0);
320 fLastFunc->SetParameter(3,n);
321 fLastFunc->SetParameter(4,T);
322 fLastFunc->SetParameter(5,norm);
323 fLastFunc->SetParLimits(1,0.01,1);
324 fLastFunc->SetParLimits(2,0.01,100);
325 fLastFunc->SetParLimits(3,0.01,100);
326 fLastFunc->SetParLimits(4,0.01,100);
327 fLastFunc->SetParNames("mass","p0star","pt0","n","T","norm");
328 fLastFunc->SetNpx(5000);
329 fLastFunc->SetLineWidth(fLineWidth);
330 return fLastFunc;
331
4c0d7bc7 332 break;
333 case kOneOverMtdNdmt:
334 AliFatal("Not Implemented");
335 // return GetUA1dNdmt(mass,T,n,norm,name);
336 break;
337 default:
338 AliFatal("Not implemented");
339 }
340
341 return 0;
342}
343
344
345
346
347// ________________________________________________________________________
348
349// Backend (private functions and support functions for numerical integration)
350
fef6f2b3 351Double_t AliPWGFunc::StaticHistoFunc(const double * x, const double* p){
4c0d7bc7 352
353 // provides a function interpolating a histo with a spline;
354 // using double to store a pointer... This is a bad hack. To be replaced
355
356 double norm = p[0];
357
a4097895 358 TObject * h = (TObject*) Long64_t(p[1]);
4c0d7bc7 359
360// Int_t bin = h->FindBin(x[0]);
361// double value = h->GetBinContent(bin);
362
4c0d7bc7 363
58eeda84 364 // static TH1 * oldptr = 0;
365 // static TSpline3 * spl = 0;
366 // if (h!=oldptr) {
367 // FIXME: recheck static pointers
a4097895 368 TSpline3 * spl = 0;
369 if(h->InheritsFrom("TH1")) {
370 if ( ((TH1*)h)->FindBin(x[0]) > ((TH1*)h)->GetNbinsX()) return 0;
371 spl= new TSpline3((TH1*)h);
372 }
373 else if(h->InheritsFrom("TGraph")) spl= new TSpline3("fGraph",(TGraph*)h);
374 else {
fef6f2b3 375 Printf("AliPWGFunc::StaticHistoFunc: Unsupported type");
a4097895 376 return 0;
377 }
58eeda84 378 // }
4c0d7bc7 379 double value = spl->Eval(x[0]);
380 delete spl;
381
382 return value*norm;
383
384}
385
fef6f2b3 386Double_t AliPWGFunc::StaticUA1Func(const double * x, const double* p) {
4c0d7bc7 387
388
389 // "mass","p0star","pt0","n","T","norm"
390 Double_t mass = p[0];
391 Double_t p0star = p[1];
392 Double_t pt0 = p[2];
393 Double_t n = p[3];
8ffc8c1c 394 Double_t temp = p[4];
4c0d7bc7 395 Double_t norm = p[5];
396
397 Double_t xx = x[0];
398
fef6f2b3 399 static AliPWGFunc * self = new AliPWGFunc;
4c0d7bc7 400 static TF1 * fPLaw = self->GetPowerLawdNdptTimesPt(pt0, n, norm, "fLocalPLawUA1");
8ffc8c1c 401 static TF1 * fPMTExp = self->GetMTExpdNdptTimesPt (mass, temp, norm, "fLocalMTexpUA1");
4c0d7bc7 402
403 fPLaw->SetParameters(norm,pt0,n);
8ffc8c1c 404 fPMTExp->SetParameters(1,temp);
4c0d7bc7 405
406
407 Double_t normMT =fPMTExp->Eval(p0star) >0 ? fPLaw->Eval(p0star) / fPMTExp->Eval(p0star) * fPMTExp->GetParameter(0) : 1;
408 fPMTExp->SetParameter(0,normMT);
409
410
411 if (TMath::Abs(fPMTExp->Eval(p0star) - fPLaw->Eval(p0star)) > 0.0001 ) {
fef6f2b3 412 Printf("AliPWGFunc::StaticUA1Func - Wrong norm") ;
413 Printf(" p0* %f NMT: %f N: %f PL: %f MT: %f", p0star, normMT, norm, fPLaw->Eval(p0star), fPMTExp->Eval(p0star));
414 }
415
416 if (xx > p0star) return fPLaw->Eval(xx);
417 return fPMTExp->Eval(xx);
418
419
420
421}
422Double_t AliPWGFunc::StaticUA1FuncOneOverPt(const double * x, const double* p) {
423 // UA1 func
424
425 // "mass","p0star","pt0","n","T","norm"
426 Double_t mass = p[0];
427 Double_t p0star = p[1];
428 Double_t pt0 = p[2];
429 Double_t n = p[3];
430 Double_t temp = p[4];
431 Double_t norm = p[5];
432
433 Double_t xx = x[0];
434
435 static AliPWGFunc * self = new AliPWGFunc;
436 static TF1 * fPLaw = self->GetPowerLawdNdpt(pt0, n, norm, "fLocalPLawUA1");
437 static TF1 * fPMTExp = self->GetMTExpdNdpt (mass, temp, norm, "fLocalMTexpUA1");
438
439 fPLaw->SetParameters(norm,pt0,n);
440 fPMTExp->SetParameters(1,temp);
441
442
443 Double_t normMT =fPMTExp->Eval(p0star) >0 ? fPLaw->Eval(p0star) / fPMTExp->Eval(p0star) * fPMTExp->GetParameter(0) : 1;
444 fPMTExp->SetParameter(0,normMT);
445
446
447 if (TMath::Abs(fPMTExp->Eval(p0star) - fPLaw->Eval(p0star)) > 0.0001 ) {
448 Printf("AliPWGFunc::StaticUA1Func - Wrong norm") ;
4c0d7bc7 449 Printf(" p0* %f NMT: %f N: %f PL: %f MT: %f", p0star, normMT, norm, fPLaw->Eval(p0star), fPMTExp->Eval(p0star));
450 }
451
452 if (xx > p0star) return fPLaw->Eval(xx);
453 return fPMTExp->Eval(xx);
454
455
456
457}
458
459
fef6f2b3 460Double_t AliPWGFunc::IntegrandBG(const double * x, const double* p){
4c0d7bc7 461 // integrand for boltzman-gibbs blast wave
a4097895 462 // x[0] -> r (radius)
463 // p[0] -> mass
464 // p[1] -> pT (transverse momentum)
465 // p[2] -> beta_max (surface velocity)
466 // p[3] -> T (freezout temperature)
467 // p[4] -> n (velocity profile)
468
4c0d7bc7 469
470 double x0 = x[0];
471
a4097895 472 double mass = p[0];
473 double pT = p[1];
474 double beta_max = p[2];
475 double temp = p[3];
476 Double_t n = p[4];
477
478 // Keep beta within reasonable limits
479 Double_t beta = beta_max * TMath::Power(x0, n);
480 if (beta > 0.9999999999999999) beta = 0.9999999999999999;
481
4c0d7bc7 482 double mT = TMath::Sqrt(mass*mass+pT*pT);
483
a4097895 484 double rho0 = TMath::ATanH(beta);
8ffc8c1c 485 double arg00 = pT*TMath::SinH(rho0)/temp;
a4097895 486 if (arg00 > 700.) arg00 = 700.; // avoid FPE
8ffc8c1c 487 double arg01 = mT*TMath::CosH(rho0)/temp;
488 double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01);
4c0d7bc7 489
a4097895 490 // printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", x0, pT, beta_max, temp, n, mT, beta, rho0, arg00, arg01);
491
4c0d7bc7 492 return f0;
493}
494
495
496
fef6f2b3 497Double_t AliPWGFunc::StaticBGdNdPt(const double * x, const double* p) {
4c0d7bc7 498
499 // implementation of BGBW (1/pt dNdpt)
500
501 double pT = x[0];;
502
503
a4097895 504 double mass = p[0];
505 double beta = p[1];
8ffc8c1c 506 double temp = p[2];
a4097895 507 double n = p[3];
508 double norm = p[4];
4c0d7bc7 509
510 static TF1 * fIntBG = 0;
511 if(!fIntBG)
a4097895 512 fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 5);
4c0d7bc7 513
a4097895 514 fIntBG->SetParameters(mass, pT, beta, temp,n);
4c0d7bc7 515 double result = fIntBG->Integral(0,1);
a4097895 516 // printf ("[%4.4f], Int :%f\n", pT, result);
517 return result*norm;//*1e30;;
4c0d7bc7 518
519}
520
fef6f2b3 521Double_t AliPWGFunc::StaticBGdNdPtTimesPt(const double * x, const double* p) {
4c0d7bc7 522 // BGBW dNdpt implementation
523 return x[0]*StaticBGdNdPt(x,p);
524}
525
4fd89c20 526Double_t AliPWGFunc::StaticBGdNdMtTimesMt(const double * x, const double* p) {
527 // BGBW dNdpt implementation
528 // X0 is mt here
529 Double_t pt = TMath::Sqrt(x[0]*x[0]-p[0]*p[0]);
530 return pt*StaticBGdNdPt(&pt,p);
531}
532
4c0d7bc7 533
fef6f2b3 534TF1 * AliPWGFunc::GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t temp,
a4097895 535 Double_t n, Double_t norm, const char * name){
4c0d7bc7 536
537 // BGBW 1/pt dNdpt
538
a4097895 539 fLastFunc = new TF1 (name, StaticBGdNdPt, 0.0, 10, 5);
540 fLastFunc->SetParameters(mass,beta,temp,n,norm);
541 fLastFunc->FixParameter(0,mass);
542 fLastFunc->SetParNames("mass", "#beta", "T", "n", "norm");
4c0d7bc7 543 fLastFunc->SetLineWidth(fLineWidth);
544 return fLastFunc;
545
546}
547
548
549//_____________________________________________________________________
550// Tsallis
551
fef6f2b3 552Double_t AliPWGFunc::IntegrandTsallis(const double * x, const double* p){
4c0d7bc7 553
554 // integrand for numerical integration (tsallis)
555
556 Double_t r = x[0];
557 Double_t phi = x[1];
558 Double_t y = x[2];
559
560 Double_t mass = p[0];
561 Double_t pt = p[1];
562 Double_t beta = p[2];
8ffc8c1c 563 Double_t temp = p[3];
4c0d7bc7 564 Double_t q = p[4];
565
566 Double_t mt = TMath::Sqrt(mass*mass+pt*pt);
567
568 Double_t rho = TMath::ATanH(beta*r); // TODO: implement different velocity profiles
569
570 Double_t res = mt*
571 r*TMath::CosH(y) *TMath::Power( (
8ffc8c1c 572 1+(q-1)/temp * (
4c0d7bc7 573 mt*TMath::CosH(y)*TMath::CosH(rho) -
574 pt*TMath::SinH(rho)*TMath::Cos(phi)
575 )
576 ),
577 -1/(q-1)
578 );
579
580
581 return res;
582}
583
584
585
fef6f2b3 586Double_t AliPWGFunc::StaticTsallisdNdPt(const double * x, const double* p) {
4c0d7bc7 587
588 // tsallis BW implementation 1/pt dNdpt
589
590 double pT = x[0];;
591
592
593 double mass = p[0];
594 double beta = p[1];
8ffc8c1c 595 double temp = p[2];
4c0d7bc7 596 double q = p[3];
597
598 Double_t ymax = p[5];
599
600
601 static TF3 * fInt = 0;
602 if(!fInt){
603 fInt = new TF3 ("fIntTsa", IntegrandTsallis, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
604// fInt->SetNpx(10000);
605// fInt->SetNpy(10000);
606// fInt->SetNpz(10000);
607 }
608
8ffc8c1c 609 fInt->SetParameters(mass, pT, beta, temp, q);
4c0d7bc7 610 double result = fInt->Integral(0,1, -TMath::Pi(), TMath::Pi(), -ymax, ymax);
611 // double result = fInt->Integral(0,1, -2, 2, -ymax, ymax);
612
613 return result*p[4];//*1e30;;
614
615}
616
fef6f2b3 617Double_t AliPWGFunc::StaticTsallisdNdPtTimesPt(const double * x, const double* p) {
4c0d7bc7 618
619 // tsallis BW , implementatio of dNdpt
620 return x[0]*StaticTsallisdNdPt(x,p);
621
622}
623
fef6f2b3 624TF1 * AliPWGFunc::GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t temp, Double_t q,
4c0d7bc7 625 Double_t norm, Double_t ymax,const char * name){
626
627
628 // tsallis BW, 1/pt dNdpt
629
630 fLastFunc = new TF1 (name, StaticTsallisdNdPt, 0.0, 10, 6);
8ffc8c1c 631 fLastFunc->SetParameters(mass,beta,temp,q,norm,ymax);
4c0d7bc7 632 fLastFunc->SetParLimits(1,0.0,0.99);
633 fLastFunc->SetParLimits(2,0.01,0.99);
634 fLastFunc->SetParLimits(3,1.0001,1.9);
8ffc8c1c 635 fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax");
4c0d7bc7 636 fLastFunc->SetLineWidth(fLineWidth);
637 return fLastFunc;
638
639}
640
641// Times Pt funcs
642// Boltzmann-Gibbs Blast Wave
fef6f2b3 643TF1 * AliPWGFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t n,
4c0d7bc7 644 Double_t norm, const char * name){
645
646 // BGBW, dNdpt
647
a4097895 648 fLastFunc = new TF1 (name, StaticBGdNdPtTimesPt, 0.0, 10, 5);
649 fLastFunc->SetParameters(mass,beta,temp,n,norm);
650 fLastFunc->FixParameter(0,mass);
651 fLastFunc->SetParNames("mass", "#beta", "temp", "n", "norm");
4c0d7bc7 652 fLastFunc->SetLineWidth(fLineWidth);
653 return fLastFunc;
654
655
656}
657
4fd89c20 658// Boltzmann-Gibbs Blast Wave
659TF1 * AliPWGFunc::GetBGBWdNdptTimesMt(Double_t mass, Double_t beta, Double_t temp, Double_t n,
660 Double_t norm, const char * name){
661
662 // BGBW, dNdpt
663 // 1/Mt dN/dmt
664 fLastFunc = new TF1 (name, StaticBGdNdMtTimesMt, 0.0, 10, 5);
665 fLastFunc->SetParameters(mass,beta,temp,n,norm);
666 fLastFunc->FixParameter(0,mass);
667 fLastFunc->SetParNames("mass", "#beta", "temp", "n", "norm");
668 fLastFunc->SetLineWidth(fLineWidth);
669 return fLastFunc;
670
4c0d7bc7 671
4fd89c20 672}
4c0d7bc7 673
fef6f2b3 674TF1 * AliPWGFunc::GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t q,
4c0d7bc7 675 Double_t norm, Double_t ymax, const char * name){
676
677// Tsallis blast wave, dNdpt
678
679 fLastFunc = new TF1 (name, StaticTsallisdNdPtTimesPt, 0.0, 10, 6);
8ffc8c1c 680 fLastFunc->SetParameters(mass,beta,temp,q,norm,ymax);
681 fLastFunc->SetParNames("mass", "#beta", "temp", "q", "norm", "ymax");
4c0d7bc7 682 fLastFunc->SetLineWidth(fLineWidth);
683 return fLastFunc;
684
685
686}
687
688
689
fef6f2b3 690TF1 * AliPWGFunc::GetMTExpdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
4c0d7bc7 691
692 // Simple exponential in 1/mt*MT, as a function of dNdpt
693 char formula[500];
29b53d49 694 snprintf(formula,500,"[0]*x*exp(-sqrt(x**2+%f**2)/[1])", mass);
4c0d7bc7 695 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 696 fLastFunc->SetParameters(norm, temp);
4c0d7bc7 697 fLastFunc->SetParLimits(1, 0.01, 10);
698 fLastFunc->SetParNames("norm", "T");
699 fLastFunc->SetLineWidth(fLineWidth);
700 return fLastFunc;
701
702
703}
704
fef6f2b3 705TF1 * AliPWGFunc::GetBoseEinsteindNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
a4097895 706
707 // Bose einstein distribution as a function of dNdpt
708 char formula[500];
709 snprintf(formula,500,"[0]*x*1./(exp(sqrt(x**2+%f**2)/[1])-1)", mass);
710 fLastFunc=new TF1(name,formula,0,10);
711 fLastFunc->SetParameters(norm, temp);
712 fLastFunc->SetParLimits(1, 0.01, 10);
713 fLastFunc->SetParNames("norm", "T");
714 fLastFunc->SetLineWidth(fLineWidth);
715 return fLastFunc;
716
717
718}
719
fef6f2b3 720TF1 * AliPWGFunc::GetFermiDiracdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
a4097895 721
722 // Bose einstein distribution as a function of dNdpt
723 char formula[500];
724 snprintf(formula,500,"[0]*x*1./(exp(sqrt(x**2+%f**2)/[1])+1)", mass);
725 fLastFunc=new TF1(name,formula,0,10);
726 fLastFunc->SetParameters(norm, temp);
727 fLastFunc->SetParLimits(1, 0.01, 10);
728 fLastFunc->SetParNames("norm", "T");
729 fLastFunc->SetLineWidth(fLineWidth);
730 return fLastFunc;
731
732
733}
734
735
4c0d7bc7 736
fef6f2b3 737TF1 * AliPWGFunc::GetPTExpdNdptTimesPt(Double_t temp, Double_t norm, const char * name){
4c0d7bc7 738
739 // Simple exponential in 1/pt*dNdpT, as a function of dNdpt
740 char formula[500];
29b53d49 741 snprintf(formula,500,"[0]*x*exp(-x/[1])");
4c0d7bc7 742 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 743 fLastFunc->SetParameters(norm, temp);
4c0d7bc7 744 fLastFunc->SetParLimits(1, 0.01, 10);
745 fLastFunc->SetParNames("norm", "T");
746 fLastFunc->SetLineWidth(fLineWidth);
747 return fLastFunc;
748
749
750}
751
752
fef6f2b3 753TF1 * AliPWGFunc::GetBoltzmanndNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
4c0d7bc7 754 // Boltzmann (exp in 1/mt*dNdmT times mt) as a function of dNdpt
755 char formula[500];
29b53d49 756 snprintf(formula,500,"[0]*x*sqrt(x**2+%f**2)*exp(-sqrt(x**2+%f**2)/[1])", mass,mass);
4c0d7bc7 757 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 758 fLastFunc->SetParameters(norm, temp);
4c0d7bc7 759 fLastFunc->SetParLimits(1, 0.01, 10);
760 fLastFunc->SetParNames("norm", "T");
761 fLastFunc->SetLineWidth(fLineWidth);
762 return fLastFunc;
763
764
765}
766
767
768// Tsallis (no BW, a la CMS)
fef6f2b3 769// TF1 * AliPWGFunc::GetTsallisdNdptTimesPt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name){
4c0d7bc7 770
771// char formula[500];
772// // sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2)-%f)),(-1/([2]-1)))", mass, mass); //CMS
773// sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass); // STAR
774// //sprintf(formula,"[0]*x*sqrt(x**2+%f**2)*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass,mass); // STAR * mt
775// fLastFunc=new TF1(name,formula,0,10);
776// fLastFunc->SetParameters(norm, T, q);
777// fLastFunc->SetParLimits(1, 0.001, 10);
778// fLastFunc->SetParNames("norm", "T", "q");
779// fLastFunc->SetLineWidth(fLineWidth);
780// return fLastFunc;
781
782
783// }
784
785
fef6f2b3 786TF1 * AliPWGFunc::GetLevidNdptTimesPt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
4c0d7bc7 787
788 // Levi function, dNdpt
789 char formula[500];
790
29b53d49 791 snprintf(formula,500,"( x*[0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2]) )^(-[1])");
4c0d7bc7 792 // sprintf(formula,"( x*[0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (sqrt([3]*[3]+x*x))/([1]*[2]) )^(-[1])");
793 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 794 fLastFunc->SetParameters(norm, n, temp,mass);
4c0d7bc7 795 fLastFunc->SetParLimits(2, 0.01, 10);
796 fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
797 fLastFunc->FixParameter(3,mass);
798 fLastFunc->SetLineWidth(fLineWidth);
799 return fLastFunc;
800
801
802}
803
fef6f2b3 804TF1 * AliPWGFunc::GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name){
4c0d7bc7 805
806 // PowerLaw function, dNdpt
807 char formula[500];
808
29b53d49 809 snprintf(formula,500,"x*[0]*( 1 + x/[1] )^(-[2])");
4c0d7bc7 810 fLastFunc=new TF1(name,formula,0,10);
811 fLastFunc->SetParameters(norm, pt0, n);
812 fLastFunc->SetParLimits(1, 0.01, 10);
813 //fLastFunc->SetParLimits(2, 0.01, 50);
814 fLastFunc->SetParNames("norm", "pt0", "n");
815 fLastFunc->SetLineWidth(fLineWidth);
816 return fLastFunc;
817
818
819}
820
fef6f2b3 821TF1 * AliPWGFunc::GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name){
4c0d7bc7 822
823 // PowerLaw function, 1/pt dNdpt
824 char formula[500];
825
29b53d49 826 snprintf(formula,500," [0]*( 1 + x/[1] )^(-[2])");
4c0d7bc7 827 fLastFunc=new TF1(name,formula,0,10);
828 fLastFunc->SetParameters(norm, pt0, n);
829 // fLastFunc->SetParLimits(2, 0.01, 10);
830 fLastFunc->SetParNames("norm", "pt0", "n");
831 fLastFunc->SetLineWidth(fLineWidth);
832 return fLastFunc;
833
834
835}
836
837
fef6f2b3 838TF1 * AliPWGFunc::GetLevidNdpt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name){
4c0d7bc7 839
840 // Levi function, dNdpt
841 char formula[500];
842
29b53d49 843 snprintf(formula,500,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2]) )^(-[1])");
4c0d7bc7 844 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 845 fLastFunc->SetParameters(norm, n, temp,mass);
4c0d7bc7 846 fLastFunc->SetParLimits(2, 0.01, 10);
847 fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
848 fLastFunc->FixParameter(3,mass);
849 fLastFunc->SetLineWidth(fLineWidth);
850 return fLastFunc;
851
852
853}
854
fef6f2b3 855TF1 * AliPWGFunc::GetLevidNdmt(Double_t mass, Double_t temp, Double_t n, Double_t norm, const char * name, VarType_t var){
4c0d7bc7 856
0e019ec6 857 // Levi function, 1/mt dNdmt
4c0d7bc7 858 char formula[500];
0e019ec6 859 if (var == kOneOverMtdNdmt)
29b53d49 860 snprintf(formula,500,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x -[3])/([1]*[2]) )^(-[1])");
0e019ec6 861 else if (var == kdNdmt)
29b53d49 862 snprintf(formula,500,"( x*[0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x-[3])/([1]*[2]) )^(-[1])");
0e019ec6 863 if (var == kOneOverMtdNdmtMinusM)
29b53d49 864 snprintf(formula,500,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (x)/([1]*[2]) )^(-[1])");
0e019ec6 865
866 //sprintf(formula,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + x/([1]*[2]) )^(-[1])");
4c0d7bc7 867 // sprintf(formula,"[0] * ( 1 + x/([1]*[2]) )^(-[1])");
868 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 869 fLastFunc->SetParameters(norm, n, temp,mass);
4c0d7bc7 870 fLastFunc->SetParLimits(2, 0.01, 10);
871 fLastFunc->SetParNames("norm", "n", "T", "mass");
872 fLastFunc->FixParameter(3,mass);
873 fLastFunc->SetLineWidth(fLineWidth);
874 return fLastFunc;
875
876
877}
878
879
880
0e019ec6 881
4c0d7bc7 882// Test Function
fef6f2b3 883Double_t AliPWGFunc::IntegrandTest(const double * x, const double* p){
4c0d7bc7 884
885 // test function
886
887 Double_t y = x[0];
888
889 Double_t mass = p[0];
890 Double_t pt = p[1];
8ffc8c1c 891 Double_t temp = p[2];
4c0d7bc7 892
893 Double_t mt = TMath::Sqrt(mass*mass+pt*pt);
894
8ffc8c1c 895 return mt*TMath::CosH(y)*TMath::Exp(-mt*TMath::CosH(y)/temp);
4c0d7bc7 896
897}
898
fef6f2b3 899Double_t AliPWGFunc::StaticTest(const double * x, const double* p) {
4c0d7bc7 900
901 // test function
902
903 double pT = x[0];;
904
905
906 double mass = p[0];
8ffc8c1c 907 double temp = p[1];
4c0d7bc7 908 Double_t ymax = p[3];
909
910
911 static TF3 * fIntTest = 0;
912 if(!fIntTest){
913 fIntTest = new TF3 ("fIntTest", IntegrandTest, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
914 // fInt->SetNpx(10000);
915 }
916
8ffc8c1c 917 fIntTest->SetParameters(mass, pT, temp);
4c0d7bc7 918 double result = fIntTest->Integral(-ymax, ymax);
919
920 return result*p[2];//*1e30;;
921
922}
923
fef6f2b3 924TF1 * AliPWGFunc::GetTestFunc(Double_t mass, Double_t temp, Double_t norm, Double_t ymax, const char * name){
4c0d7bc7 925
926 // test function
927
928 fLastFunc = new TF1 (name, StaticTest, 0.0, 10, 4);
8ffc8c1c 929 fLastFunc->SetParameters(mass,temp,norm,ymax);
4c0d7bc7 930 fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
931 fLastFunc->SetLineWidth(fLineWidth);
932 return fLastFunc;
933
934}
935
936
937//___________________________________________________________
938
939
fef6f2b3 940TF1 * AliPWGFunc::GetMTExpdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
4c0d7bc7 941 // Simple exp in 1/mt dNdmt, as a function of dNdpt
942 // mt scaling
943 char formula[500];
29b53d49 944 snprintf(formula,500,"[0]*exp(-sqrt(x**2+%f**2)/[1])", mass);
4c0d7bc7 945 fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 946 fLastFunc->SetParameters(norm, temp);
4c0d7bc7 947 fLastFunc->SetParLimits(1, 0.01, 10);
948 fLastFunc->SetParNames("norm", "T");
949 fLastFunc->SetLineWidth(fLineWidth);
950 return fLastFunc;
951}
952
4fd89c20 953
954TF1 * AliPWGFunc::GetMTExpdNdmt(Double_t mass, Double_t temp, Double_t norm, const char * name, VarType_t var){
955
956 // Levi function, 1/mt dNdmt1/mt dNdmt,
957 char formula[500];
958 if (var == kOneOverMtdNdmt)
959 snprintf(formula,500,"[0] * exp (-x/[1]) + %f ", mass);
960 else if (var == kdNdmt)
961 snprintf(formula,500,"[0] * x * exp (-x/[1]) + %f ", mass);
962 if (var == kOneOverMtdNdmtMinusM)
38a77948 963 snprintf(formula,500,"[0] * exp (-x/[1])");
4fd89c20 964
965 //sprintf(formula,"( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + x/([1]*[2]) )^(-[1])");
966 // sprintf(formula,"[0] * ( 1 + x/([1]*[2]) )^(-[1])");
967 fLastFunc=new TF1(name,formula,0,10);
968 fLastFunc->SetParameters(norm, temp);
969 fLastFunc->SetParLimits(1, 0.01, 10);
970 fLastFunc->SetParNames("norm", "T");
971 fLastFunc->SetLineWidth(fLineWidth);
972 return fLastFunc;
973
974
975}
976
977
978
fef6f2b3 979TF1 * AliPWGFunc::GetBoseEinsteindNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
a4097895 980 // bose einstein
981 char formula[500];
982 snprintf(formula,500,"[0]*1./(exp(sqrt(x**2+%f**2)/[1])-1)", mass);
983 fLastFunc=new TF1(name,formula,0,10);
984 fLastFunc->SetParameters(norm, temp);
985 fLastFunc->SetParLimits(1, 0.01, 10);
986 fLastFunc->SetParNames("norm", "T");
987 fLastFunc->SetLineWidth(fLineWidth);
988 return fLastFunc;
989}
990
fef6f2b3 991TF1 * AliPWGFunc::GetFermiDiracdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
a4097895 992 // bose einstein
993 char formula[500];
994 snprintf(formula,500,"[0]*1./(exp(sqrt(x**2+%f**2)/[1])+1)", mass);
995 fLastFunc=new TF1(name,formula,0,10);
996 fLastFunc->SetParameters(norm, temp);
997 fLastFunc->SetParLimits(1, 0.01, 10);
998 fLastFunc->SetParNames("norm", "T");
999 fLastFunc->SetLineWidth(fLineWidth);
1000 return fLastFunc;
1001}
1002
4c0d7bc7 1003
1004// // Simple tsallis (a la CMS)
fef6f2b3 1005// TF1 * AliPWGFunc::GetTsallisdNdpt(Double_t mass, Double_t temp, Double_t q, Double_t norm, const char * name){
4c0d7bc7 1006
1007// char formula[500];
1008// sprintf(formula,"[0]*sqrt(x**2+%f**2)*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass,mass);
1009// fLastFunc=new TF1(name,formula,0,10);
8ffc8c1c 1010// fLastFunc->SetParameters(norm, temp, q);
4c0d7bc7 1011// fLastFunc->SetParLimits(1, 0.01, 10);
1012// fLastFunc->SetParNames("norm", "T", "q");
1013// fLastFunc->SetLineWidth(fLineWidth);
1014// return fLastFunc;
1015// }