]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenMUONlib.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
First cleanup of Loaders
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONlib.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
88cb7938 16/* $Id$ */
4c039060 17
53904666 18// Library class for particle pt and y distributions used for
19// muon spectrometer simulations.
20// To be used with AliGenParam.
21// The following particle typed can be simulated:
22// pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
23//
24// andreas.morsch@cern.ch
25//
26
65fb704d 27#include "TMath.h"
28#include "TRandom.h"
29
fe4da5cc 30#include "AliGenMUONlib.h"
5c3fd7ea 31
fe4da5cc 32ClassImp(AliGenMUONlib)
33//
34// Pions
75e0cc59 35Double_t AliGenMUONlib::PtPion(const Double_t *px, const Double_t* /*dummy*/)
fe4da5cc 36{
37//
38// PT-PARAMETERIZATION CDF, PRL 61(88) 1819
39// POWER LAW FOR PT > 500 MEV
40// MT SCALING BELOW (T=160 MEV)
41//
d90f80fd 42 const Double_t kp0 = 1.3;
43 const Double_t kxn = 8.28;
44 const Double_t kxlim=0.5;
45 const Double_t kt=0.160;
46 const Double_t kxmpi=0.139;
47 const Double_t kb=1.;
fe4da5cc 48 Double_t y, y1, xmpi2, ynorm, a;
49 Double_t x=*px;
50 //
d90f80fd 51 y1=TMath::Power(kp0/(kp0+kxlim),kxn);
52 xmpi2=kxmpi*kxmpi;
53 ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
fe4da5cc 54 a=ynorm/y1;
d90f80fd 55 if (x > kxlim)
56 y=a*TMath::Power(kp0/(kp0+x),kxn);
fe4da5cc 57 else
d90f80fd 58 y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
fe4da5cc 59 return y*x;
60}
753690b0 61//
62// y-distribution
63//
75e0cc59 64Double_t AliGenMUONlib::YPion( const Double_t *py, const Double_t */*dummy*/)
753690b0 65{
d90f80fd 66// Pion y
2280e6af 67 Double_t y=TMath::Abs(*py);
68/*
d90f80fd 69 const Double_t ka = 7000.;
70 const Double_t kdy = 4.;
d90f80fd 71 Double_t ex = y*y/(2*kdy*kdy);
72 return ka*TMath::Exp(-ex);
2280e6af 73*/
74 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
75
753690b0 76}
77// particle composition
78//
65fb704d 79Int_t AliGenMUONlib::IpPion(TRandom *ran)
753690b0 80{
d90f80fd 81// Pion composition
65fb704d 82 if (ran->Rndm() < 0.5) {
753690b0 83 return 211;
84 } else {
85 return -211;
86 }
87}
fe4da5cc 88
89//____________________________________________________________
90//
91// Mt-scaling
92
93Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
94{
95 // SCALING EN MASSE PAR RAPPORT A PTPI
96 // MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
d90f80fd 97 const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
fe4da5cc 98 // VALUE MESON/PI AT 5 GEV
d90f80fd 99 const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
fe4da5cc 100 np--;
d90f80fd 101 Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
102 Double_t fmax2=f5/kfmax[np];
fe4da5cc 103 // PIONS
104 Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
105 Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
d90f80fd 106 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
fe4da5cc 107 return fmtscal*ptpion;
108}
109//
753690b0 110// kaon
111//
112// pt-distribution
113//____________________________________________________________
75e0cc59 114Double_t AliGenMUONlib::PtKaon( const Double_t *px, const Double_t */*dummy*/)
753690b0 115{
d90f80fd 116// Kaon pT
753690b0 117 return PtScal(*px,2);
118}
119
120// y-distribution
fe4da5cc 121//____________________________________________________________
75e0cc59 122Double_t AliGenMUONlib::YKaon( const Double_t *py, const Double_t */*dummy*/)
fe4da5cc 123{
d90f80fd 124// Kaon y
2280e6af 125 Double_t y=TMath::Abs(*py);
126/*
d90f80fd 127 const Double_t ka = 1000.;
128 const Double_t kdy = 4.;
fe4da5cc 129 //
d90f80fd 130 Double_t ex = y*y/(2*kdy*kdy);
131 return ka*TMath::Exp(-ex);
2280e6af 132*/
133
134 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
753690b0 135}
136
137// particle composition
138//
65fb704d 139Int_t AliGenMUONlib::IpKaon(TRandom *ran)
753690b0 140{
d90f80fd 141// Kaon composition
65fb704d 142 if (ran->Rndm() < 0.5) {
753690b0 143 return 321;
144 } else {
145 return -321;
146 }
fe4da5cc 147}
753690b0 148
fe4da5cc 149// J/Psi
150//
151//
152// pt-distribution
153//____________________________________________________________
75e0cc59 154Double_t AliGenMUONlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/)
fe4da5cc 155{
d90f80fd 156// J/Psi pT
157 const Double_t kpt0 = 4.;
158 const Double_t kxn = 3.6;
fe4da5cc 159 Double_t x=*px;
160 //
d90f80fd 161 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
162 return x/TMath::Power(pass1,kxn);
fe4da5cc 163}
05932df6 164
75e0cc59 165Double_t AliGenMUONlib::PtJpsiCDFscaled( const Double_t *px, const Double_t */*dummy*/)
0e137c25 166{
167// J/Psi pT
9e9ae065 168//
169// PbPb 5.5 TeV
170// scaled from CDF data at 2 TeV
171// see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
172
173 const Double_t kpt0 = 5.100;
174 const Double_t kxn = 4.102;
175 Double_t x=*px;
176 //
177 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
178 return x/TMath::Power(pass1,kxn);
179}
180
75e0cc59 181Double_t AliGenMUONlib::PtJpsiCDFscaledold( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 182{
183// J/Psi pT
184//
185// PbPb 5.5 TeV
186// scaled from CDF "old" data at 2 TeV
187
48416d65 188 const Double_t kpt0 = 4.703;
189 const Double_t kxn = 3.826;
0e137c25 190 Double_t x=*px;
191 //
192 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
193 return x/TMath::Power(pass1,kxn);
194}
195
75e0cc59 196Double_t AliGenMUONlib::PtJpsiCDFscaledPP( const Double_t *px, const Double_t */*dummy*/)
9ff13849 197{
198// J/Psi pT
199//
200// pp 14 TeV
9ff13849 201// scaled from CDF data at 2 TeV
202
9e9ae065 203 const Double_t kpt0 = 5.630;
204 const Double_t kxn = 4.071;
205 Double_t x=*px;
206 //
207 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
208 return x/TMath::Power(pass1,kxn);
209}
210
75e0cc59 211Double_t AliGenMUONlib::PtJpsiCDFscaledPPold( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 212{
213// J/Psi pT
214//
215// pp 14 TeV
216// scaled from CDF "old" data at 2 TeV
217
9ff13849 218 const Double_t kpt0 = 5.355;
219 const Double_t kxn = 3.821;
220 Double_t x=*px;
221 //
222 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
223 return x/TMath::Power(pass1,kxn);
224}
225
75e0cc59 226Double_t AliGenMUONlib::PtJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 227{
228// J/Psi pT
229//
230// pp 10 TeV
231// scaled from CDF data at 2 TeV
232
233 const Double_t kpt0 = 5.334;
234 const Double_t kxn = 4.071;
235 Double_t x=*px;
236 //
237 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
238 return x/TMath::Power(pass1,kxn);
239}
240
75e0cc59 241Double_t AliGenMUONlib::PtJpsiFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
4ca8d070 242{
243 return 1.;
244}
245
75e0cc59 246Double_t AliGenMUONlib::PtJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
05932df6 247{
1af7144e 248// J/Psi pT spectrum
05932df6 249//
250// R. Vogt 2002
251// PbPb 5.5 TeV
252// MRST HO
253// mc = 1.4 GeV, pt-kick 1 GeV
254//
1af7144e 255 Float_t x = px[0];
256 Float_t c[8] = {
257 -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
258 -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
05932df6 259 };
1af7144e 260
3d905dd7 261 Double_t y;
36349df1 262 if (x < 10.) {
3d905dd7 263 Int_t j;
264 y = c[j = 7];
265 while (j > 0) y = y * x +c[--j];
266 y = x * TMath::Exp(y);
267 } else {
268 y = 0.;
269 }
1af7144e 270 return y;
05932df6 271}
17d28ba5 272
75e0cc59 273Double_t AliGenMUONlib::PtJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
17d28ba5 274{
275// J/Psi pT spectrum
276// B -> J/Psi X
277 Double_t x0 = 4.0384;
278 Double_t n = 3.0288;
279
280 Double_t x = px[0];
281 Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
282
283 return y;
284}
285
286
75e0cc59 287Double_t AliGenMUONlib::PtJpsiPP( const Double_t *px, const Double_t */*dummy*/)
bb6e81ac 288{
289// J/Psi pT spectrum
290//
291// R. Vogt 2002
292// pp 14 TeV
293// MRST HO
294// mc = 1.4 GeV, pt-kick 1 GeV
295//
296 Float_t x = px[0];
297 Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
298
299 Double_t y;
300 if (x < 10.) {
301 Int_t j;
302 y = c[j = 3];
303 while (j > 0) y = y * x +c[--j];
304 y = x * TMath::Exp(y);
305 } else {
306 y = 0.;
307 }
308 return y;
309}
310
fe4da5cc 311//
312// y-distribution
313//____________________________________________________________
75e0cc59 314Double_t AliGenMUONlib::YJpsi(const Double_t *py, const Double_t */*dummy*/)
fe4da5cc 315{
d90f80fd 316// J/psi y
317 const Double_t ky0 = 4.;
318 const Double_t kb=1.;
fe4da5cc 319 Double_t yj;
320 Double_t y=TMath::Abs(*py);
321 //
d90f80fd 322 if (y < ky0)
323 yj=kb;
fe4da5cc 324 else
d90f80fd 325 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
fe4da5cc 326 return yj;
327}
05932df6 328
75e0cc59 329Double_t AliGenMUONlib::YJpsiFlat( const Double_t */*py*/, const Double_t */*dummy*/ )
4ca8d070 330{
331 return 1.;
332}
333
05932df6 334
75e0cc59 335Double_t AliGenMUONlib::YJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
05932df6 336{
337
338//
339// J/Psi y
340//
341//
342// R. Vogt 2002
343// PbPb 5.5 TeV
344// MRST HO
345// mc = 1.4 GeV, pt-kick 1 GeV
346//
1af7144e 347 Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
348 Double_t x = TMath::Abs(px[0]);
349 Double_t y;
350
351 if (x < 4.) {
352 y = 31.754;
353 } else if (x < 6) {
354 Int_t j;
355 y = c[j = 4];
356 while (j > 0) y = y * x + c[--j];
357 } else {
358 y =0.;
359 }
360
361 return y;
05932df6 362}
363
75e0cc59 364Double_t AliGenMUONlib::YJpsiCDFscaled( const Double_t *px, const Double_t* dummy)
0e137c25 365{
366 // J/Psi y
367 return AliGenMUONlib::YJpsiPbPb(px, dummy);
368}
369
75e0cc59 370Double_t AliGenMUONlib::YJpsiCDFscaledold( const Double_t *px, const Double_t* dummy)
9e9ae065 371{
372 // J/Psi y
373 return AliGenMUONlib::YJpsiPbPb(px, dummy);
374}
375
75e0cc59 376Double_t AliGenMUONlib::YJpsiCDFscaledPP( const Double_t *px, const Double_t* dummy)
9ff13849 377{
378 // J/Psi y
379 return AliGenMUONlib::YJpsiPP(px, dummy);
380}
0e137c25 381
75e0cc59 382Double_t AliGenMUONlib::YJpsiCDFscaledPPold( const Double_t *px, const Double_t* dummy)
9e9ae065 383{
384 // J/Psi y
385 return AliGenMUONlib::YJpsiPP(px, dummy);
386}
387
75e0cc59 388Double_t AliGenMUONlib::YJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 389{
390
391//
392// J/Psi y
393//
394// pp 10 TeV
395// scaled from YJpsiPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD.
396// see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
397//
398
399 Double_t c[5] = {2.46681e+01, 8.91486e+01, -3.21227e+01, 3.63075e+00, -1.32047e-01};
400
401 Double_t x = TMath::Abs(px[0]);
402 Double_t y;
403
404 if (x < 3.2) {
405 y = 98.523 - 1.3664 * x * x;
406 } else if (x < 7.5) {
407 Int_t j;
408 y = c[j = 4];
409 while (j > 0) y = y * x + c[--j];
410 } else {
411 y =0.;
412 }
413
414 if(y<0) y=0;
415
416 return y;
417}
418
75e0cc59 419Double_t AliGenMUONlib::YJpsiPP( const Double_t *px, const Double_t */*dummy*/)
bb6e81ac 420{
421
422//
423// J/Psi y
424//
425//
426// R. Vogt 2002
427// pp 14 TeV
428// MRST HO
429// mc = 1.4 GeV, pt-kick 1 GeV
430//
431
432 Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
433 Double_t x = TMath::Abs(px[0]);
434 Double_t y;
435
436 if (x < 2.5) {
437 y = 96.455 - 0.8483 * x * x;
438 } else if (x < 7.9) {
439 Int_t j;
440 y = c[j = 4];
441 while (j > 0) y = y * x + c[--j];
442 } else {
443 y =0.;
444 }
445
446 return y;
447}
448
75e0cc59 449Double_t AliGenMUONlib::YJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
17d28ba5 450{
451
452//
453// J/Psi from B->J/Psi X
454//
455//
456
457
458 Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
459
460 Double_t x = TMath::Abs(px[0]);
461 Double_t y;
462
463 if (x > 6.) {
464 y = 0.;
465 } else {
466 Int_t j;
467 y = c[j = 6];
468 while (j > 0) y = y * x + c[--j];
469 }
470
471 return y;
472}
473
474
475
fe4da5cc 476// particle composition
477//
65fb704d 478Int_t AliGenMUONlib::IpJpsi(TRandom *)
fe4da5cc 479{
d90f80fd 480// J/Psi composition
88cb7938 481 return 443;
fe4da5cc 482}
88e5db43 483Int_t AliGenMUONlib::IpPsiP(TRandom *)
484{
485// Psi prime composition
486 return 100443;
487}
0ad09590 488Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
489{
490// J/Psi composition
491 Int_t ip;
492 Float_t r = gRandom->Rndm();
493 if (r < 0.98) {
494 ip = 443;
495 } else {
496 ip = 100443;
497 }
498 return ip;
499}
500
fe4da5cc 501
88e5db43 502
fe4da5cc 503// Upsilon
504//
505//
506// pt-distribution
507//____________________________________________________________
75e0cc59 508Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
fe4da5cc 509{
d90f80fd 510// Upsilon pT
511 const Double_t kpt0 = 5.3;
512 const Double_t kxn = 2.5;
fe4da5cc 513 Double_t x=*px;
514 //
d90f80fd 515 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
516 return x/TMath::Power(pass1,kxn);
fe4da5cc 517}
05932df6 518
75e0cc59 519Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
0e137c25 520{
521// Upsilon pT
48416d65 522 const Double_t kpt0 = 7.753;
523 const Double_t kxn = 3.042;
0e137c25 524 Double_t x=*px;
525 //
526 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
527 return x/TMath::Power(pass1,kxn);
528}
529
75e0cc59 530Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
9ff13849 531{
532// Upsilon pT
533//
534// pp 14 TeV
535//
536// scaled from CDF data at 2 TeV
537
538 const Double_t kpt0 = 8.610;
539 const Double_t kxn = 3.051;
540 Double_t x=*px;
541 //
542 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
543 return x/TMath::Power(pass1,kxn);
544}
545
75e0cc59 546Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
4ca8d070 547{
548 return 1.;
549}
550
75e0cc59 551Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
05932df6 552{
553
554//
555// Upsilon pT
556//
557//
558// R. Vogt 2002
559// PbPb 5.5 TeV
560// MRST HO
561// mc = 1.4 GeV, pt-kick 1 GeV
562//
1af7144e 563 Float_t x = px[0];
564 Double_t c[8] = {
565 -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
566 -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
567 };
3d905dd7 568 Double_t y;
569 if (x < 10.) {
570 Int_t j;
571 y = c[j = 7];
572 while (j > 0) y = y * x +c[--j];
573 y = x * TMath::Exp(y);
574 } else {
575 y = 0.;
576 }
1af7144e 577 return y;
05932df6 578}
579
75e0cc59 580Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
bb6e81ac 581{
582
583//
584// Upsilon pT
585//
586//
587// R. Vogt 2002
588// pp 14 TeV
589// MRST HO
590// mc = 1.4 GeV, pt-kick 1 GeV
591//
592 Float_t x = px[0];
593 Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
594 -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
595
596 Double_t y;
597 if (x < 10.) {
598 Int_t j;
599 y = c[j = 7];
600 while (j > 0) y = y * x +c[--j];
601 y = x * TMath::Exp(y);
602 } else {
603 y = 0.;
604 }
605 return y;
606}
607
fe4da5cc 608//
609// y-distribution
610//
611//____________________________________________________________
75e0cc59 612Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
fe4da5cc 613{
d90f80fd 614// Upsilon y
615 const Double_t ky0 = 3.;
616 const Double_t kb=1.;
fe4da5cc 617 Double_t yu;
618 Double_t y=TMath::Abs(*py);
619 //
d90f80fd 620 if (y < ky0)
621 yu=kb;
fe4da5cc 622 else
d90f80fd 623 yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
fe4da5cc 624 return yu;
625}
05932df6 626
627
75e0cc59 628Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
05932df6 629{
630
631//
632// Upsilon y
633//
634//
635// R. Vogt 2002
636// PbPb 5.5 TeV
637// MRST HO
638// mc = 1.4 GeV, pt-kick 1 GeV
639//
640
1af7144e 641 Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
642 -2.99753e-09, 1.28895e-05};
643
644 Double_t x = px[0];
645 if (TMath::Abs(x) > 5.55) return 0.;
646 Int_t j;
647 Double_t y = c[j = 6];
648 while (j > 0) y = y * x +c[--j];
649 return y;
05932df6 650}
651
75e0cc59 652Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
0e137c25 653{
654 // Upsilon y
655 return AliGenMUONlib::YUpsilonPbPb(px, dummy);
656
4ca8d070 657}
9ff13849 658
75e0cc59 659Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
9ff13849 660{
661 // Upsilon y
662 return AliGenMUONlib::YUpsilonPP(px, dummy);
663
664}
665
75e0cc59 666Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
4ca8d070 667{
668 // Upsilon y
669 return 1.;
670
0e137c25 671}
672
75e0cc59 673Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
bb6e81ac 674{
675
676//
677// Upsilon y
678//
679//
680// R. Vogt 2002
681// p p 14. TeV
682// MRST HO
683// mc = 1.4 GeV, pt-kick 1 GeV
684//
685 Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
686 -6.20539e-10, 1.29943e-05};
687
688 Double_t x = px[0];
689 if (TMath::Abs(x) > 6.2) return 0.;
690 Int_t j;
691 Double_t y = c[j = 6];
692 while (j > 0) y = y * x +c[--j];
693 return y;
694}
695
fe4da5cc 696// particle composition
697//
65fb704d 698Int_t AliGenMUONlib::IpUpsilon(TRandom *)
fe4da5cc 699{
d90f80fd 700// y composition
88cb7938 701 return 553;
fe4da5cc 702}
88e5db43 703Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
704{
705// y composition
706 return 100553;
707}
708Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
709{
710// y composition
711 return 200553;
712}
0ad09590 713Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
714{
715// y composition
716 Int_t ip;
717 Float_t r = gRandom->Rndm();
718
719 if (r < 0.712) {
720 ip = 553;
721 } else if (r < 0.896) {
722 ip = 100553;
723 } else {
724 ip = 200553;
725 }
726 return ip;
727}
fe4da5cc 728
88e5db43 729
fe4da5cc 730//
731// Phi
732//
733//
734// pt-distribution (by scaling of pion distribution)
735//____________________________________________________________
75e0cc59 736Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
fe4da5cc 737{
d90f80fd 738// Phi pT
fe4da5cc 739 return PtScal(*px,7);
740}
741// y-distribution
75e0cc59 742Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
fe4da5cc 743{
d90f80fd 744// Phi y
745 Double_t *dum=0;
746 return YJpsi(px,dum);
fe4da5cc 747}
748// particle composition
749//
65fb704d 750Int_t AliGenMUONlib::IpPhi(TRandom *)
fe4da5cc 751{
d90f80fd 752// Phi composition
89512a3b 753 return 333;
754}
755
756//
757// omega
758//
759//
760// pt-distribution (by scaling of pion distribution)
761//____________________________________________________________
75e0cc59 762Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
89512a3b 763{
764// Omega pT
765 return PtScal(*px,5);
766}
767// y-distribution
75e0cc59 768Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
89512a3b 769{
770// Omega y
771 Double_t *dum=0;
772 return YJpsi(px,dum);
773}
774// particle composition
775//
776Int_t AliGenMUONlib::IpOmega(TRandom *)
777{
778// Omega composition
779 return 223;
780}
781
782
783//
784// Eta
785//
786//
787// pt-distribution (by scaling of pion distribution)
788//____________________________________________________________
75e0cc59 789Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
89512a3b 790{
791// Eta pT
792 return PtScal(*px,3);
793}
794// y-distribution
75e0cc59 795Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
89512a3b 796{
797// Eta y
798 Double_t *dum=0;
799 return YJpsi(px,dum);
800}
801// particle composition
802//
803Int_t AliGenMUONlib::IpEta(TRandom *)
804{
805// Eta composition
806 return 221;
fe4da5cc 807}
808
809//
810// Charm
811//
812//
813// pt-distribution
814//____________________________________________________________
75e0cc59 815Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
fe4da5cc 816{
d90f80fd 817// Charm pT
bd0276a8 818 const Double_t kpt0 = 2.25;
819 const Double_t kxn = 3.17;
fe4da5cc 820 Double_t x=*px;
821 //
bd0276a8 822 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
823 return x/TMath::Power(pass1,kxn);
824}
825
75e0cc59 826Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
bd0276a8 827{
828// Charm pT
829 const Double_t kpt0 = 2.12;
830 const Double_t kxn = 2.78;
bd0276a8 831 Double_t x=*px;
832 //
833 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
d90f80fd 834 return x/TMath::Power(pass1,kxn);
fe4da5cc 835}
75e0cc59 836Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 837{
838// FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
839// PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
840// i=0,1,2; j=0,1,2; k=0,1,...,6
841// dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
842// http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
843// calculations for the following inputs:
844// Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
845// for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
846// for j=0,1 & 2 respectively;
847// factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
848// with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
849// for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
850// (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
851// June 2008, Smbat.Grigoryan@cern.ch
852
853// Charm pT
854// Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
855// for pp collisions at 14 TeV with one c-cbar pair per event.
856// Corresponding NLO total cross section is 5.68 mb
857
858
859 const Double_t kpt0 = 2.2930;
860 const Double_t kxn = 3.1196;
861 Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
862 Double_t x=*px;
863 //
864 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
865 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
866}
75e0cc59 867Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 868{
869// Charm pT
870// Corresponding NLO total cross section is 6.06 mb
871 const Double_t kpt0 = 2.8669;
872 const Double_t kxn = 3.1044;
873 Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
874 Double_t x=*px;
875 //
876 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
877 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
878}
75e0cc59 879Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 880{
881// Charm pT
882// Corresponding NLO total cross section is 6.06 mb
883 const Double_t kpt0 = 1.8361;
884 const Double_t kxn = 3.2966;
885 Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
886 Double_t x=*px;
887 //
888 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
889 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
890}
75e0cc59 891Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 892{
893// Charm pT
894// Corresponding NLO total cross section is 7.69 mb
895 const Double_t kpt0 = 2.1280;
896 const Double_t kxn = 3.1397;
897 Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
898 Double_t x=*px;
899 //
900 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
901 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
902}
75e0cc59 903Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 904{
905// Charm pT
906// Corresponding NLO total cross section is 4.81 mb
907 const Double_t kpt0 = 2.4579;
908 const Double_t kxn = 3.1095;
909 Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
910 Double_t x=*px;
911 //
912 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
913 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
914}
75e0cc59 915Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 916{
917// Charm pT
918// Corresponding NLO total cross section is 14.09 mb
919 const Double_t kpt0 = 2.1272;
920 const Double_t kxn = 3.1904;
921 Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
922 Double_t x=*px;
923 //
924 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
925 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
926}
75e0cc59 927Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 928{
929// Charm pT
930// Corresponding NLO total cross section is 1.52 mb
931 const Double_t kpt0 = 2.8159;
932 const Double_t kxn = 3.0857;
933 Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
934 Double_t x=*px;
935 //
936 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
937 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
938}
75e0cc59 939Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 940{
941// Charm pT
942// Corresponding NLO total cross section is 3.67 mb
943 const Double_t kpt0 = 2.7297;
944 const Double_t kxn = 3.3019;
945 Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
946 Double_t x=*px;
947 //
948 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
949 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
950}
75e0cc59 951Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 952{
953// Charm pT
954// Corresponding NLO total cross section is 3.38 mb
955 const Double_t kpt0 = 2.3894;
956 const Double_t kxn = 3.1075;
957 Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
958 Double_t x=*px;
959 //
960 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
961 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
962}
75e0cc59 963Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 964{
965// Charm pT
966// Corresponding NLO total cross section is 10.37 mb
967 const Double_t kpt0 = 2.0187;
968 const Double_t kxn = 3.3011;
969 Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
970 Double_t x=*px;
971 //
972 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
973 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
974}
75e0cc59 975Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 976{
977// Charm pT
978// Corresponding NLO total cross section is 7.22 mb
979 const Double_t kpt0 = 2.1089;
980 const Double_t kxn = 3.1848;
981 Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
982 Double_t x=*px;
983 //
984 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
985 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
986}
987
fe4da5cc 988// y-distribution
75e0cc59 989Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
fe4da5cc 990{
ec772ba2 991// Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
992// Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
993// shadowing + kt broadening
994
995 Double_t x=px[0];
996 Double_t c[2]={-2.42985e-03,-2.31001e-04};
997 Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
998 Double_t ycharm;
999
1000 if (TMath::Abs(x)>8) {
1001 ycharm=0.;
1002 }
1003 else {
1004 ycharm=TMath::Power(y,3);
1005 }
1006
1007 return ycharm;
fe4da5cc 1008}
75e0cc59 1009Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1010{
1011// FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1012// YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
1013// i=0,1,2; j=0,1,2; k=0,1,...,6
1014// dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
1015// http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
1016// calculations for the following inputs:
1017// Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
1018// for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
1019// for j=0,1 & 2 respectively;
1020// factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
1021// with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
1022// k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1023// (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1024// June 2008, Smbat.Grigoryan@cern.ch
1025
1026// Charm y
1027// Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
1028// for pp collisions at 14 TeV with one c-cbar pair per event.
1029// Corresponding NLO total cross section is 5.68 mb
1030
1031 Double_t x=px[0];
1032 Double_t c[2]={7.0909e-03,6.1967e-05};
1033 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1034 Double_t ycharm;
1035
1036 if (TMath::Abs(x)>9) {
1037 ycharm=0.;
1038 }
1039 else {
1040 ycharm=TMath::Power(y,3);
1041 }
1042
1043 return ycharm;
1044}
75e0cc59 1045Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1046{
1047// Charm y
1048// Corresponding NLO total cross section is 6.06 mb
1049 Double_t x=px[0];
1050 Double_t c[2]={6.9707e-03,6.0971e-05};
1051 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1052 Double_t ycharm;
1053
1054 if (TMath::Abs(x)>9) {
1055 ycharm=0.;
1056 }
1057 else {
1058 ycharm=TMath::Power(y,3);
1059 }
1060
1061 return ycharm;
1062}
75e0cc59 1063Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1064{
1065// Charm y
1066// Corresponding NLO total cross section is 6.06 mb
1067 Double_t x=px[0];
1068 Double_t c[2]={7.1687e-03,6.5303e-05};
1069 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1070 Double_t ycharm;
1071
1072 if (TMath::Abs(x)>9) {
1073 ycharm=0.;
1074 }
1075 else {
1076 ycharm=TMath::Power(y,3);
1077 }
1078
1079 return ycharm;
1080}
75e0cc59 1081Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1082{
1083// Charm y
1084// Corresponding NLO total cross section is 7.69 mb
1085 Double_t x=px[0];
1086 Double_t c[2]={5.9090e-03,7.1854e-05};
1087 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1088 Double_t ycharm;
1089
1090 if (TMath::Abs(x)>9) {
1091 ycharm=0.;
1092 }
1093 else {
1094 ycharm=TMath::Power(y,3);
1095 }
1096
1097 return ycharm;
1098}
75e0cc59 1099Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1100{
1101// Charm y
1102// Corresponding NLO total cross section is 4.81 mb
1103 Double_t x=px[0];
1104 Double_t c[2]={8.0882e-03,5.5872e-05};
1105 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1106 Double_t ycharm;
1107
1108 if (TMath::Abs(x)>9) {
1109 ycharm=0.;
1110 }
1111 else {
1112 ycharm=TMath::Power(y,3);
1113 }
1114
1115 return ycharm;
1116}
75e0cc59 1117Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1118{
1119// Charm y
1120// Corresponding NLO total cross section is 14.09 mb
1121 Double_t x=px[0];
1122 Double_t c[2]={7.2520e-03,6.2691e-05};
1123 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1124 Double_t ycharm;
1125
1126 if (TMath::Abs(x)>9) {
1127 ycharm=0.;
1128 }
1129 else {
1130 ycharm=TMath::Power(y,3);
1131 }
1132
1133 return ycharm;
1134}
75e0cc59 1135Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1136{
1137// Charm y
1138// Corresponding NLO total cross section is 1.52 mb
1139 Double_t x=px[0];
1140 Double_t c[2]={1.1040e-04,1.4498e-04};
1141 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1142 Double_t ycharm;
1143
1144 if (TMath::Abs(x)>9) {
1145 ycharm=0.;
1146 }
1147 else {
1148 ycharm=TMath::Power(y,3);
1149 }
1150
1151 return ycharm;
1152}
75e0cc59 1153Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1154{
1155// Charm y
1156// Corresponding NLO total cross section is 3.67 mb
1157 Double_t x=px[0];
1158 Double_t c[2]={-3.1328e-03,1.8270e-04};
1159 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1160 Double_t ycharm;
1161
1162 if (TMath::Abs(x)>9) {
1163 ycharm=0.;
1164 }
1165 else {
1166 ycharm=TMath::Power(y,3);
1167 }
1168
1169 return ycharm;
1170}
75e0cc59 1171Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1172{
1173// Charm y
1174// Corresponding NLO total cross section is 3.38 mb
1175 Double_t x=px[0];
1176 Double_t c[2]={7.0865e-03,6.2532e-05};
1177 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1178 Double_t ycharm;
1179
1180 if (TMath::Abs(x)>9) {
1181 ycharm=0.;
1182 }
1183 else {
1184 ycharm=TMath::Power(y,3);
1185 }
1186
1187 return ycharm;
1188}
75e0cc59 1189Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1190{
1191// Charm y
1192// Corresponding NLO total cross section is 10.37 mb
1193 Double_t x=px[0];
1194 Double_t c[2]={7.7070e-03,5.3533e-05};
1195 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1196 Double_t ycharm;
1197
1198 if (TMath::Abs(x)>9) {
1199 ycharm=0.;
1200 }
1201 else {
1202 ycharm=TMath::Power(y,3);
1203 }
1204
1205 return ycharm;
1206}
75e0cc59 1207Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1208{
1209// Charm y
1210// Corresponding NLO total cross section is 7.22 mb
1211 Double_t x=px[0];
1212 Double_t c[2]={7.9195e-03,5.3823e-05};
1213 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1214 Double_t ycharm;
1215
1216 if (TMath::Abs(x)>9) {
1217 ycharm=0.;
1218 }
1219 else {
1220 ycharm=TMath::Power(y,3);
1221 }
1222
1223 return ycharm;
1224}
fe4da5cc 1225
ec772ba2 1226
65fb704d 1227Int_t AliGenMUONlib::IpCharm(TRandom *ran)
fe4da5cc 1228{
d90f80fd 1229// Charm composition
65fb704d 1230 Float_t random;
fe4da5cc 1231 Int_t ip;
1232// 411,421,431,4122
65fb704d 1233 random = ran->Rndm();
ec772ba2 1234// Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
1235// >>>>> cf. tab 4 p 11
1236
1237 if (random < 0.30) {
1238 ip=421;
1239 } else if (random < 0.60) {
1240 ip=-421;
1241 } else if (random < 0.70) {
1242 ip=411;
1243 } else if (random < 0.80) {
1244 ip=-411;
1245 } else if (random < 0.86) {
1246 ip=431;
1247 } else if (random < 0.92) {
1248 ip=-431;
1249 } else if (random < 0.96) {
1250 ip=4122;
fe4da5cc 1251 } else {
ec772ba2 1252 ip=-4122;
fe4da5cc 1253 }
fe4da5cc 1254
1255 return ip;
1256}
1257
fe4da5cc 1258//
1259// Beauty
1260//
1261//
1262// pt-distribution
1263//____________________________________________________________
75e0cc59 1264Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
fe4da5cc 1265{
d90f80fd 1266// Beauty pT
bd0276a8 1267 const Double_t kpt0 = 6.53;
1268 const Double_t kxn = 3.59;
1269 Double_t x=*px;
1270 //
1271 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1272 return x/TMath::Power(pass1,kxn);
1273}
1274
75e0cc59 1275Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
bd0276a8 1276{
1277// Beauty pT
1278 const Double_t kpt0 = 6.14;
1279 const Double_t kxn = 2.93;
fe4da5cc 1280 Double_t x=*px;
1281 //
d90f80fd 1282 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1283 return x/TMath::Power(pass1,kxn);
fe4da5cc 1284}
75e0cc59 1285Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1286{
1287// FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1288// PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
1289// i=0,1,2; j=0,1,2; k=0,1,...,6
1290// dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
1291// http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
1292// calculations for the following inputs:
1293// Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
1294// for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
1295// for j=0,1 & 2 respectively;
1296// factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
1297// with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
1298// k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1299// (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1300// June 2008, Smbat.Grigoryan@cern.ch
1301
1302// Beauty pT
1303// Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
1304// for pp collisions at 14 TeV with one b-bbar pair per event.
1305// Corresponding NLO total cross section is 0.494 mb
1306
1307 const Double_t kpt0 = 8.0575;
1308 const Double_t kxn = 3.1921;
1309 Double_t x=*px;
1310 //
1311 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1312 return x/TMath::Power(pass1,kxn);
1313}
75e0cc59 1314Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1315{
1316// Beauty pT
1317// Corresponding NLO total cross section is 0.445 mb
1318 const Double_t kpt0 = 8.6239;
1319 const Double_t kxn = 3.2911;
1320 Double_t x=*px;
1321 //
1322 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1323 return x/TMath::Power(pass1,kxn);
1324}
75e0cc59 1325Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1326{
1327// Beauty pT
1328// Corresponding NLO total cross section is 0.445 mb
1329 const Double_t kpt0 = 7.3367;
1330 const Double_t kxn = 3.0692;
1331 Double_t x=*px;
1332 //
1333 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1334 return x/TMath::Power(pass1,kxn);
1335}
75e0cc59 1336Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1337{
1338// Beauty pT
1339// Corresponding NLO total cross section is 0.518 mb
1340 const Double_t kpt0 = 7.6409;
1341 const Double_t kxn = 3.1364;
1342 Double_t x=*px;
1343 //
1344 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1345 return x/TMath::Power(pass1,kxn);
1346}
75e0cc59 1347Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1348{
1349// Beauty pT
1350// Corresponding NLO total cross section is 0.384 mb
1351 const Double_t kpt0 = 8.4948;
1352 const Double_t kxn = 3.2546;
1353 Double_t x=*px;
1354 //
1355 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1356 return x/TMath::Power(pass1,kxn);
1357}
75e0cc59 1358Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1359{
1360// Beauty pT
1361// Corresponding NLO total cross section is 0.648 mb
1362 const Double_t kpt0 = 7.6631;
1363 const Double_t kxn = 3.1621;
1364 Double_t x=*px;
1365 //
1366 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1367 return x/TMath::Power(pass1,kxn);
1368}
75e0cc59 1369Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1370{
1371// Beauty pT
1372// Corresponding NLO total cross section is 0.294 mb
1373 const Double_t kpt0 = 8.7245;
1374 const Double_t kxn = 3.2213;
1375 Double_t x=*px;
1376 //
1377 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1378 return x/TMath::Power(pass1,kxn);
1379}
75e0cc59 1380Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1381{
1382// Beauty pT
1383// Corresponding NLO total cross section is 0.475 mb
1384 const Double_t kpt0 = 8.5296;
1385 const Double_t kxn = 3.2187;
1386 Double_t x=*px;
1387 //
1388 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1389 return x/TMath::Power(pass1,kxn);
1390}
75e0cc59 1391Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1392{
1393// Beauty pT
1394// Corresponding NLO total cross section is 0.324 mb
1395 const Double_t kpt0 = 7.9440;
1396 const Double_t kxn = 3.1614;
1397 Double_t x=*px;
1398 //
1399 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1400 return x/TMath::Power(pass1,kxn);
1401}
75e0cc59 1402Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1403{
1404// Beauty pT
1405// Corresponding NLO total cross section is 0.536 mb
1406 const Double_t kpt0 = 8.2408;
1407 const Double_t kxn = 3.3029;
1408 Double_t x=*px;
1409 //
1410 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1411 return x/TMath::Power(pass1,kxn);
1412}
75e0cc59 1413Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1414{
1415// Beauty pT
1416// Corresponding NLO total cross section is 0.420 mb
1417 const Double_t kpt0 = 7.8041;
1418 const Double_t kxn = 3.2094;
1419 Double_t x=*px;
1420 //
1421 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1422 return x/TMath::Power(pass1,kxn);
1423}
1424
fe4da5cc 1425// y-distribution
75e0cc59 1426Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
fe4da5cc 1427{
ec772ba2 1428// Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
1429// Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
1430// shadowing + kt broadening
1431
1432 Double_t x=px[0];
1433 Double_t c[2]={-1.27590e-02,-2.42731e-04};
1434 Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
1435 Double_t ybeauty;
1436
1437 if (TMath::Abs(x)>6) {
1438 ybeauty=0.;
1439 }
1440 else {
1441 ybeauty=TMath::Power(y,3);
1442 }
1443
1444 return ybeauty;
fe4da5cc 1445}
75e0cc59 1446Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1447{
1448// FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1449// YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
1450// i=0,1,2; j=0,1,2; k=0,1,...,6
1451// dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
1452// http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
1453// calculations for the following inputs:
1454// Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
1455// for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
1456// for j=0,1 & 2 respectively;
1457// factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
1458// with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
1459// for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1460// (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1461// June 2008, Smbat.Grigoryan@cern.ch
1462
1463// Beauty y
1464// Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
1465// for pp collisions at 14 TeV with one b-bbar pair per event.
1466// Corresponding NLO total cross section is 0.494 mb
1467
fe4da5cc 1468
9e9ae065 1469 Double_t x=px[0];
1470 Double_t c[2]={1.2350e-02,9.2667e-05};
1471 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1472 Double_t ybeauty;
1473
1474 if (TMath::Abs(x)>7.6) {
1475 ybeauty=0.;
1476 }
1477 else {
1478 ybeauty=TMath::Power(y,3);
1479 }
1480
1481 return ybeauty;
1482}
75e0cc59 1483Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1484{
1485// Beauty y
1486// Corresponding NLO total cross section is 0.445 mb
1487 Double_t x=px[0];
1488 Double_t c[2]={1.2292e-02,9.1847e-05};
1489 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1490 Double_t ybeauty;
1491
1492 if (TMath::Abs(x)>7.6) {
1493 ybeauty=0.;
1494 }
1495 else {
1496 ybeauty=TMath::Power(y,3);
1497 }
1498
1499 return ybeauty;
1500}
75e0cc59 1501Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1502{
1503// Beauty y
1504// Corresponding NLO total cross section is 0.445 mb
1505 Double_t x=px[0];
1506 Double_t c[2]={1.2436e-02,9.3709e-05};
1507 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1508 Double_t ybeauty;
1509
1510 if (TMath::Abs(x)>7.6) {
1511 ybeauty=0.;
1512 }
1513 else {
1514 ybeauty=TMath::Power(y,3);
1515 }
1516
1517 return ybeauty;
1518}
75e0cc59 1519Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1520{
1521// Beauty y
1522// Corresponding NLO total cross section is 0.518 mb
1523 Double_t x=px[0];
1524 Double_t c[2]={1.1714e-02,1.0068e-04};
1525 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1526 Double_t ybeauty;
1527
1528 if (TMath::Abs(x)>7.6) {
1529 ybeauty=0.;
1530 }
1531 else {
1532 ybeauty=TMath::Power(y,3);
1533 }
1534
1535 return ybeauty;
1536}
75e0cc59 1537Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1538{
1539// Beauty y
1540// Corresponding NLO total cross section is 0.384 mb
1541 Double_t x=px[0];
1542 Double_t c[2]={1.2944e-02,8.5500e-05};
1543 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1544 Double_t ybeauty;
1545
1546 if (TMath::Abs(x)>7.6) {
1547 ybeauty=0.;
1548 }
1549 else {
1550 ybeauty=TMath::Power(y,3);
1551 }
1552
1553 return ybeauty;
1554}
75e0cc59 1555Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1556{
1557// Beauty y
1558// Corresponding NLO total cross section is 0.648 mb
1559 Double_t x=px[0];
1560 Double_t c[2]={1.2455e-02,9.2713e-05};
1561 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1562 Double_t ybeauty;
1563
1564 if (TMath::Abs(x)>7.6) {
1565 ybeauty=0.;
1566 }
1567 else {
1568 ybeauty=TMath::Power(y,3);
1569 }
1570
1571 return ybeauty;
1572}
75e0cc59 1573Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1574{
1575// Beauty y
1576// Corresponding NLO total cross section is 0.294 mb
1577 Double_t x=px[0];
1578 Double_t c[2]={1.0897e-02,1.1878e-04};
1579 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1580 Double_t ybeauty;
1581
1582 if (TMath::Abs(x)>7.6) {
1583 ybeauty=0.;
1584 }
1585 else {
1586 ybeauty=TMath::Power(y,3);
1587 }
1588
1589 return ybeauty;
1590}
75e0cc59 1591Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1592{
1593// Beauty y
1594// Corresponding NLO total cross section is 0.475 mb
1595 Double_t x=px[0];
1596 Double_t c[2]={1.0912e-02,1.1858e-04};
1597 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1598 Double_t ybeauty;
1599
1600 if (TMath::Abs(x)>7.6) {
1601 ybeauty=0.;
1602 }
1603 else {
1604 ybeauty=TMath::Power(y,3);
1605 }
1606
1607 return ybeauty;
1608}
75e0cc59 1609Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1610{
1611// Beauty y
1612// Corresponding NLO total cross section is 0.324 mb
1613 Double_t x=px[0];
1614 Double_t c[2]={1.2378e-02,9.2490e-05};
1615 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1616 Double_t ybeauty;
1617
1618 if (TMath::Abs(x)>7.6) {
1619 ybeauty=0.;
1620 }
1621 else {
1622 ybeauty=TMath::Power(y,3);
1623 }
1624
1625 return ybeauty;
1626}
75e0cc59 1627Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1628{
1629// Beauty y
1630// Corresponding NLO total cross section is 0.536 mb
1631 Double_t x=px[0];
1632 Double_t c[2]={1.2886e-02,8.2912e-05};
1633 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1634 Double_t ybeauty;
1635
1636 if (TMath::Abs(x)>7.6) {
1637 ybeauty=0.;
1638 }
1639 else {
1640 ybeauty=TMath::Power(y,3);
1641 }
1642
1643 return ybeauty;
1644}
75e0cc59 1645Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
9e9ae065 1646{
1647// Beauty y
1648// Corresponding NLO total cross section is 0.420 mb
1649 Double_t x=px[0];
1650 Double_t c[2]={1.3106e-02,8.0115e-05};
1651 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1652 Double_t ybeauty;
1653
1654 if (TMath::Abs(x)>7.6) {
1655 ybeauty=0.;
1656 }
1657 else {
1658 ybeauty=TMath::Power(y,3);
1659 }
1660
1661 return ybeauty;
1662}
ec772ba2 1663
65fb704d 1664Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
fe4da5cc 1665{
d90f80fd 1666// Beauty Composition
65fb704d 1667 Float_t random;
fe4da5cc 1668 Int_t ip;
ec772ba2 1669 random = ran->Rndm();
1670
1671// Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
1672// >>>>> cf. tab 4 p 11
1673
1674 if (random < 0.20) {
1675 ip=511;
1676 } else if (random < 0.40) {
1677 ip=-511;
1678 } else if (random < 0.605) {
1679 ip=521;
1680 } else if (random < 0.81) {
1681 ip=-521;
1682 } else if (random < 0.87) {
1683 ip=531;
1684 } else if (random < 0.93) {
1685 ip=-531;
1686 } else if (random < 0.965) {
1687 ip=5122;
fe4da5cc 1688 } else {
ec772ba2 1689 ip=-5122;
fe4da5cc 1690 }
fe4da5cc 1691
ec772ba2 1692 return ip;
fe4da5cc 1693}
1694
ec772ba2 1695
75e0cc59 1696typedef Double_t (*GenFunc) (const Double_t*, const Double_t*);
53904666 1697GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname) const
fe4da5cc 1698{
d90f80fd 1699// Return pointer to pT parameterisation
05932df6 1700 TString sname = TString(tname);
fe4da5cc 1701 GenFunc func;
753690b0 1702 switch (param)
fe4da5cc 1703 {
34f60c01 1704 case kPhi:
fe4da5cc 1705 func=PtPhi;
1706 break;
89512a3b 1707 case kOmega:
1708 func=PtOmega;
1709 break;
1710 case kEta:
1711 func=PtEta;
1712 break;
0ad09590 1713 case kJpsiFamily:
88e5db43 1714 case kPsiP:
34f60c01 1715 case kJpsi:
bb6e81ac 1716 if (sname == "Vogt" || sname == "Vogt PbPb") {
05932df6 1717 func=PtJpsiPbPb;
bb6e81ac 1718 } else if (sname == "Vogt pp") {
1719 func=PtJpsiPP;
0e137c25 1720 } else if (sname == "CDF scaled") {
1721 func=PtJpsiCDFscaled;
9ff13849 1722 } else if (sname == "CDF pp") {
1723 func=PtJpsiCDFscaledPP;
9e9ae065 1724 } else if (sname == "CDF pp 10") {
1725 func=PtJpsiCDFscaledPP10;
1726 } else if (sname == "CDF scaled old") {
1727 func=PtJpsiCDFscaledold;
1728 } else if (sname == "CDF pp old") {
1729 func=PtJpsiCDFscaledPPold;
4ca8d070 1730 } else if (sname == "Flat") {
1731 func=PtJpsiFlat;
05932df6 1732 } else {
1733 func=PtJpsi;
1734 }
fe4da5cc 1735 break;
17d28ba5 1736 case kJpsiFromB:
1737 func = PtJpsiBPbPb;
1738 break;
0ad09590 1739 case kUpsilonFamily:
88e5db43 1740 case kUpsilonP:
1741 case kUpsilonPP:
34f60c01 1742 case kUpsilon:
bb6e81ac 1743 if (sname == "Vogt" || sname == "Vogt PbPb") {
05932df6 1744 func=PtUpsilonPbPb;
bb6e81ac 1745 } else if (sname == "Vogt pp") {
1746 func=PtUpsilonPP;
0e137c25 1747 } else if (sname == "CDF scaled") {
1748 func=PtUpsilonCDFscaled;
9ff13849 1749 } else if (sname == "CDF pp") {
1750 func=PtUpsilonCDFscaledPP;
4ca8d070 1751 } else if (sname == "Flat") {
1752 func=PtUpsilonFlat;
05932df6 1753 } else {
1754 func=PtUpsilon;
1755 }
0ad09590 1756 break;
34f60c01 1757 case kCharm:
9e9ae065 1758 if (sname == "F0M0S0 pp") {
1759 func=PtCharmF0M0S0PP;
1760 } else if (sname == "F1M0S0 pp") {
1761 func=PtCharmF1M0S0PP;
1762 } else if (sname == "F2M0S0 pp") {
1763 func=PtCharmF2M0S0PP;
1764 } else if (sname == "F0M1S0 pp") {
1765 func=PtCharmF0M1S0PP;
1766 } else if (sname == "F0M2S0 pp") {
1767 func=PtCharmF0M2S0PP;
1768 } else if (sname == "F0M0S1 pp") {
1769 func=PtCharmF0M0S1PP;
1770 } else if (sname == "F0M0S2 pp") {
1771 func=PtCharmF0M0S2PP;
1772 } else if (sname == "F0M0S3 pp") {
1773 func=PtCharmF0M0S3PP;
1774 } else if (sname == "F0M0S4 pp") {
1775 func=PtCharmF0M0S4PP;
1776 } else if (sname == "F0M0S5 pp") {
1777 func=PtCharmF0M0S5PP;
1778 } else if (sname == "F0M0S6 pp") {
1779 func=PtCharmF0M0S6PP;
1780 } else if (sname == "central") {
bd0276a8 1781 func=PtCharmCentral;
1782 } else {
1783 func=PtCharm;
1784 }
fe4da5cc 1785 break;
34f60c01 1786 case kBeauty:
9e9ae065 1787 if (sname == "F0M0S0 pp") {
1788 func=PtBeautyF0M0S0PP;
1789 } else if (sname == "F1M0S0 pp") {
1790 func=PtBeautyF1M0S0PP;
1791 } else if (sname == "F2M0S0 pp") {
1792 func=PtBeautyF2M0S0PP;
1793 } else if (sname == "F0M1S0 pp") {
1794 func=PtBeautyF0M1S0PP;
1795 } else if (sname == "F0M2S0 pp") {
1796 func=PtBeautyF0M2S0PP;
1797 } else if (sname == "F0M0S1 pp") {
1798 func=PtBeautyF0M0S1PP;
1799 } else if (sname == "F0M0S2 pp") {
1800 func=PtBeautyF0M0S2PP;
1801 } else if (sname == "F0M0S3 pp") {
1802 func=PtBeautyF0M0S3PP;
1803 } else if (sname == "F0M0S4 pp") {
1804 func=PtBeautyF0M0S4PP;
1805 } else if (sname == "F0M0S5 pp") {
1806 func=PtBeautyF0M0S5PP;
1807 } else if (sname == "F0M0S6 pp") {
1808 func=PtBeautyF0M0S6PP;
1809 } else if (sname == "central") {
bd0276a8 1810 func=PtBeautyCentral;
1811 } else {
1812 func=PtBeauty;
1813 }
fe4da5cc 1814 break;
34f60c01 1815 case kPion:
753690b0 1816 func=PtPion;
1817 break;
34f60c01 1818 case kKaon:
753690b0 1819 func=PtKaon;
1820 break;
6b450aa0 1821 case kChic0:
1822 func=PtChic0;
00935af2 1823 break;
6b450aa0 1824 case kChic1:
1825 func=PtChic1;
00935af2 1826 break;
6b450aa0 1827 case kChic2:
1828 func=PtChic2;
00935af2 1829 break;
6b450aa0 1830 case kChic:
1831 func=PtChic;
00935af2 1832 break;
119b35c7 1833 default:
1834 func=0;
1835 printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
fe4da5cc 1836 }
1837 return func;
1838}
1839
53904666 1840GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
fe4da5cc 1841{
ac3faee4 1842 //
1843 // Return pointer to y- parameterisation
1844 //
05932df6 1845 TString sname = TString(tname);
fe4da5cc 1846 GenFunc func;
753690b0 1847 switch (param)
fe4da5cc 1848 {
34f60c01 1849 case kPhi:
fe4da5cc 1850 func=YPhi;
1851 break;
89512a3b 1852 case kEta:
1853 func=YEta;
1854 break;
1855 case kOmega:
1856 func=YOmega;
1857 break;
0ad09590 1858 case kJpsiFamily:
88e5db43 1859 case kPsiP:
34f60c01 1860 case kJpsi:
bb6e81ac 1861 if (sname == "Vogt" || sname == "Vogt PbPb") {
05932df6 1862 func=YJpsiPbPb;
bb6e81ac 1863 } else if (sname == "Vogt pp"){
1864 func=YJpsiPP;
0e137c25 1865 } else if (sname == "CDF scaled") {
1866 func=YJpsiCDFscaled;
9ff13849 1867 } else if (sname == "CDF pp") {
1868 func=YJpsiCDFscaledPP;
9e9ae065 1869 } else if (sname == "CDF pp 10") {
1870 func=YJpsiCDFscaledPP10;
1871 } else if (sname == "CDF scaled old") {
1872 func=YJpsiCDFscaledold;
1873 } else if (sname == "CDF pp old") {
1874 func=YJpsiCDFscaledPPold;
4ca8d070 1875 } else if (sname == "Flat") {
1876 func=YJpsiFlat;
05932df6 1877 } else {
1878 func=YJpsi;
1879 }
17d28ba5 1880 break;
1881 case kJpsiFromB:
1882 func = YJpsiBPbPb;
fe4da5cc 1883 break;
0ad09590 1884 case kUpsilonFamily:
88e5db43 1885 case kUpsilonP:
1886 case kUpsilonPP:
34f60c01 1887 case kUpsilon:
bb6e81ac 1888 if (sname == "Vogt" || sname == "Vogt PbPb") {
05932df6 1889 func=YUpsilonPbPb;
bb6e81ac 1890 } else if (sname == "Vogt pp") {
1891 func = YUpsilonPP;
0e137c25 1892 } else if (sname == "CDF scaled") {
1893 func=YUpsilonCDFscaled;
9ff13849 1894 } else if (sname == "CDF pp") {
1895 func=YUpsilonCDFscaledPP;
4ca8d070 1896 } else if (sname == "Flat") {
1897 func=YUpsilonFlat;
05932df6 1898 } else {
1899 func=YUpsilon;
1900 }
fe4da5cc 1901 break;
34f60c01 1902 case kCharm:
9e9ae065 1903 if (sname == "F0M0S0 pp") {
1904 func=YCharmF0M0S0PP;
1905 } else if (sname == "F1M0S0 pp") {
1906 func=YCharmF1M0S0PP;
1907 } else if (sname == "F2M0S0 pp") {
1908 func=YCharmF2M0S0PP;
1909 } else if (sname == "F0M1S0 pp") {
1910 func=YCharmF0M1S0PP;
1911 } else if (sname == "F0M2S0 pp") {
1912 func=YCharmF0M2S0PP;
1913 } else if (sname == "F0M0S1 pp") {
1914 func=YCharmF0M0S1PP;
1915 } else if (sname == "F0M0S2 pp") {
1916 func=YCharmF0M0S2PP;
1917 } else if (sname == "F0M0S3 pp") {
1918 func=YCharmF0M0S3PP;
1919 } else if (sname == "F0M0S4 pp") {
1920 func=YCharmF0M0S4PP;
1921 } else if (sname == "F0M0S5 pp") {
1922 func=YCharmF0M0S5PP;
1923 } else if (sname == "F0M0S6 pp") {
1924 func=YCharmF0M0S6PP;
1925 } else {
1926 func=YCharm;
1927 }
fe4da5cc 1928 break;
34f60c01 1929 case kBeauty:
9e9ae065 1930 if (sname == "F0M0S0 pp") {
1931 func=YBeautyF0M0S0PP;
1932 } else if (sname == "F1M0S0 pp") {
1933 func=YBeautyF1M0S0PP;
1934 } else if (sname == "F2M0S0 pp") {
1935 func=YBeautyF2M0S0PP;
1936 } else if (sname == "F0M1S0 pp") {
1937 func=YBeautyF0M1S0PP;
1938 } else if (sname == "F0M2S0 pp") {
1939 func=YBeautyF0M2S0PP;
1940 } else if (sname == "F0M0S1 pp") {
1941 func=YBeautyF0M0S1PP;
1942 } else if (sname == "F0M0S2 pp") {
1943 func=YBeautyF0M0S2PP;
1944 } else if (sname == "F0M0S3 pp") {
1945 func=YBeautyF0M0S3PP;
1946 } else if (sname == "F0M0S4 pp") {
1947 func=YBeautyF0M0S4PP;
1948 } else if (sname == "F0M0S5 pp") {
1949 func=YBeautyF0M0S5PP;
1950 } else if (sname == "F0M0S6 pp") {
1951 func=YBeautyF0M0S6PP;
1952 } else {
1953 func=YBeauty;
1954 }
fe4da5cc 1955 break;
34f60c01 1956 case kPion:
753690b0 1957 func=YPion;
1958 break;
34f60c01 1959 case kKaon:
753690b0 1960 func=YKaon;
1961 break;
6b450aa0 1962 case kChic0:
1963 func=YChic0;
00935af2 1964 break;
6b450aa0 1965 case kChic1:
1966 func=YChic1;
00935af2 1967 break;
6b450aa0 1968 case kChic2:
1969 func=YChic2;
00935af2 1970 break;
6b450aa0 1971 case kChic:
1972 func=YChic;
00935af2 1973 break;
119b35c7 1974 default:
1975 func=0;
1976 printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
fe4da5cc 1977 }
1978 return func;
1979}
00935af2 1980
1981//
1982// Chi
1983//
1984//
1985// pt-distribution
1986//____________________________________________________________
75e0cc59 1987Double_t AliGenMUONlib::PtChic0( const Double_t *px, const Double_t */*dummy*/)
00935af2 1988{
1989// Chi_c1 pT
1990 const Double_t kpt0 = 4.;
1991 const Double_t kxn = 3.6;
1992 Double_t x=*px;
1993 //
1994 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1995 return x/TMath::Power(pass1,kxn);
1996}
75e0cc59 1997Double_t AliGenMUONlib::PtChic1( const Double_t *px, const Double_t */*dummy*/)
00935af2 1998{
1999// Chi_c1 pT
2000 const Double_t kpt0 = 4.;
2001 const Double_t kxn = 3.6;
2002 Double_t x=*px;
2003 //
2004 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2005 return x/TMath::Power(pass1,kxn);
2006}
75e0cc59 2007Double_t AliGenMUONlib::PtChic2( const Double_t *px, const Double_t */*dummy*/)
00935af2 2008{
2009// Chi_c2 pT
2010 const Double_t kpt0 = 4.;
2011 const Double_t kxn = 3.6;
2012 Double_t x=*px;
2013 //
2014 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2015 return x/TMath::Power(pass1,kxn);
2016}
75e0cc59 2017Double_t AliGenMUONlib::PtChic( const Double_t *px, const Double_t */*dummy*/)
00935af2 2018{
2019// Chi_c family pT
2020 const Double_t kpt0 = 4.;
2021 const Double_t kxn = 3.6;
2022 Double_t x=*px;
2023 //
2024 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2025 return x/TMath::Power(pass1,kxn);
2026}
2027
2028//
2029// y-distribution
2030//____________________________________________________________
75e0cc59 2031Double_t AliGenMUONlib::YChic0(const Double_t *py, const Double_t */*dummy*/)
00935af2 2032{
2033// Chi-1c y
2034 const Double_t ky0 = 4.;
2035 const Double_t kb=1.;
2036 Double_t yj;
2037 Double_t y=TMath::Abs(*py);
2038 //
2039 if (y < ky0)
2040 yj=kb;
2041 else
2042 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2043 return yj;
2044}
2045
75e0cc59 2046Double_t AliGenMUONlib::YChic1(const Double_t *py, const Double_t */*dummy*/)
00935af2 2047{
2048// Chi-1c y
2049 const Double_t ky0 = 4.;
2050 const Double_t kb=1.;
2051 Double_t yj;
2052 Double_t y=TMath::Abs(*py);
2053 //
2054 if (y < ky0)
2055 yj=kb;
2056 else
2057 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2058 return yj;
2059}
2060
75e0cc59 2061Double_t AliGenMUONlib::YChic2(const Double_t *py, const Double_t */*dummy*/)
00935af2 2062{
2063// Chi-2c y
2064 const Double_t ky0 = 4.;
2065 const Double_t kb=1.;
2066 Double_t yj;
2067 Double_t y=TMath::Abs(*py);
2068 //
2069 if (y < ky0)
2070 yj=kb;
2071 else
2072 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2073 return yj;
2074}
2075
75e0cc59 2076Double_t AliGenMUONlib::YChic(const Double_t *py, const Double_t */*dummy*/)
00935af2 2077{
2078// Chi_c family y
2079 const Double_t ky0 = 4.;
2080 const Double_t kb=1.;
2081 Double_t yj;
2082 Double_t y=TMath::Abs(*py);
2083 //
2084 if (y < ky0)
2085 yj=kb;
2086 else
2087 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2088 return yj;
2089}
2090
2091// particle composition
2092//
6b450aa0 2093Int_t AliGenMUONlib::IpChic0(TRandom *)
00935af2 2094{
2095// Chi composition
2096 return 10441;
2097}
2098//
6b450aa0 2099Int_t AliGenMUONlib::IpChic1(TRandom *)
00935af2 2100{
2101// Chi composition
2102 return 20443;
2103}
6b450aa0 2104Int_t AliGenMUONlib::IpChic2(TRandom *)
00935af2 2105{
2106// Chi_c2 prime composition
2107 return 445;
2108}
6b450aa0 2109Int_t AliGenMUONlib::IpChic(TRandom *)
00935af2 2110{
2111// Chi composition
2112 Int_t ip;
2113 Float_t r = gRandom->Rndm();
2114 if (r < 0.001) {
2115 ip = 10441;
2116 } else if( r < 0.377 ) {
2117 ip = 20443;
2118 } else {
2119 ip = 445;
2120 }
2121 return ip;
2122}
2123
2124
2125//_____________________________________________________________
2126
65fb704d 2127typedef Int_t (*GenFuncIp) (TRandom *);
198bb1c7 2128GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* /*tname*/) const
fe4da5cc 2129{
d90f80fd 2130// Return pointer to particle type parameterisation
fe4da5cc 2131 GenFuncIp func;
753690b0 2132 switch (param)
fe4da5cc 2133 {
34f60c01 2134 case kPhi:
fe4da5cc 2135 func=IpPhi;
2136 break;
89512a3b 2137 case kEta:
2138 func=IpEta;
2139 break;
2140 case kOmega:
2141 func=IpOmega;
2142 break;
0ad09590 2143 case kJpsiFamily:
2144 func=IpJpsiFamily;
2145 break;
88e5db43 2146 case kPsiP:
2147 func=IpPsiP;
2148 break;
34f60c01 2149 case kJpsi:
17d28ba5 2150 case kJpsiFromB:
fe4da5cc 2151 func=IpJpsi;
2152 break;
34f60c01 2153 case kUpsilon:
fe4da5cc 2154 func=IpUpsilon;
2155 break;
0ad09590 2156 case kUpsilonFamily:
2157 func=IpUpsilonFamily;
2158 break;
88e5db43 2159 case kUpsilonP:
2160 func=IpUpsilonP;
2161 break;
2162 case kUpsilonPP:
2163 func=IpUpsilonPP;
2164 break;
34f60c01 2165 case kCharm:
fe4da5cc 2166 func=IpCharm;
2167 break;
34f60c01 2168 case kBeauty:
fe4da5cc 2169 func=IpBeauty;
2170 break;
34f60c01 2171 case kPion:
753690b0 2172 func=IpPion;
2173 break;
34f60c01 2174 case kKaon:
753690b0 2175 func=IpKaon;
2176 break;
6b450aa0 2177 case kChic0:
2178 func=IpChic0;
00935af2 2179 break;
6b450aa0 2180 case kChic1:
2181 func=IpChic1;
00935af2 2182 break;
6b450aa0 2183 case kChic2:
2184 func=IpChic2;
00935af2 2185 break;
6b450aa0 2186 case kChic:
2187 func=IpChic;
00935af2 2188 break;
119b35c7 2189 default:
2190 func=0;
2191 printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
fe4da5cc 2192 }
2193 return func;
2194}
2195
2196
753690b0 2197
05932df6 2198Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
2199 Float_t dx,
2200 Int_t n, Int_t no)
2201{
2202//
2203// Neville's alorithm for interpolation
2204//
2205// x: x-value
2206// y: Input array
2207// x0: minimum x
2208// dx: step size
2209// n: number of data points
2210// no: order of polynom
2211//
2212 Float_t* c = new Float_t[n];
2213 Float_t* d = new Float_t[n];
2214 Int_t m, i;
2215 for (i = 0; i < n; i++) {
2216 c[i] = y[i];
2217 d[i] = y[i];
2218 }
2219
2220 Int_t ns = int((x - x0)/dx);
2221
2222 Float_t y1 = y[ns];
2223 ns--;
2224 for (m = 0; m < no; m++) {
2225 for (i = 0; i < n-m; i++) {
2226 Float_t ho = x0 + Float_t(i) * dx - x;
2227 Float_t hp = x0 + Float_t(i+m+1) * dx - x;
2228 Float_t w = c[i+1] - d[i];
2229 Float_t den = ho-hp;
2230 den = w/den;
2231 d[i] = hp * den;
2232 c[i] = ho * den;
2233 }
2234 Float_t dy;
2235
2236 if (2*ns < (n-m-1)) {
2237 dy = c[ns+1];
2238 } else {
2239 dy = d[ns--];
2240 }
2241 y1 += dy;}
2242 delete[] c;
2243 delete[] d;
2244
2245 return y1;
2246}
2247
753690b0 2248