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