Small fix Sarah
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONlib.cxx
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
16 /* $Id$ */
17
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
27 #include "TMath.h"
28 #include "TRandom.h"
29 #include "TDatabasePDG.h"
30
31 #include "AliGenMUONlib.h"
32
33 ClassImp(AliGenMUONlib)
34 //
35 //  Pions
36 Double_t AliGenMUONlib::PtPion(const Double_t *px, const Double_t* /*dummy*/)
37 {
38 //
39 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
40 //     POWER LAW FOR PT > 500 MEV
41 //     MT SCALING BELOW (T=160 MEV)
42 //
43   const Double_t kp0 = 1.3;
44   const Double_t kxn = 8.28;
45   const Double_t kxlim=0.5;
46   const Double_t kt=0.160;
47   const Double_t kxmpi=0.139;
48   const Double_t kb=1.;
49   Double_t y, y1, xmpi2, ynorm, a;
50   Double_t x=*px;
51   //
52   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
53   xmpi2=kxmpi*kxmpi;
54   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
55   a=ynorm/y1;
56   if (x > kxlim)
57     y=a*TMath::Power(kp0/(kp0+x),kxn);
58   else
59     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
60   return y*x;
61 }
62 //
63 // y-distribution
64 //
65 Double_t AliGenMUONlib::YPion( const Double_t *py, const Double_t */*dummy*/)
66 {
67 // Pion y
68   Double_t y=TMath::Abs(*py);
69 /*
70   const Double_t ka    = 7000.;
71   const Double_t kdy   = 4.;
72   Double_t ex = y*y/(2*kdy*kdy);
73   return ka*TMath::Exp(-ex);
74 */
75   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
76   
77 }
78 //                 particle composition
79 //
80 Int_t AliGenMUONlib::IpPion(TRandom *ran)
81 {
82 // Pion composition 
83     if (ran->Rndm() < 0.5) {
84         return  211;
85     } else {
86         return -211;
87     }
88 }
89
90 //____________________________________________________________
91 //
92 // Mt-scaling
93
94 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
95 {
96   //    SCALING EN MASSE PAR RAPPORT A PTPI
97   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
98   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
99   //     VALUE MESON/PI AT 5 GEV
100   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
101   np--;
102   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
103   Double_t fmax2=f5/kfmax[np];
104   // PIONS
105   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
106   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
107                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
108   return fmtscal*ptpion;
109 }
110 //
111 // kaon
112 //
113 //                pt-distribution
114 //____________________________________________________________
115 Double_t AliGenMUONlib::PtKaon( const Double_t *px, const Double_t */*dummy*/)
116 {
117 // Kaon pT
118   return PtScal(*px,2);
119 }
120
121 // y-distribution
122 //____________________________________________________________
123 Double_t AliGenMUONlib::YKaon( const Double_t *py, const Double_t */*dummy*/)
124 {
125 // Kaon y
126   Double_t y=TMath::Abs(*py);
127 /*
128   const Double_t ka    = 1000.;
129   const Double_t kdy   = 4.;
130   //
131   Double_t ex = y*y/(2*kdy*kdy);
132   return ka*TMath::Exp(-ex);
133 */
134
135   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
136 }
137
138 //                 particle composition
139 //
140 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
141 {
142 // Kaon composition
143     if (ran->Rndm() < 0.5) {
144         return  321;
145     } else {
146         return -321;
147     }
148 }
149
150 //                    J/Psi 
151 //
152 //
153 //                pt-distribution
154 //____________________________________________________________
155 Double_t AliGenMUONlib::PtJpsiPPdummy(Double_t x, Double_t energy)
156 {
157 // J/Psi pT
158 // pp
159 // from the fit of RHIC, CDF & LHC data, see arXiv:1103.2394
160 //
161   const Double_t kpt0 = 1.04*TMath::Power(energy,0.101);
162   const Double_t kxn  = 3.9;
163   //
164   Double_t pass1 = 1.+0.363*(x/kpt0)*(x/kpt0);
165   return x/TMath::Power(pass1,kxn);
166 }
167
168 Double_t AliGenMUONlib::PtJpsiPP7000(const Double_t *px, const Double_t */*dummy*/)
169 {
170 // J/Psi pT
171 // pp 7 TeV
172 //
173   return PtJpsiPPdummy(*px,7000);
174 }
175
176 Double_t AliGenMUONlib::PtJpsiPP2760(const Double_t *px, const Double_t */*dummy*/)
177 {
178 // J/Psi pT
179 // pp 2.76 TeV
180 //
181   return PtJpsiPPdummy(*px,2760);
182 }
183
184 Double_t AliGenMUONlib::PtJpsiPP4400(const Double_t *px, const Double_t */*dummy*/)
185 {
186 // J/Psi pT
187 // pp 4.4 TeV
188 //
189   return PtJpsiPPdummy(*px,4400);
190 }
191
192 Double_t AliGenMUONlib::PtJpsiPP5030(const Double_t *px, const Double_t */*dummy*/)
193 {
194 // J/Psi pT
195 // pp 5.03 TeV
196 //
197   return PtJpsiPPdummy(*px,5030);
198 }
199
200 Double_t AliGenMUONlib::PtJpsiPP8800(const Double_t *px, const Double_t */*dummy*/)
201 {
202 // J/Psi pT
203 // pp 8.8 TeV
204 //
205   return PtJpsiPPdummy(*px,8800);
206 }
207
208 Double_t AliGenMUONlib::PtJpsiPbPb2760ShFdummy(Double_t x, Int_t n)
209 {
210 // J/Psi shadowing factor vs pT for PbPb min. bias and 11 centr. bins (in 2.5<y<4)
211 //
212 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.66 in 4pi
213 // S.Grigoryan, details presented at the PWG3-Muon meeting (05.10.2011)
214 // https://indico.cern.ch/conferenceDisplay.py?confId=157367
215 //
216   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
217                            0.428, 0.317, 0.231, 0.156};
218   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
219                            0.106, 0.041, 0.013, 0.002};
220   const Double_t c1[7] = {1.6077e+00, 7.6300e-02,-7.1880e-03, 3.4067e-04,
221                           -9.2776e-06,1.5138e-07, 1.4652e-09}; 
222   const Double_t c2[7] = {6.2047e-01, 5.7653e-02,-4.1414e-03, 1.0301e-04, 
223                           9.6205e-07,-7.4098e-08, 5.0946e-09}; 
224   Double_t y1, y2;
225   Int_t j;
226   y1 = c1[j = 6]; y2 = c2[6];
227   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
228   
229   y1 /= 1.+c1[6]*TMath::Power(x,6);
230   y2 /= 1.+c2[6]*TMath::Power(x,6);
231   //  
232   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
233   if(y1<0) y1=0;
234   return y1;
235 }
236
237 Double_t AliGenMUONlib::PtJpsiPbPb2760(const Double_t *px, const Double_t *dummy)
238 {
239 // J/Psi pT
240 // PbPb 2.76 TeV, minimum bias 0-100 %
241 //
242   return PtJpsiPbPb2760ShFdummy(*px, 0) * PtJpsiPP2760(px, dummy);
243 }
244
245 Double_t AliGenMUONlib::PtJpsiPbPb2760c1(const Double_t *px, const Double_t *dummy)
246 {
247 // J/Psi pT
248 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
249 //
250   return PtJpsiPbPb2760ShFdummy(*px, 1) * PtJpsiPP2760(px, dummy);
251 }
252
253 Double_t AliGenMUONlib::PtJpsiPbPb2760c2(const Double_t *px, const Double_t *dummy)
254 {
255 // J/Psi pT
256 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
257 //
258   return PtJpsiPbPb2760ShFdummy(*px, 2) * PtJpsiPP2760(px, dummy);
259 }
260
261 Double_t AliGenMUONlib::PtJpsiPbPb2760c3(const Double_t *px, const Double_t *dummy)
262 {
263 // J/Psi pT
264 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
265 //
266   return PtJpsiPbPb2760ShFdummy(*px, 3) * PtJpsiPP2760(px, dummy);
267 }
268
269 Double_t AliGenMUONlib::PtJpsiPbPb2760c4(const Double_t *px, const Double_t *dummy)
270 {
271 // J/Psi pT
272 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
273 //
274   return PtJpsiPbPb2760ShFdummy(*px, 4) * PtJpsiPP2760(px, dummy);
275 }
276
277 Double_t AliGenMUONlib::PtJpsiPbPb2760c5(const Double_t *px, const Double_t *dummy)
278 {
279 // J/Psi pT
280 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
281 //
282   return PtJpsiPbPb2760ShFdummy(*px, 5) * PtJpsiPP2760(px, dummy);
283 }
284
285 Double_t AliGenMUONlib::PtJpsiPbPb2760c6(const Double_t *px, const Double_t *dummy)
286 {
287 // J/Psi pT
288 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
289 //
290   return PtJpsiPbPb2760ShFdummy(*px, 6) * PtJpsiPP2760(px, dummy);
291 }
292
293 Double_t AliGenMUONlib::PtJpsiPbPb2760c7(const Double_t *px, const Double_t *dummy)
294 {
295 // J/Psi pT
296 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
297 //
298   return PtJpsiPbPb2760ShFdummy(*px, 7) * PtJpsiPP2760(px, dummy);
299 }
300
301 Double_t AliGenMUONlib::PtJpsiPbPb2760c8(const Double_t *px, const Double_t *dummy)
302 {
303 // J/Psi pT
304 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
305 //
306   return PtJpsiPbPb2760ShFdummy(*px, 8) * PtJpsiPP2760(px, dummy);
307 }
308
309 Double_t AliGenMUONlib::PtJpsiPbPb2760c9(const Double_t *px, const Double_t *dummy)
310 {
311 // J/Psi pT
312 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
313 //
314   return PtJpsiPbPb2760ShFdummy(*px, 9) * PtJpsiPP2760(px, dummy);
315 }
316
317 Double_t AliGenMUONlib::PtJpsiPbPb2760c10(const Double_t *px, const Double_t *dummy)
318 {
319 // J/Psi pT
320 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
321 //
322   return PtJpsiPbPb2760ShFdummy(*px, 10) * PtJpsiPP2760(px, dummy);
323 }
324
325 Double_t AliGenMUONlib::PtJpsiPbPb2760c11(const Double_t *px, const Double_t *dummy)
326 {
327 // J/Psi pT
328 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
329 //
330   return PtJpsiPbPb2760ShFdummy(*px, 11) * PtJpsiPP2760(px, dummy);
331 }
332
333 Double_t AliGenMUONlib::PtJpsiPPb5030ShFdummy(Double_t x, Int_t n)
334 {
335 // J/Psi shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
336 //
337 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.81 in 4pi
338 //
339   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
340   const Double_t c[7] = {6.675e-01, 1.808e-02, 2.721e-03,-7.793e-04, 7.504e-05,-3.884e-06, 5.759e-07}; 
341   Double_t y;
342   Int_t j;
343   y = c[j = 6];
344   while (j > 0) y  = y * x + c[--j];
345   y /= 1 + c[6]*TMath::Power(x,6);
346   //  
347   return 1 + (y-1)*f[n];
348 }
349
350 Double_t AliGenMUONlib::PtJpsiPPb5030(const Double_t *px, const Double_t *dummy)
351 {
352 // J/Psi pT
353 // pPb 5.03 TeV, minimum bias 0-100 %
354 //
355   return PtJpsiPPb5030ShFdummy(*px, 0) * PtJpsiPP5030(px, dummy);
356 }
357
358 Double_t AliGenMUONlib::PtJpsiPPb5030c1(const Double_t *px, const Double_t *dummy)
359 {
360 // J/Psi pT
361 // pPb 5.03 TeV, 1st centrality bin 0-20 %
362 //
363   return PtJpsiPPb5030ShFdummy(*px, 1) * PtJpsiPP5030(px, dummy);
364 }
365
366 Double_t AliGenMUONlib::PtJpsiPPb5030c2(const Double_t *px, const Double_t *dummy)
367 {
368 // J/Psi pT
369 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
370 //
371   return PtJpsiPPb5030ShFdummy(*px, 2) * PtJpsiPP5030(px, dummy);
372 }
373
374 Double_t AliGenMUONlib::PtJpsiPPb5030c3(const Double_t *px, const Double_t *dummy)
375 {
376 // J/Psi pT
377 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
378 //
379   return PtJpsiPPb5030ShFdummy(*px, 3) * PtJpsiPP5030(px, dummy);
380 }
381
382 Double_t AliGenMUONlib::PtJpsiPPb5030c4(const Double_t *px, const Double_t *dummy)
383 {
384 // J/Psi pT
385 // pPb 5.03 TeV, 4th centrality bin 60-100 %
386 //
387   return PtJpsiPPb5030ShFdummy(*px, 4) * PtJpsiPP5030(px, dummy);
388 }
389
390 Double_t AliGenMUONlib::PtJpsiPbP5030ShFdummy(Double_t x, Int_t n)
391 {
392 // J/Psi shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
393 //
394 // Pbp 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.81 in 4pi
395 //
396   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
397   const Double_t c[7] = {8.966e-01, 3.498e-02, 6.637e-03,-1.765e-03, 1.240e-04,-2.086e-06, 4.062e-07};
398   Double_t y;
399   Int_t j;
400   y = c[j = 6];
401   while (j > 0) y  = y * x + c[--j];
402   y /= 1 + c[6]*TMath::Power(x,6);
403   //  
404   return 1 + (y-1)*f[n];
405 }
406
407 Double_t AliGenMUONlib::PtJpsiPbP5030(const Double_t *px, const Double_t *dummy)
408 {
409 // J/Psi pT
410 // Pbp 5.03 TeV, minimum bias 0-100 %
411 //
412   return PtJpsiPbP5030ShFdummy(*px, 0) * PtJpsiPP5030(px, dummy);
413 }
414
415 Double_t AliGenMUONlib::PtJpsiPbP5030c1(const Double_t *px, const Double_t *dummy)
416 {
417 // J/Psi pT
418 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
419 //
420   return PtJpsiPbP5030ShFdummy(*px, 1) * PtJpsiPP5030(px, dummy);
421 }
422
423 Double_t AliGenMUONlib::PtJpsiPbP5030c2(const Double_t *px, const Double_t *dummy)
424 {
425 // J/Psi pT
426 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
427 //
428   return PtJpsiPbP5030ShFdummy(*px, 2) * PtJpsiPP5030(px, dummy);
429 }
430
431 Double_t AliGenMUONlib::PtJpsiPbP5030c3(const Double_t *px, const Double_t *dummy)
432 {
433 // J/Psi pT
434 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
435 //
436   return PtJpsiPbP5030ShFdummy(*px, 3) * PtJpsiPP5030(px, dummy);
437 }
438
439 Double_t AliGenMUONlib::PtJpsiPbP5030c4(const Double_t *px, const Double_t *dummy)
440 {
441 // J/Psi pT
442 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
443 //
444   return PtJpsiPbP5030ShFdummy(*px, 4) * PtJpsiPP5030(px, dummy);
445 }
446
447 Double_t AliGenMUONlib::PtJpsiPPb8800ShFdummy(Double_t x, Int_t n)
448 {
449 // J/Psi shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
450 //
451 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
452 //
453   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
454   const Double_t c[7] = {6.4922e-01, 6.4857e-03, 4.7268e-03,-9.5075e-04, 
455                          8.4075e-05,-4.2006e-06, 4.9970e-07};
456   Double_t y;
457   Int_t j;
458   y = c[j = 6];
459   while (j > 0) y  = y * x + c[--j];
460   y /= 1 + c[6]*TMath::Power(x,6);
461   //  
462   return 1 + (y-1)*f[n];
463 }
464
465 Double_t AliGenMUONlib::PtJpsiPPb8800(const Double_t *px, const Double_t *dummy)
466 {
467 // J/Psi pT
468 // pPb 8.8 TeV, minimum bias 0-100 %
469 //
470   return PtJpsiPPb8800ShFdummy(*px, 0) * PtJpsiPP8800(px, dummy);
471 }
472
473 Double_t AliGenMUONlib::PtJpsiPPb8800c1(const Double_t *px, const Double_t *dummy)
474 {
475 // J/Psi pT
476 // pPb 8.8 TeV, 1st centrality bin 0-20 %
477 //
478   return PtJpsiPPb8800ShFdummy(*px, 1) * PtJpsiPP8800(px, dummy);
479 }
480
481 Double_t AliGenMUONlib::PtJpsiPPb8800c2(const Double_t *px, const Double_t *dummy)
482 {
483 // J/Psi pT
484 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
485 //
486   return PtJpsiPPb8800ShFdummy(*px, 2) * PtJpsiPP8800(px, dummy);
487 }
488
489 Double_t AliGenMUONlib::PtJpsiPPb8800c3(const Double_t *px, const Double_t *dummy)
490 {
491 // J/Psi pT
492 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
493 //
494   return PtJpsiPPb8800ShFdummy(*px, 3) * PtJpsiPP8800(px, dummy);
495 }
496
497 Double_t AliGenMUONlib::PtJpsiPPb8800c4(const Double_t *px, const Double_t *dummy)
498 {
499 // J/Psi pT
500 // pPb 8.8 TeV, 4th centrality bin 60-100 %
501 //
502   return PtJpsiPPb8800ShFdummy(*px, 4) * PtJpsiPP8800(px, dummy);
503 }
504
505 Double_t AliGenMUONlib::PtJpsiPbP8800ShFdummy(Double_t x, Int_t n)
506 {
507 // J/Psi shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
508 //
509 // Pbp 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
510 //
511   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
512   const Double_t c[7] = {8.7562e-01, 2.1944e-02, 7.8509e-03,-1.3979e-03, 
513                          3.8513e-05, 4.2008e-06, 1.7088e-06};
514   Double_t y;
515   Int_t j;
516   y = c[j = 6];
517   while (j > 0) y  = y * x + c[--j];
518   y /= 1 + c[6]*TMath::Power(x,6);
519   //  
520   return 1 + (y-1)*f[n];
521 }
522
523 Double_t AliGenMUONlib::PtJpsiPbP8800(const Double_t *px, const Double_t *dummy)
524 {
525 // J/Psi pT
526 // Pbp 8.8 TeV, minimum bias 0-100 %
527 //
528   return PtJpsiPbP8800ShFdummy(*px, 0) * PtJpsiPP8800(px, dummy);
529 }
530
531 Double_t AliGenMUONlib::PtJpsiPbP8800c1(const Double_t *px, const Double_t *dummy)
532 {
533 // J/Psi pT
534 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
535 //
536   return PtJpsiPbP8800ShFdummy(*px, 1) * PtJpsiPP8800(px, dummy);
537 }
538
539 Double_t AliGenMUONlib::PtJpsiPbP8800c2(const Double_t *px, const Double_t *dummy)
540 {
541 // J/Psi pT
542 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
543 //
544   return PtJpsiPbP8800ShFdummy(*px, 2) * PtJpsiPP8800(px, dummy);
545 }
546
547 Double_t AliGenMUONlib::PtJpsiPbP8800c3(const Double_t *px, const Double_t *dummy)
548 {
549 // J/Psi pT
550 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
551 //
552   return PtJpsiPbP8800ShFdummy(*px, 3) * PtJpsiPP8800(px, dummy);
553 }
554
555 Double_t AliGenMUONlib::PtJpsiPbP8800c4(const Double_t *px, const Double_t *dummy)
556 {
557 // J/Psi pT
558 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
559 //
560   return PtJpsiPbP8800ShFdummy(*px, 4) * PtJpsiPP8800(px, dummy);
561 }
562
563 Double_t AliGenMUONlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/)
564 {
565 // J/Psi pT
566   const Double_t kpt0 = 4.;
567   const Double_t kxn  = 3.6;
568   Double_t x=*px;
569   //
570   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
571   return x/TMath::Power(pass1,kxn);
572 }
573
574 Double_t AliGenMUONlib::PtJpsiCDFscaled( const Double_t *px, const Double_t */*dummy*/)
575 {
576 // J/Psi pT
577 //
578 // PbPb 5.5 TeV
579 // scaled from CDF data at 2 TeV
580 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
581
582   const Double_t kpt0 = 5.100;
583   const Double_t kxn  = 4.102;
584   Double_t x=*px;
585   //
586   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
587   return x/TMath::Power(pass1,kxn);
588 }
589
590 Double_t AliGenMUONlib::PtJpsiCDFscaledPP( const Double_t *px, const Double_t */*dummy*/)
591 {
592 // J/Psi pT
593 //
594 // pp 14 TeV
595 // scaled from CDF data at 2 TeV
596
597   const Double_t kpt0 = 5.630;
598   const Double_t kxn  = 4.071;
599   Double_t x=*px;
600   //
601   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
602   return x/TMath::Power(pass1,kxn);
603 }
604
605 Double_t AliGenMUONlib::PtJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
606 {
607 // J/Psi pT
608 //
609 // pp 10 TeV
610 // scaled from CDF data at 2 TeV
611
612   const Double_t kpt0 = 5.334;
613   const Double_t kxn  = 4.071;
614   Double_t x=*px;
615   //
616   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
617   return x/TMath::Power(pass1,kxn);
618 }
619
620 Double_t AliGenMUONlib::PtJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
621 {
622 // J/Psi pT
623 //
624 // pp 8.8 TeV
625 // scaled from CDF data at 2 TeV
626 //
627   const Double_t kpt0 = 5.245;
628   const Double_t kxn  = 4.071;
629   Double_t x=*px;
630   //
631   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
632   return x/TMath::Power(pass1,kxn);
633 }
634
635 Double_t AliGenMUONlib::PtJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
636 {
637 // J/Psi pT
638 //
639 // pp 7 TeV
640 // scaled from CDF data at 2 TeV
641
642   const Double_t kpt0 = 5.072;
643   const Double_t kxn  = 4.071;
644   Double_t x=*px;
645   //
646   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
647   return x/TMath::Power(pass1,kxn);
648 }
649
650 Double_t AliGenMUONlib::PtJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
651 {
652 // J/Psi pT
653 //
654 // pp 3.94 TeV
655 // scaled from CDF data at 2 TeV
656 //
657   const Double_t kpt0 = 4.647;
658   const Double_t kxn  = 4.071;
659   Double_t x=*px;
660   //
661   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
662   return x/TMath::Power(pass1,kxn);
663 }
664
665 Double_t AliGenMUONlib::PtJpsiCDFscaledPP3( const Double_t *px, const Double_t */*dummy*/)
666 {
667 // J/Psi pT
668 //
669 // pp 2.76 TeV
670 // scaled from CDF data at 1.9 TeV
671 //
672   const Double_t kpt0 = 4.435;
673   const Double_t kxn  = 4.071;
674   Double_t x=*px;
675   //
676   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
677   return x/TMath::Power(pass1,kxn);
678 }
679
680 Double_t AliGenMUONlib::PtJpsiCDFscaledPP2( const Double_t *px, const Double_t */*dummy*/)
681 {
682 // J/Psi pT
683 //
684 // pp 1.96 TeV
685 // fit of the CDF data at 1.96 TeV
686 //
687   const Double_t kpt0 = 4.233;
688   const Double_t kxn  = 4.071;
689   Double_t x=*px;
690   //
691   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
692   return x/TMath::Power(pass1,kxn);
693 }
694
695 Double_t AliGenMUONlib::PtJpsiCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
696 {
697 // J/Psi pT
698 //
699 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
700 //
701   Double_t c[5] = {6.42774e-01, 1.86168e-02, -6.77296e-04, 8.93512e-06, 1.31586e-07};
702   Double_t x=*px;
703   Double_t y;
704   Int_t j;
705   y = c[j = 4];
706   while (j > 0) y  = y * x + c[--j];
707   //  
708   Double_t d = 1.+c[4]*TMath::Power(x,4);
709   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
710 }
711
712 Double_t AliGenMUONlib::PtJpsiCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
713 {
714 // J/Psi pT
715 //
716 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
717 //
718   Double_t c[5] = {8.58557e-01, 5.39791e-02, -4.75180e-03, 2.49463e-04, 5.52396e-05};
719   Double_t x=*px;
720   Double_t y;
721   Int_t j;
722   y = c[j = 4];
723   while (j > 0) y  = y * x + c[--j];
724   //  
725   Double_t d = 1.+c[4]*TMath::Power(x,4);
726   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
727 }
728
729 Double_t AliGenMUONlib::PtJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
730 {
731 // J/Psi pT
732 //
733 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
734 //
735   Double_t c[5] = {6.01022e-01, 4.70988e-02, -2.27917e-03, 3.09885e-05, 1.31955e-06};
736   Double_t x=*px;
737   Double_t y;
738   Int_t j;
739   y = c[j = 4];
740   while (j > 0) y  = y * x + c[--j];
741   //  
742   Double_t d = 1.+c[4]*TMath::Power(x,4);
743   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP4(px,dummy);
744 }
745
746 Double_t AliGenMUONlib::PtJpsiFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
747 {
748   return 1.;
749 }
750
751 Double_t AliGenMUONlib::PtJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
752 {
753 // J/Psi pT spectrum
754 //
755 // R. Vogt 2002
756 // PbPb 5.5 TeV
757 // MRST HO
758 // mc = 1.4 GeV, pt-kick 1 GeV
759 //
760     Float_t x = px[0];
761     Float_t c[8] = {
762         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00, 
763         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
764     };
765     
766     Double_t y;
767     if (x < 10.) {
768         Int_t j;
769         y = c[j = 7];
770         while (j > 0) y  = y * x +c[--j];
771         y = x * TMath::Exp(y);
772     } else {
773         y = 0.;
774     }
775     return y;
776 }
777
778 Double_t AliGenMUONlib::PtJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
779 {
780 // J/Psi pT spectrum
781 // B -> J/Psi X
782     Double_t x0 =   4.0384;
783     Double_t  n =   3.0288;
784     
785     Double_t x = px[0];
786     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
787     
788     return y;
789 }
790
791
792 Double_t AliGenMUONlib::PtJpsiPP( const Double_t *px, const Double_t */*dummy*/)
793 {
794 // J/Psi pT spectrum
795 //
796 // R. Vogt 2002
797 // pp 14 TeV
798 // MRST HO
799 // mc = 1.4 GeV, pt-kick 1 GeV
800 //
801     Float_t x = px[0];
802     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
803  
804     Double_t y;
805     if (x < 10.) {
806         Int_t j;
807         y = c[j = 3];
808         while (j > 0) y  = y * x +c[--j];
809         y = x * TMath::Exp(y);
810     } else {
811         y = 0.;
812     }
813     return y;
814 }
815
816 //
817 //               y-distribution
818 //____________________________________________________________
819 Double_t AliGenMUONlib::YJpsiPPdummy(Double_t x, Double_t energy)
820 {
821 // J/Psi y
822 // pp
823 // from the fit of RHIC + LHC data, see arXiv:1103.2394
824 //
825     x = x/TMath::Log(energy/3.097);
826     x = x*x;
827     Double_t y = TMath::Exp(-x/0.4/0.4/2);
828     if(x > 1) y=0;
829     return y;
830 }
831
832 Double_t AliGenMUONlib::YJpsiPPpoly(Double_t x, Double_t energy)
833 {
834 // J/Psi y
835 // pp
836 // from the fit of RHIC + LHC data, see arXiv:1103.2394
837 //
838     x = x/TMath::Log(energy/3.097);
839     x = x*x;
840     Double_t y = 1 - 6.9*x*x;
841     if(y < 0) y=0;
842     return y;
843 }
844
845 Double_t AliGenMUONlib::YJpsiPP7000(const Double_t *px, const Double_t */*dummy*/)
846 {
847 // J/Psi y
848 // pp 7 TeV
849 //
850   return YJpsiPPdummy(*px, 7000);
851 }
852
853 Double_t AliGenMUONlib::YJpsiPP2760(const Double_t *px, const Double_t */*dummy*/)
854 {
855 // J/Psi y
856 // pp 2.76 TeV
857 //
858   return YJpsiPPdummy(*px, 2760);
859 }
860
861 Double_t AliGenMUONlib::YJpsiPP4400(const Double_t *px, const Double_t */*dummy*/)
862 {
863 // J/Psi y
864 // pp 4.4 TeV
865 //
866   return YJpsiPPdummy(*px, 4400);
867 }
868
869 Double_t AliGenMUONlib::YJpsiPP5030(const Double_t *px, const Double_t */*dummy*/)
870 {
871 // J/Psi y
872 // pp 5.03 TeV
873 //
874   return YJpsiPPdummy(*px, 5030);
875 }
876
877 Double_t AliGenMUONlib::YJpsiPP8800(const Double_t *px, const Double_t */*dummy*/)
878 {
879 // J/Psi y
880 // pp 8.8 TeV
881 //
882   return YJpsiPPdummy(*px, 8800);
883 }
884
885 Double_t AliGenMUONlib::YJpsiPPpoly7000(const Double_t *px, const Double_t */*dummy*/)
886 {
887 // J/Psi y
888 // pp 7 TeV
889 //
890   return YJpsiPPpoly(*px, 7000);
891 }
892
893 Double_t AliGenMUONlib::YJpsiPPpoly2760(const Double_t *px, const Double_t */*dummy*/)
894 {
895 // J/Psi y
896 // pp 2.76 TeV
897 //
898   return YJpsiPPpoly(*px, 2760);
899 }
900
901 Double_t AliGenMUONlib::YJpsiPbPb2760ShFdummy(Double_t x, Int_t n)
902 {
903 // J/Psi shadowing factor vs y for PbPb min. bias and 11 centr. bins
904 //
905 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.66 in 4pi
906 //
907   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
908                            0.428, 0.317, 0.231, 0.156};
909   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
910                            0.106, 0.041, 0.013, 0.002};
911   const Double_t c1[5] = {1.5591e+00, 7.5827e-03, 2.0676e-03,-1.1717e-04, 1.5237e-06}; 
912   const Double_t c2[5] = {6.0861e-01, 4.8854e-03, 1.3685e-03,-7.9182e-05, 1.0475e-06}; 
913
914   x = x*x;
915   Double_t y1, y2;
916   Int_t j;
917   y1 = c1[j = 4]; y2 = c2[4];
918   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
919   
920   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
921   if(y1<0) y1=0;
922   return y1;
923 }
924
925 Double_t AliGenMUONlib::YJpsiPbPb2760(const Double_t *px, const Double_t *dummy)
926 {
927 // J/Psi y
928 // PbPb 2.76 TeV, minimum bias 0-100 %
929 //
930   return YJpsiPbPb2760ShFdummy(*px, 0) * YJpsiPP2760(px, dummy);
931 }
932
933 Double_t AliGenMUONlib::YJpsiPbPb2760c1(const Double_t *px, const Double_t *dummy)
934 {
935 // J/Psi y
936 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
937 //
938   return YJpsiPbPb2760ShFdummy(*px, 1) * YJpsiPP2760(px, dummy);
939 }
940
941 Double_t AliGenMUONlib::YJpsiPbPb2760c2(const Double_t *px, const Double_t *dummy)
942 {
943 // J/Psi y
944 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
945 //
946   return YJpsiPbPb2760ShFdummy(*px, 2) * YJpsiPP2760(px, dummy);
947 }
948
949 Double_t AliGenMUONlib::YJpsiPbPb2760c3(const Double_t *px, const Double_t *dummy)
950 {
951 // J/Psi y
952 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
953 //
954   return YJpsiPbPb2760ShFdummy(*px, 3) * YJpsiPP2760(px, dummy);
955 }
956
957 Double_t AliGenMUONlib::YJpsiPbPb2760c4(const Double_t *px, const Double_t *dummy)
958 {
959 // J/Psi y
960 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
961 //
962   return YJpsiPbPb2760ShFdummy(*px, 4) * YJpsiPP2760(px, dummy);
963 }
964
965 Double_t AliGenMUONlib::YJpsiPbPb2760c5(const Double_t *px, const Double_t *dummy)
966 {
967 // J/Psi y
968 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
969 //
970   return YJpsiPbPb2760ShFdummy(*px, 5) * YJpsiPP2760(px, dummy);
971 }
972
973 Double_t AliGenMUONlib::YJpsiPbPb2760c6(const Double_t *px, const Double_t *dummy)
974 {
975 // J/Psi y
976 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
977 //
978   return YJpsiPbPb2760ShFdummy(*px, 6) * YJpsiPP2760(px, dummy);
979 }
980
981 Double_t AliGenMUONlib::YJpsiPbPb2760c7(const Double_t *px, const Double_t *dummy)
982 {
983 // J/Psi y
984 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
985 //
986   return YJpsiPbPb2760ShFdummy(*px, 7) * YJpsiPP2760(px, dummy);
987 }
988
989 Double_t AliGenMUONlib::YJpsiPbPb2760c8(const Double_t *px, const Double_t *dummy)
990 {
991 // J/Psi y
992 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
993 //
994   return YJpsiPbPb2760ShFdummy(*px, 8) * YJpsiPP2760(px, dummy);
995 }
996
997 Double_t AliGenMUONlib::YJpsiPbPb2760c9(const Double_t *px, const Double_t *dummy)
998 {
999 // J/Psi y
1000 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
1001 //
1002   return YJpsiPbPb2760ShFdummy(*px, 9) * YJpsiPP2760(px, dummy);
1003 }
1004
1005 Double_t AliGenMUONlib::YJpsiPbPb2760c10(const Double_t *px, const Double_t *dummy)
1006 {
1007 // J/Psi y
1008 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
1009 //
1010   return YJpsiPbPb2760ShFdummy(*px, 10) * YJpsiPP2760(px, dummy);
1011 }
1012
1013 Double_t AliGenMUONlib::YJpsiPbPb2760c11(const Double_t *px, const Double_t *dummy)
1014 {
1015 // J/Psi y
1016 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
1017 //
1018   return YJpsiPbPb2760ShFdummy(*px, 11) * YJpsiPP2760(px, dummy);
1019 }
1020
1021 Double_t AliGenMUONlib::YJpsiPP5030dummy(Double_t px)
1022 {
1023   return AliGenMUONlib::YJpsiPP5030(&px, (Double_t*) 0);
1024 }
1025
1026 Double_t AliGenMUONlib::YJpsiPPb5030ShFdummy(Double_t x, Int_t n)
1027 {
1028 // J/Psi shadowing factor vs y for pPb min. bias and 4 centr. bins
1029 //
1030 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.81 in 4pi
1031 //
1032   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1033   const Double_t c[7] = {7.641e-01, 1.611e-02, 4.109e-03, 2.818e-03, 3.359e-04,-6.376e-05,-9.717e-06};
1034   Double_t y;
1035   Int_t j;
1036   y = c[j = 6];
1037   while (j > 0) y = y * x + c[--j];
1038   if(y<0) y=0;
1039   //
1040   return 1 +(y-1)*f[n];
1041 }
1042
1043 Double_t AliGenMUONlib::YJpsiPPb5030(const Double_t *px, const Double_t */*dummy*/)
1044 {
1045 // J/Psi y
1046 // pPb 5.03 TeV, minimum bias 0-100 %
1047 //
1048   Double_t x = px[0] + 0.47;              // rapidity shift
1049   return YJpsiPPb5030ShFdummy(x, 0) * YJpsiPP5030dummy(x);
1050 }
1051
1052 Double_t AliGenMUONlib::YJpsiPPb5030c1(const Double_t *px, const Double_t */*dummy*/)
1053 {
1054 // J/Psi y
1055 // pPb 5.03 TeV, 1st centrality bin 0-20 %
1056 //
1057   Double_t x = px[0] + 0.47;              // rapidity shift
1058   return YJpsiPPb5030ShFdummy(x, 1) * YJpsiPP5030dummy(x);
1059 }
1060
1061 Double_t AliGenMUONlib::YJpsiPPb5030c2(const Double_t *px, const Double_t */*dummy*/)
1062 {
1063 // J/Psi y
1064 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
1065 //
1066   Double_t x = px[0] + 0.47;              // rapidity shift
1067   return YJpsiPPb5030ShFdummy(x, 2) * YJpsiPP5030dummy(x);
1068 }
1069
1070 Double_t AliGenMUONlib::YJpsiPPb5030c3(const Double_t *px, const Double_t */*dummy*/)
1071 {
1072 // J/Psi y
1073 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
1074 //
1075   Double_t x = px[0] + 0.47;              // rapidity shift
1076   return YJpsiPPb5030ShFdummy(x, 3) * YJpsiPP5030dummy(x);
1077 }
1078
1079 Double_t AliGenMUONlib::YJpsiPPb5030c4(const Double_t *px, const Double_t */*dummy*/)
1080 {
1081 // J/Psi y
1082 // pPb 5.03 TeV, 4th centrality bin 60-100 %
1083 //
1084   Double_t x = px[0] + 0.47;              // rapidity shift
1085   return YJpsiPPb5030ShFdummy(x, 4) * YJpsiPP5030dummy(x);
1086 }
1087
1088 Double_t AliGenMUONlib::YJpsiPbP5030(const Double_t *px, const Double_t */*dummy*/)
1089 {
1090 // J/Psi y
1091 // Pbp 5.03 TeV, minimum bias 0-100 %
1092 //
1093   Double_t x = -px[0] + 0.47;              // rapidity shift
1094   return YJpsiPPb5030ShFdummy(x, 0) * YJpsiPP5030dummy(x);
1095 }
1096
1097 Double_t AliGenMUONlib::YJpsiPbP5030c1(const Double_t *px, const Double_t */*dummy*/)
1098 {
1099 // J/Psi y
1100 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
1101 //
1102   Double_t x = -px[0] + 0.47;              // rapidity shift
1103   return YJpsiPPb5030ShFdummy(x, 1) * YJpsiPP5030dummy(x);
1104 }
1105
1106 Double_t AliGenMUONlib::YJpsiPbP5030c2(const Double_t *px, const Double_t */*dummy*/)
1107 {
1108 // J/Psi y
1109 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
1110 //
1111   Double_t x = -px[0] + 0.47;              // rapidity shift
1112   return YJpsiPPb5030ShFdummy(x, 2) * YJpsiPP5030dummy(x);
1113 }
1114
1115 Double_t AliGenMUONlib::YJpsiPbP5030c3(const Double_t *px, const Double_t */*dummy*/)
1116 {
1117 // J/Psi y
1118 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
1119 //
1120   Double_t x = -px[0] + 0.47;              // rapidity shift
1121   return YJpsiPPb5030ShFdummy(x, 3) * YJpsiPP5030dummy(x);
1122 }
1123
1124 Double_t AliGenMUONlib::YJpsiPbP5030c4(const Double_t *px, const Double_t */*dummy*/)
1125 {
1126 // J/Psi y
1127 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
1128 //
1129   Double_t x = -px[0] + 0.47;              // rapidity shift
1130   return YJpsiPPb5030ShFdummy(x, 4) * YJpsiPP5030dummy(x);
1131 }
1132
1133 Double_t AliGenMUONlib::YJpsiPP8800dummy(Double_t px)
1134 {
1135     return AliGenMUONlib::YJpsiPP8800(&px, (Double_t*) 0);
1136 }
1137
1138 Double_t AliGenMUONlib::YJpsiPPb8800ShFdummy(Double_t x, Int_t n)
1139 {
1140 // J/Psi shadowing factor vs y for pPb min. bias and 4 centr. bins
1141 //
1142 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
1143 //
1144     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1145     const Double_t c[7] = {7.4372e-01, 2.3299e-02, 2.8678e-03, 1.9595e-03, 
1146                            3.2849e-04,-4.0547e-05,-7.9732e-06}; 
1147     Double_t y;
1148     Int_t j;
1149     y = c[j = 6];
1150     while (j > 0) y = y * x + c[--j];
1151     if(y<0) y=0;
1152     //
1153     return 1 +(y-1)*f[n];
1154 }
1155
1156 Double_t AliGenMUONlib::YJpsiPPb8800(const Double_t *px, const Double_t */*dummy*/)
1157 {
1158 // J/Psi y
1159 // pPb 8.8 TeV, minimum bias 0-100 %
1160 //
1161   Double_t x = px[0] + 0.47;              // rapidity shift
1162   return YJpsiPPb8800ShFdummy(x, 0) * YJpsiPP8800dummy(x);
1163 }
1164
1165 Double_t AliGenMUONlib::YJpsiPPb8800c1(const Double_t *px, const Double_t */*dummy*/)
1166 {
1167 // J/Psi y
1168 // pPb 8.8 TeV, 1st centrality bin 0-20 %
1169 //
1170   Double_t x = px[0] + 0.47;              // rapidity shift
1171   return YJpsiPPb8800ShFdummy(x, 1) * YJpsiPP8800dummy(x);
1172 }
1173
1174 Double_t AliGenMUONlib::YJpsiPPb8800c2(const Double_t *px, const Double_t */*dummy*/)
1175 {
1176 // J/Psi y
1177 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
1178 //
1179   Double_t x = px[0] + 0.47;              // rapidity shift
1180   return YJpsiPPb8800ShFdummy(x, 2) * YJpsiPP8800dummy(x);
1181 }
1182
1183 Double_t AliGenMUONlib::YJpsiPPb8800c3(const Double_t *px, const Double_t */*dummy*/)
1184 {
1185 // J/Psi y
1186 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
1187 //
1188   Double_t x = px[0] + 0.47;              // rapidity shift
1189   return YJpsiPPb8800ShFdummy(x, 3) * YJpsiPP8800dummy(x);
1190 }
1191
1192 Double_t AliGenMUONlib::YJpsiPPb8800c4(const Double_t *px, const Double_t */*dummy*/)
1193 {
1194 // J/Psi y
1195 // pPb 8.8 TeV, 4th centrality bin 60-100 %
1196 //
1197   Double_t x = px[0] + 0.47;              // rapidity shift
1198   return YJpsiPPb8800ShFdummy(x, 4) * YJpsiPP8800dummy(x);
1199 }
1200
1201 Double_t AliGenMUONlib::YJpsiPbP8800(const Double_t *px, const Double_t */*dummy*/)
1202 {
1203 // J/Psi y
1204 // Pbp 8.8 TeV, minimum bias 0-100 %
1205 //
1206   Double_t x = -px[0] + 0.47;              // rapidity shift
1207   return YJpsiPPb8800ShFdummy(x, 0) * YJpsiPP8800dummy(x);
1208 }
1209
1210 Double_t AliGenMUONlib::YJpsiPbP8800c1(const Double_t *px, const Double_t */*dummy*/)
1211 {
1212 // J/Psi y
1213 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
1214 //
1215   Double_t x = -px[0] + 0.47;              // rapidity shift
1216   return YJpsiPPb8800ShFdummy(x, 1) * YJpsiPP8800dummy(x);
1217 }
1218
1219 Double_t AliGenMUONlib::YJpsiPbP8800c2(const Double_t *px, const Double_t */*dummy*/)
1220 {
1221 // J/Psi y
1222 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
1223 //
1224   Double_t x = -px[0] + 0.47;              // rapidity shift
1225   return YJpsiPPb8800ShFdummy(x, 2) * YJpsiPP8800dummy(x);
1226 }
1227
1228 Double_t AliGenMUONlib::YJpsiPbP8800c3(const Double_t *px, const Double_t */*dummy*/)
1229 {
1230 // J/Psi y
1231 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
1232 //
1233   Double_t x = -px[0] + 0.47;              // rapidity shift
1234   return YJpsiPPb8800ShFdummy(x, 3) * YJpsiPP8800dummy(x);
1235 }
1236
1237 Double_t AliGenMUONlib::YJpsiPbP8800c4(const Double_t *px, const Double_t */*dummy*/)
1238 {
1239 // J/Psi y
1240 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
1241 //
1242   Double_t x = -px[0] + 0.47;              // rapidity shift
1243   return YJpsiPPb8800ShFdummy(x, 4) * YJpsiPP8800dummy(x);
1244 }
1245
1246 Double_t AliGenMUONlib::YJpsi(const Double_t *py, const Double_t */*dummy*/)
1247 {
1248 // J/psi y
1249   const Double_t ky0 = 4.;
1250   const Double_t kb=1.;
1251   Double_t yj;
1252   Double_t y=TMath::Abs(*py);
1253   //
1254   if (y < ky0)
1255     yj=kb;
1256   else
1257     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
1258   return yj;
1259 }
1260
1261 Double_t AliGenMUONlib::YJpsiFlat( const Double_t */*py*/, const Double_t */*dummy*/ )
1262 {
1263   return 1.;
1264 }
1265
1266
1267 Double_t AliGenMUONlib::YJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
1268 {
1269
1270 //
1271 // J/Psi y
1272 //
1273 //
1274 // R. Vogt 2002
1275 // PbPb 5.5 TeV
1276 // MRST HO
1277 // mc = 1.4 GeV, pt-kick 1 GeV
1278 //
1279     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
1280     Double_t x = TMath::Abs(px[0]);
1281     Double_t y;
1282     
1283     if (x < 4.) {
1284         y = 31.754;
1285     } else if (x < 6) {
1286         Int_t j;
1287         y = c[j = 4];
1288         while (j > 0) y  = y * x + c[--j];
1289     } else {
1290         y =0.;
1291     }
1292     
1293     return y;
1294 }
1295
1296 Double_t AliGenMUONlib::YJpsiCDFscaled( const Double_t *px, const Double_t* dummy)
1297 {
1298     // J/Psi y 
1299     return AliGenMUONlib::YJpsiPbPb(px, dummy);
1300 }
1301
1302 Double_t AliGenMUONlib::YJpsiCDFscaledPP( const Double_t *px, const Double_t* dummy)
1303 {
1304     // J/Psi y 
1305     return AliGenMUONlib::YJpsiPP(px, dummy);
1306 }
1307
1308 Double_t AliGenMUONlib::YJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
1309 {
1310 // J/Psi y
1311 //
1312 // pp 10 TeV
1313 // scaled from YJpsiPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
1314 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
1315 //
1316
1317     Double_t c[5] = {2.46681e+01, 8.91486e+01, -3.21227e+01, 3.63075e+00, -1.32047e-01};
1318
1319     Double_t x = TMath::Abs(px[0]);
1320     Double_t y;
1321
1322     if (x < 3.2) {
1323         y = 98.523 - 1.3664 * x * x;
1324     } else if (x < 7.5) {
1325         Int_t j;
1326         y = c[j = 4];
1327         while (j > 0) y  = y * x + c[--j];
1328     } else {
1329         y =0.;
1330     }
1331
1332     if(y<0) y=0;
1333
1334     return y;
1335 }
1336
1337 Double_t AliGenMUONlib::YJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
1338 {
1339 // J/Psi y
1340 //
1341 // pp 8.8 TeV
1342 // rescaling of YJpsiPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD 
1343 //
1344     Double_t c[5] = {3.33882e+02, -1.30980e+02, 2.59082e+01, -3.08935e+00, 1.56375e-01};
1345     Double_t x = TMath::Abs(px[0]);
1346     Double_t y;
1347
1348     if (x < 3.7) {
1349         y = 99.236 - 1.5498 * x * x;
1350     } else if (x < 7.4) {
1351         Int_t j;
1352         y = c[j = 4];
1353         while (j > 0) y  = y * x + c[--j];
1354     } else {
1355         y =0.;
1356     }
1357
1358     if(y<0) y=0;
1359
1360     return y;
1361 }
1362
1363 Double_t AliGenMUONlib::YJpsiCDFscaledPP9dummy(Double_t px)
1364 {
1365     return AliGenMUONlib::YJpsiCDFscaledPP9(&px, (Double_t*) 0);
1366 }
1367
1368 Double_t AliGenMUONlib::YJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
1369 {
1370 // J/Psi y
1371 //
1372 // pp 7 TeV
1373 // scaled from YJpsiPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
1374 //
1375
1376     Double_t c[5] = {6.71181e+02, -3.69240e+02, 8.89644e+01, -1.04937e+01, 4.80959e-01};
1377
1378     Double_t x = TMath::Abs(px[0]);
1379     Double_t y;
1380
1381     if (x < 4.0) {
1382         y = 100.78 - 1.8353 * x * x;
1383     } else if (x < 7.3) {
1384         Int_t j;
1385         y = c[j = 4];
1386         while (j > 0) y  = y * x + c[--j];
1387     } else {
1388         y =0.;
1389     }
1390
1391     if(y<0) y=0;
1392
1393     return y;
1394 }
1395
1396 Double_t AliGenMUONlib::YJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
1397 {
1398 // J/Psi y
1399 //
1400 // pp 3.94 TeV
1401 // rescaling of YJpsiPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD 
1402 //
1403     Double_t c[5] = {4.00785e+02, -1.41159e+01, -3.28599e+01, 5.53048e+00, -2.45151e-01};
1404     Double_t x = TMath::Abs(px[0]);
1405     Double_t y;
1406
1407     if (x < 5.5) {
1408         y = 107.389 - 2.7454 * x * x;
1409     } else if (x < 7.0) {
1410         Int_t j;
1411         y = c[j = 4];
1412         while (j > 0) y  = y * x + c[--j];
1413     } else {
1414         y =0.;
1415     }
1416
1417     if(y<0) y=0;
1418
1419     return y;
1420 }
1421
1422 Double_t AliGenMUONlib::YJpsiCDFscaledPP3( const Double_t *px, const Double_t *dummy)
1423 {
1424 // J/Psi y 
1425     return AliGenMUONlib::YJpsiPP2760(px, dummy);
1426 }
1427
1428 Double_t AliGenMUONlib::YJpsiCDFscaledPP2( const Double_t *px, const Double_t */*dummy*/)
1429 {
1430 // J/Psi y
1431 // pp 1.96 TeV
1432 //
1433   return YJpsiPPdummy(*px, 1960);
1434 }
1435
1436 Double_t AliGenMUONlib::YJpsiPP( const Double_t *px, const Double_t */*dummy*/)
1437 {
1438
1439 //
1440 // J/Psi y
1441 //
1442 //
1443 // R. Vogt 2002
1444 // pp 14  TeV
1445 // MRST HO
1446 // mc = 1.4 GeV, pt-kick 1 GeV
1447 //
1448
1449     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
1450     Double_t x = TMath::Abs(px[0]);
1451     Double_t y;
1452     
1453     if (x < 2.5) {
1454         y = 96.455 - 0.8483 * x * x;
1455     } else if (x < 7.9) {
1456         Int_t j;
1457         y = c[j = 4];
1458         while (j > 0) y  = y * x + c[--j];
1459     } else {
1460         y =0.;
1461     }
1462     
1463     return y;
1464 }
1465
1466 Double_t AliGenMUONlib::YJpsiCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
1467 {
1468 // J/Psi y
1469 //
1470 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
1471 //
1472     Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
1473                      -4.16509e-05,-7.62709e-06}; 
1474     Double_t y;
1475     Double_t x = px[0] + 0.47;              // rapidity shift
1476     Int_t j;
1477     y = c[j = 6];
1478     while (j > 0) y = y * x + c[--j];
1479     if(y<0) y=0;
1480
1481     return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
1482 }
1483
1484 Double_t AliGenMUONlib::YJpsiCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
1485 {
1486 // J/Psi y
1487 //
1488 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
1489 //
1490     Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
1491                      -4.16509e-05,-7.62709e-06}; 
1492     Double_t y;
1493     Double_t x = -px[0] + 0.47;              // rapidity shift
1494     Int_t j;
1495     y = c[j = 6];
1496     while (j > 0) y = y * x + c[--j];
1497     if(y<0) y=0;
1498
1499     return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
1500 }
1501
1502 Double_t AliGenMUONlib::YJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
1503 {
1504 // J/Psi y
1505 //
1506 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
1507 //
1508     Double_t c[4] = {5.95228e-01, 9.45069e-03, 2.44710e-04, -1.32894e-05}; 
1509     Double_t x = px[0]*px[0];
1510     Double_t y;
1511     Int_t j;
1512     y = c[j = 3];
1513     while (j > 0) y  = y * x + c[--j];
1514     if(y<0) y=0;
1515
1516     return y * AliGenMUONlib::YJpsiCDFscaledPP4(px,dummy);
1517 }
1518
1519 Double_t AliGenMUONlib::YJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
1520 {
1521
1522 //
1523 // J/Psi from B->J/Psi X
1524 //
1525 //
1526     
1527
1528     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
1529     
1530     Double_t x = TMath::Abs(px[0]);
1531     Double_t y;
1532     
1533     if (x > 6.) {
1534         y = 0.;
1535     } else {
1536         Int_t j;
1537         y = c[j = 6];
1538         while (j > 0) y  = y * x + c[--j];
1539     } 
1540     
1541     return y;
1542 }
1543
1544
1545
1546 //                 particle composition
1547 //
1548 Int_t AliGenMUONlib::IpJpsi(TRandom *)
1549 {
1550 // J/Psi composition
1551     return 443;
1552 }
1553 Int_t AliGenMUONlib::IpPsiP(TRandom *)
1554 {
1555 // Psi prime composition
1556     return 100443;
1557 }
1558 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
1559 {
1560 // J/Psi composition
1561   Int_t ip;
1562   Float_t r = gRandom->Rndm();
1563   if (r < 0.98) {
1564     ip = 443;
1565   } else {
1566     ip = 100443;
1567   }
1568   return ip;
1569 }
1570
1571
1572
1573 //                      Upsilon
1574 //
1575 //
1576 //                  pt-distribution
1577 //____________________________________________________________
1578 Double_t AliGenMUONlib::PtUpsilonPPdummy(Double_t x, Double_t energy)
1579 {
1580 // Upsilon pT
1581 // pp
1582 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
1583 //
1584   const Double_t kpt0 = 1.96*TMath::Power(energy,0.095);
1585   const Double_t kxn  = 3.4;
1586   //
1587   Double_t pass1 = 1.+0.471*(x/kpt0)*(x/kpt0);
1588   return x/TMath::Power(pass1,kxn);
1589 }
1590
1591 Double_t AliGenMUONlib::PtUpsilonPP7000(const Double_t *px, const Double_t */*dummy*/)
1592 {
1593 // Upsilon pT
1594 // pp 7 TeV
1595 //
1596   return PtUpsilonPPdummy(*px,7000);
1597 }
1598
1599 Double_t AliGenMUONlib::PtUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
1600 {
1601 // Upsilon pT
1602 // pp 2.76 TeV
1603 //
1604   return PtUpsilonPPdummy(*px,2760);
1605 }
1606
1607 Double_t AliGenMUONlib::PtUpsilonPP4400(const Double_t *px, const Double_t */*dummy*/)
1608 {
1609 // Upsilon pT
1610 // pp 4.4 TeV
1611 //
1612   return PtUpsilonPPdummy(*px,4400);
1613 }
1614
1615 Double_t AliGenMUONlib::PtUpsilonPP5030(const Double_t *px, const Double_t */*dummy*/)
1616 {
1617 // Upsilon pT
1618 // pp 5.03 TeV
1619 //
1620   return PtUpsilonPPdummy(*px,5030);
1621 }
1622
1623 Double_t AliGenMUONlib::PtUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
1624 {
1625 // Upsilon pT
1626 // pp 8.8 TeV
1627 //
1628   return PtUpsilonPPdummy(*px,8800);
1629 }
1630
1631 Double_t AliGenMUONlib::PtUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
1632 {
1633 // Usilon shadowing factor vs pT for PbPb min. bias and 11 centr. bins (in 2.5<y<4)
1634 //
1635 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
1636 //
1637   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
1638                            0.428, 0.317, 0.231, 0.156};
1639   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
1640                            0.106, 0.041, 0.013, 0.002};
1641   const Double_t c1[7] = {1.9089e+00, 1.2969e-03, 8.9786e-05,-5.3062e-06,
1642                           -1.0046e-06,6.1446e-08, 1.0885e-09};
1643   const Double_t c2[7] = {8.8423e-01,-8.7488e-05, 5.9857e-04,-5.7959e-05, 
1644                           2.0059e-06,-2.7343e-08, 6.6053e-10};
1645   Double_t y1, y2;
1646   Int_t j;
1647   y1 = c1[j = 6]; y2 = c2[6];
1648   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
1649   
1650   y1 /= 1.+c1[6]*TMath::Power(x,6);
1651   y2 /= 1.+c2[6]*TMath::Power(x,6);
1652   //  
1653   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
1654   if(y1<0) y1=0;
1655   return y1;
1656 }
1657
1658 Double_t AliGenMUONlib::PtUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
1659 {
1660 // Upsilon pT
1661 // PbPb 2.76 TeV, minimum bias 0-100 %
1662 //
1663   return PtUpsilonPbPb2760ShFdummy(*px, 0) * PtUpsilonPP2760(px, dummy);
1664 }
1665
1666 Double_t AliGenMUONlib::PtUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
1667 {
1668 // Upsilon pT
1669 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
1670 //
1671   return PtUpsilonPbPb2760ShFdummy(*px, 1) * PtUpsilonPP2760(px, dummy);
1672 }
1673
1674 Double_t AliGenMUONlib::PtUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
1675 {
1676 // Upsilon pT
1677 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
1678 //
1679   return PtUpsilonPbPb2760ShFdummy(*px, 2) * PtUpsilonPP2760(px, dummy);
1680 }
1681
1682 Double_t AliGenMUONlib::PtUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
1683 {
1684 // Upsilon pT
1685 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
1686 //
1687   return PtUpsilonPbPb2760ShFdummy(*px, 3) * PtUpsilonPP2760(px, dummy);
1688 }
1689
1690 Double_t AliGenMUONlib::PtUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
1691 {
1692 // Upsilon pT
1693 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
1694 //
1695   return PtUpsilonPbPb2760ShFdummy(*px, 4) * PtUpsilonPP2760(px, dummy);
1696 }
1697
1698 Double_t AliGenMUONlib::PtUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
1699 {
1700 // Upsilon pT
1701 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
1702 //
1703   return PtUpsilonPbPb2760ShFdummy(*px, 5) * PtUpsilonPP2760(px, dummy);
1704 }
1705
1706 Double_t AliGenMUONlib::PtUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
1707 {
1708 // Upsilon pT
1709 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
1710 //
1711   return PtUpsilonPbPb2760ShFdummy(*px, 6) * PtUpsilonPP2760(px, dummy);
1712 }
1713
1714 Double_t AliGenMUONlib::PtUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
1715 {
1716 // Upsilon pT
1717 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
1718 //
1719   return PtUpsilonPbPb2760ShFdummy(*px, 7) * PtUpsilonPP2760(px, dummy);
1720 }
1721
1722 Double_t AliGenMUONlib::PtUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
1723 {
1724 // Upsilon pT
1725 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
1726 //
1727   return PtUpsilonPbPb2760ShFdummy(*px, 8) * PtUpsilonPP2760(px, dummy);
1728 }
1729
1730 Double_t AliGenMUONlib::PtUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
1731 {
1732 // Upsilon pT
1733 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
1734 //
1735   return PtUpsilonPbPb2760ShFdummy(*px, 9) * PtUpsilonPP2760(px, dummy);
1736 }
1737
1738 Double_t AliGenMUONlib::PtUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
1739 {
1740 // Upsilon pT
1741 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
1742 //
1743   return PtUpsilonPbPb2760ShFdummy(*px, 10) * PtUpsilonPP2760(px, dummy);
1744 }
1745
1746 Double_t AliGenMUONlib::PtUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
1747 {
1748 // Upsilon pT
1749 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
1750 //
1751   return PtUpsilonPbPb2760ShFdummy(*px, 11) * PtUpsilonPP2760(px, dummy);
1752 }
1753
1754 Double_t AliGenMUONlib::PtUpsilonPPb5030ShFdummy(Double_t x, Int_t n)
1755 {
1756 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1757 //
1758 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
1759 //
1760   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1761   const Double_t c[5] = {8.069e-01, 1.407e-04, 4.372e-04,-2.797e-05, 4.405e-06};
1762   Double_t y;
1763   Int_t j;
1764   y = c[j = 4];
1765   while (j > 0) y  = y * x + c[--j];
1766   y /= 1 + c[4]*TMath::Power(x,4);
1767   //  
1768   return 1 + (y-1)*f[n];
1769 }
1770
1771 Double_t AliGenMUONlib::PtUpsilonPPb5030(const Double_t *px, const Double_t *dummy)
1772 {
1773 // Upsilon pT
1774 // pPb 5.03 TeV, minimum bias 0-100 %
1775 //
1776   return PtUpsilonPPb5030ShFdummy(*px, 0) * PtUpsilonPP5030(px, dummy);
1777 }
1778
1779 Double_t AliGenMUONlib::PtUpsilonPPb5030c1(const Double_t *px, const Double_t *dummy)
1780 {
1781 // Upsilon pT
1782 // pPb 5.03 TeV, 1st centrality bin 0-20 %
1783 //
1784   return PtUpsilonPPb5030ShFdummy(*px, 1) * PtUpsilonPP5030(px, dummy);
1785 }
1786
1787 Double_t AliGenMUONlib::PtUpsilonPPb5030c2(const Double_t *px, const Double_t *dummy)
1788 {
1789 // Upsilon pT
1790 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
1791 //
1792   return PtUpsilonPPb5030ShFdummy(*px, 2) * PtUpsilonPP5030(px, dummy);
1793 }
1794
1795 Double_t AliGenMUONlib::PtUpsilonPPb5030c3(const Double_t *px, const Double_t *dummy)
1796 {
1797 // Upsilon pT
1798 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
1799 //
1800   return PtUpsilonPPb5030ShFdummy(*px, 3) * PtUpsilonPP5030(px, dummy);
1801 }
1802
1803 Double_t AliGenMUONlib::PtUpsilonPPb5030c4(const Double_t *px, const Double_t *dummy)
1804 {
1805 // Upsilon pT
1806 // pPb 5.03 TeV, 4th centrality bin 60-100 %
1807 //
1808   return PtUpsilonPPb5030ShFdummy(*px, 4) * PtUpsilonPP5030(px, dummy);
1809 }
1810
1811 Double_t AliGenMUONlib::PtUpsilonPbP5030ShFdummy(Double_t x, Int_t n)
1812 {
1813 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1814 //
1815 // Pbp 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
1816 //
1817   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1818   const Double_t c[5] = {1.122, 2.565e-03,-3.025e-04, 4.113e-06, 6.140e-07};
1819   Double_t y;
1820   Int_t j;
1821   y = c[j = 4];
1822   while (j > 0) y  = y * x + c[--j];
1823   y /= 1 + c[4]*TMath::Power(x,4);
1824   //  
1825   return 1 + (y-1)*f[n];
1826 }
1827
1828 Double_t AliGenMUONlib::PtUpsilonPbP5030(const Double_t *px, const Double_t *dummy)
1829 {
1830 // Upsilon pT
1831 // Pbp 5.03 TeV, minimum bias 0-100 %
1832 //
1833   return PtUpsilonPbP5030ShFdummy(*px, 0) * PtUpsilonPP5030(px, dummy);
1834 }
1835
1836 Double_t AliGenMUONlib::PtUpsilonPbP5030c1(const Double_t *px, const Double_t *dummy)
1837 {
1838 // Upsilon pT
1839 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
1840 //
1841   return PtUpsilonPbP5030ShFdummy(*px, 1) * PtUpsilonPP5030(px, dummy);
1842 }
1843
1844 Double_t AliGenMUONlib::PtUpsilonPbP5030c2(const Double_t *px, const Double_t *dummy)
1845 {
1846 // Upsilon pT
1847 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
1848 //
1849   return PtUpsilonPbP5030ShFdummy(*px, 2) * PtUpsilonPP5030(px, dummy);
1850 }
1851
1852 Double_t AliGenMUONlib::PtUpsilonPbP5030c3(const Double_t *px, const Double_t *dummy)
1853 {
1854 // Upsilon pT
1855 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
1856 //
1857   return PtUpsilonPbP5030ShFdummy(*px, 3) * PtUpsilonPP5030(px, dummy);
1858 }
1859
1860 Double_t AliGenMUONlib::PtUpsilonPbP5030c4(const Double_t *px, const Double_t *dummy)
1861 {
1862 // Upsilon pT
1863 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
1864 //
1865   return PtUpsilonPbP5030ShFdummy(*px, 4) * PtUpsilonPP5030(px, dummy);
1866 }
1867
1868 Double_t AliGenMUONlib::PtUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
1869 {
1870 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1871 //
1872 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1873 //
1874   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1875   const Double_t c[5] = {7.6561e-01, 1.1360e-04, 4.9596e-04,-3.0287e-05, 3.7555e-06};
1876   Double_t y;
1877   Int_t j;
1878   y = c[j = 4];
1879   while (j > 0) y  = y * x + c[--j];
1880   y /= 1 + c[4]*TMath::Power(x,4);
1881   //  
1882   return 1 + (y-1)*f[n];
1883 }
1884
1885 Double_t AliGenMUONlib::PtUpsilonPPb8800(const Double_t *px, const Double_t *dummy)
1886 {
1887 // Upsilon pT
1888 // pPb 8.8 TeV, minimum bias 0-100 %
1889 //
1890   return PtUpsilonPPb8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1891 }
1892
1893 Double_t AliGenMUONlib::PtUpsilonPPb8800c1(const Double_t *px, const Double_t *dummy)
1894 {
1895 // Upsilon pT
1896 // pPb 8.8 TeV, 1st centrality bin 0-20 %
1897 //
1898   return PtUpsilonPPb8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1899 }
1900
1901 Double_t AliGenMUONlib::PtUpsilonPPb8800c2(const Double_t *px, const Double_t *dummy)
1902 {
1903 // Upsilon pT
1904 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
1905 //
1906   return PtUpsilonPPb8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1907 }
1908
1909 Double_t AliGenMUONlib::PtUpsilonPPb8800c3(const Double_t *px, const Double_t *dummy)
1910 {
1911 // Upsilon pT
1912 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
1913 //
1914   return PtUpsilonPPb8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1915 }
1916
1917 Double_t AliGenMUONlib::PtUpsilonPPb8800c4(const Double_t *px, const Double_t *dummy)
1918 {
1919 // Upsilon pT
1920 // pPb 8.8 TeV, 4th centrality bin 60-100 %
1921 //
1922   return PtUpsilonPPb8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1923 }
1924
1925 Double_t AliGenMUONlib::PtUpsilonPbP8800ShFdummy(Double_t x, Int_t n)
1926 {
1927 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1928 //
1929 // Pbp 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1930 //
1931   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1932   const Double_t c[5] = {1.0975, 3.1905e-03,-2.0477e-04, 8.5270e-06, 2.5343e-06};
1933   Double_t y;
1934   Int_t j;
1935   y = c[j = 4];
1936   while (j > 0) y  = y * x + c[--j];
1937   y /= 1 + c[4]*TMath::Power(x,4);
1938   //  
1939   return 1 + (y-1)*f[n];
1940 }
1941
1942 Double_t AliGenMUONlib::PtUpsilonPbP8800(const Double_t *px, const Double_t *dummy)
1943 {
1944 // Upsilon pT
1945 // Pbp 8.8 TeV, minimum bias 0-100 %
1946 //
1947   return PtUpsilonPbP8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1948 }
1949
1950 Double_t AliGenMUONlib::PtUpsilonPbP8800c1(const Double_t *px, const Double_t *dummy)
1951 {
1952 // Upsilon pT
1953 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
1954 //
1955   return PtUpsilonPbP8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1956 }
1957
1958 Double_t AliGenMUONlib::PtUpsilonPbP8800c2(const Double_t *px, const Double_t *dummy)
1959 {
1960 // Upsilon pT
1961 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
1962 //
1963   return PtUpsilonPbP8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1964 }
1965
1966 Double_t AliGenMUONlib::PtUpsilonPbP8800c3(const Double_t *px, const Double_t *dummy)
1967 {
1968 // Upsilon pT
1969 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
1970 //
1971   return PtUpsilonPbP8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1972 }
1973
1974 Double_t AliGenMUONlib::PtUpsilonPbP8800c4(const Double_t *px, const Double_t *dummy)
1975 {
1976 // Upsilon pT
1977 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
1978 //
1979   return PtUpsilonPbP8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1980 }
1981
1982 Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
1983 {
1984 // Upsilon pT
1985   const Double_t kpt0 = 5.3;
1986   const Double_t kxn  = 2.5;
1987   Double_t x=*px;
1988   //
1989   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1990   return x/TMath::Power(pass1,kxn);
1991 }
1992
1993 Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
1994 {
1995 // Upsilon pT
1996   const Double_t kpt0 = 7.753;
1997   const Double_t kxn  = 3.042;
1998   Double_t x=*px;
1999   //
2000   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2001   return x/TMath::Power(pass1,kxn);
2002 }
2003
2004 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
2005 {
2006 // Upsilon pT
2007 //
2008 // pp 14 TeV
2009 //
2010 // scaled from CDF data at 2 TeV
2011
2012   const Double_t kpt0 = 8.610;
2013   const Double_t kxn  = 3.051;
2014   Double_t x=*px;
2015   //
2016   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2017   return x/TMath::Power(pass1,kxn);
2018 }
2019
2020 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2021 {
2022 // Upsilon pT
2023 //
2024 // pp 10 TeV
2025 //
2026 // scaled from CDF data at 2 TeV
2027
2028   const Double_t kpt0 = 8.235;
2029   const Double_t kxn  = 3.051;
2030   Double_t x=*px;
2031   //
2032   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2033   return x/TMath::Power(pass1,kxn);
2034 }
2035
2036 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2037 {
2038 // Upsilon pT
2039 //
2040 // pp 8.8 TeV
2041 // scaled from CDF data at 2 TeV
2042 //
2043   const Double_t kpt0 = 8.048;
2044   const Double_t kxn  = 3.051;
2045   Double_t x=*px;
2046   //
2047   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2048   return x/TMath::Power(pass1,kxn);
2049 }
2050
2051 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2052 {
2053 // Upsilon pT
2054 //
2055 // pp 7 TeV
2056 //
2057 // scaled from CDF data at 2 TeV
2058
2059   const Double_t kpt0 = 7.817;
2060   const Double_t kxn  = 3.051;
2061   Double_t x=*px;
2062   //
2063   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2064   return x/TMath::Power(pass1,kxn);
2065 }
2066
2067 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2068 {
2069 // Upsilon pT
2070 //
2071 // pp 3.94 TeV
2072 // scaled from CDF data at 2 TeV
2073 //
2074   const Double_t kpt0 = 7.189;
2075   const Double_t kxn  = 3.051;
2076   Double_t x=*px;
2077   //
2078   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2079   return x/TMath::Power(pass1,kxn);
2080 }
2081
2082 Double_t AliGenMUONlib::PtUpsilonCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
2083 {
2084 // Upsilon pT
2085 //
2086 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2087 //
2088   Double_t c[5] = {7.64952e-01, 1.12501e-04, 4.96038e-04, -3.03198e-05, 3.74035e-06};
2089   Double_t x=*px;
2090   Double_t y;
2091   Int_t j;
2092   y = c[j = 4];
2093   while (j > 0) y  = y * x + c[--j];
2094   //  
2095   Double_t d = 1.+c[4]*TMath::Power(x,4);
2096   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
2097 }
2098
2099 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
2100 {
2101 // Upsilon pT
2102 //
2103 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2104 //
2105   Double_t c[5] = {1.09881e+00, 3.08329e-03, -2.00356e-04, 8.28991e-06, 2.52576e-06};
2106   Double_t x=*px;
2107   Double_t y;
2108   Int_t j;
2109   y = c[j = 4];
2110   while (j > 0) y  = y * x + c[--j];
2111   //  
2112   Double_t d = 1.+c[4]*TMath::Power(x,4);
2113   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
2114 }
2115
2116 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2117 {
2118 // Upsilon pT
2119 //
2120 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2121 //
2122   Double_t c[5] = {8.65872e-01, 2.05465e-03, 2.56063e-04, -1.65598e-05, 2.29209e-06};
2123   Double_t x=*px;
2124   Double_t y;
2125   Int_t j;
2126   y = c[j = 4];
2127   while (j > 0) y  = y * x + c[--j];
2128   //  
2129   Double_t d = 1.+c[4]*TMath::Power(x,4);
2130   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP4(px,dummy);
2131 }
2132
2133 Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
2134 {
2135   return 1.;
2136 }
2137
2138 Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2139 {
2140 //
2141 // Upsilon pT
2142 //
2143 //
2144 // R. Vogt 2002
2145 // PbPb 5.5 TeV
2146 // MRST HO
2147 // mc = 1.4 GeV, pt-kick 1 GeV
2148 //
2149     Float_t x = px[0];
2150     Double_t c[8] = {
2151         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,       
2152         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
2153     };
2154     Double_t y;
2155     if (x < 10.) {
2156         Int_t j;
2157         y = c[j = 7];
2158         while (j > 0) y  = y * x +c[--j];
2159         y = x * TMath::Exp(y);
2160     } else {
2161         y = 0.;
2162     }
2163     return y;
2164 }
2165
2166 Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2167 {
2168 //
2169 // Upsilon pT
2170 //
2171 //
2172 // R. Vogt 2002
2173 // pp 14 TeV
2174 // MRST HO
2175 // mc = 1.4 GeV, pt-kick 1 GeV
2176 //
2177     Float_t x = px[0];
2178     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,   
2179                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
2180     
2181     Double_t y;
2182     if (x < 10.) {
2183         Int_t j;
2184         y = c[j = 7];
2185         while (j > 0) y  = y * x +c[--j];
2186         y = x * TMath::Exp(y);
2187     } else {
2188         y = 0.;
2189     }
2190     return y;
2191 }
2192
2193 //
2194 //                    y-distribution
2195 //
2196 //____________________________________________________________
2197 Double_t AliGenMUONlib::YUpsilonPPdummy(Double_t x, Double_t energy)
2198 {
2199 // Upsilon y
2200 // pp
2201 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
2202 //
2203     x = x/TMath::Log(energy/9.46);
2204     x = x*x;
2205     Double_t y = TMath::Exp(-x/0.4/0.4/2);
2206     if(x > 1) y=0;
2207     return y;
2208 }
2209
2210 Double_t AliGenMUONlib::YUpsilonPPpoly(Double_t x, Double_t energy)
2211 {
2212 // Upsilon y
2213 // pp
2214 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
2215 //
2216     x = x/TMath::Log(energy/9.46);
2217     x = x*x;
2218     Double_t y = 1 - 6.9*x*x;
2219     if(y < 0) y=0;
2220     return y;
2221 }
2222
2223 Double_t AliGenMUONlib::YUpsilonPP7000(const Double_t *px, const Double_t */*dummy*/)
2224 {
2225 // Upsilon y
2226 // pp 7 TeV
2227 //
2228   return YUpsilonPPdummy(*px, 7000);
2229 }
2230
2231 Double_t AliGenMUONlib::YUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
2232 {
2233 // Upsilon y
2234 // pp 2.76 TeV
2235 //
2236   return YUpsilonPPdummy(*px, 2760);
2237 }
2238
2239 Double_t AliGenMUONlib::YUpsilonPP4400(const Double_t *px, const Double_t */*dummy*/)
2240 {
2241 // Upsilon y
2242 // pp 4.4 TeV
2243 //
2244   return YUpsilonPPdummy(*px, 4400);
2245 }
2246
2247 Double_t AliGenMUONlib::YUpsilonPP5030(const Double_t *px, const Double_t */*dummy*/)
2248 {
2249 // Upsilon y
2250 // pp 5.03 TeV
2251 //
2252   return YUpsilonPPdummy(*px, 5030);
2253 }
2254
2255 Double_t AliGenMUONlib::YUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
2256 {
2257 // Upsilon y
2258 // pp 8.8 TeV
2259 //
2260   return YUpsilonPPdummy(*px, 8800);
2261 }
2262
2263 Double_t AliGenMUONlib::YUpsilonPPpoly7000(const Double_t *px, const Double_t */*dummy*/)
2264 {
2265 // Upsilon y
2266 // pp 7 TeV
2267 //
2268   return YUpsilonPPpoly(*px, 7000);
2269 }
2270
2271 Double_t AliGenMUONlib::YUpsilonPPpoly2760(const Double_t *px, const Double_t */*dummy*/)
2272 {
2273 // Upsilon y
2274 // pp 2.76 TeV
2275 //
2276   return YUpsilonPPpoly(*px, 2760);
2277 }
2278
2279 Double_t AliGenMUONlib::YUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
2280 {
2281 // Upsilon shadowing factor vs y for PbPb min. bias and 11 centr. bins
2282 //
2283 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
2284 //
2285   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
2286                            0.428, 0.317, 0.231, 0.156};
2287   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
2288                            0.106, 0.041, 0.013, 0.002};
2289   const Double_t c1[5] = {1.8547e+00, 1.6717e-02,-2.1285e-04,-9.7662e-05, 2.5768e-06};
2290   const Double_t c2[5] = {8.6029e-01, 1.1742e-02,-2.7944e-04,-6.7973e-05, 1.8838e-06}; 
2291
2292   x = x*x;
2293   Double_t y1, y2;
2294   Int_t j;
2295   y1 = c1[j = 4]; y2 = c2[4];
2296   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
2297   
2298   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
2299   if(y1<0) y1=0;
2300   return y1;
2301 }
2302
2303 Double_t AliGenMUONlib::YUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
2304 {
2305 // Upsilon y
2306 // PbPb 2.76 TeV, minimum bias 0-100 %
2307 //
2308   return YUpsilonPbPb2760ShFdummy(*px, 0) * YUpsilonPP2760(px, dummy);
2309 }
2310
2311 Double_t AliGenMUONlib::YUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
2312 {
2313 // Upsilon y
2314 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
2315 //
2316   return YUpsilonPbPb2760ShFdummy(*px, 1) * YUpsilonPP2760(px, dummy);
2317 }
2318
2319 Double_t AliGenMUONlib::YUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
2320 {
2321 // Upsilon y
2322 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
2323 //
2324   return YUpsilonPbPb2760ShFdummy(*px, 2) * YUpsilonPP2760(px, dummy);
2325 }
2326
2327 Double_t AliGenMUONlib::YUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
2328 {
2329 // Upsilon y
2330 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
2331 //
2332   return YUpsilonPbPb2760ShFdummy(*px, 3) * YUpsilonPP2760(px, dummy);
2333 }
2334
2335 Double_t AliGenMUONlib::YUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
2336 {
2337 // Upsilon y
2338 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
2339 //
2340   return YUpsilonPbPb2760ShFdummy(*px, 4) * YUpsilonPP2760(px, dummy);
2341 }
2342
2343 Double_t AliGenMUONlib::YUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
2344 {
2345 // Upsilon y
2346 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
2347 //
2348   return YUpsilonPbPb2760ShFdummy(*px, 5) * YUpsilonPP2760(px, dummy);
2349 }
2350
2351 Double_t AliGenMUONlib::YUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
2352 {
2353 // Upsilon y
2354 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
2355 //
2356   return YUpsilonPbPb2760ShFdummy(*px, 6) * YUpsilonPP2760(px, dummy);
2357 }
2358
2359 Double_t AliGenMUONlib::YUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
2360 {
2361 // Upsilon y
2362 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
2363 //
2364   return YUpsilonPbPb2760ShFdummy(*px, 7) * YUpsilonPP2760(px, dummy);
2365 }
2366
2367 Double_t AliGenMUONlib::YUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
2368 {
2369 // Upsilon y
2370 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
2371 //
2372   return YUpsilonPbPb2760ShFdummy(*px, 8) * YUpsilonPP2760(px, dummy);
2373 }
2374
2375 Double_t AliGenMUONlib::YUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
2376 {
2377 // Upsilon y
2378 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
2379 //
2380   return YUpsilonPbPb2760ShFdummy(*px, 9) * YUpsilonPP2760(px, dummy);
2381 }
2382
2383 Double_t AliGenMUONlib::YUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
2384 {
2385 // Upsilon y
2386 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
2387 //
2388   return YUpsilonPbPb2760ShFdummy(*px, 10) * YUpsilonPP2760(px, dummy);
2389 }
2390
2391 Double_t AliGenMUONlib::YUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
2392 {
2393 // Upsilon y
2394 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
2395 //
2396   return YUpsilonPbPb2760ShFdummy(*px, 11) * YUpsilonPP2760(px, dummy);
2397 }
2398
2399 Double_t AliGenMUONlib::YUpsilonPP5030dummy(Double_t px)
2400 {
2401     return AliGenMUONlib::YUpsilonPP5030(&px, (Double_t*) 0);
2402 }
2403
2404 Double_t AliGenMUONlib::YUpsilonPPb5030ShFdummy(Double_t x, Int_t n)
2405 {
2406 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2407 //
2408 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
2409 //
2410     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2411     const Double_t c[7] = {8.885e-01, 4.620e-02, 1.158e-02, 4.959e-04,-4.422e-04,-5.345e-05, 0.};
2412     Double_t y;
2413     Int_t j;
2414     y = c[j = 6];
2415     while (j > 0) y = y * x + c[--j];
2416     if(y<0) y=0;
2417     //
2418     return 1 +(y-1)*f[n];
2419 }
2420
2421 Double_t AliGenMUONlib::YUpsilonPPb5030(const Double_t *px, const Double_t */*dummy*/)
2422 {
2423 // Upsilon y
2424 // pPb 5.03 TeV, minimum bias 0-100 %
2425 //
2426   Double_t x = px[0] + 0.47;              // rapidity shift
2427   return YUpsilonPPb5030ShFdummy(x, 0) * YUpsilonPP5030dummy(x);
2428 }
2429
2430 Double_t AliGenMUONlib::YUpsilonPPb5030c1(const Double_t *px, const Double_t */*dummy*/)
2431 {
2432 // Upsilon y
2433 // pPb 5.03 TeV, 1st centrality bin 0-20 %
2434 //
2435   Double_t x = px[0] + 0.47;              // rapidity shift
2436   return YUpsilonPPb5030ShFdummy(x, 1) * YUpsilonPP5030dummy(x);
2437 }
2438
2439 Double_t AliGenMUONlib::YUpsilonPPb5030c2(const Double_t *px, const Double_t */*dummy*/)
2440 {
2441 // Upsilon y
2442 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
2443 //
2444   Double_t x = px[0] + 0.47;              // rapidity shift
2445   return YUpsilonPPb5030ShFdummy(x, 2) * YUpsilonPP5030dummy(x);
2446 }
2447
2448 Double_t AliGenMUONlib::YUpsilonPPb5030c3(const Double_t *px, const Double_t */*dummy*/)
2449 {
2450 // Upsilon y
2451 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
2452 //
2453   Double_t x = px[0] + 0.47;              // rapidity shift
2454   return YUpsilonPPb5030ShFdummy(x, 3) * YUpsilonPP5030dummy(x);
2455 }
2456
2457 Double_t AliGenMUONlib::YUpsilonPPb5030c4(const Double_t *px, const Double_t */*dummy*/)
2458 {
2459 // Upsilon y
2460 // pPb 5.03 TeV, 4th centrality bin 60-100 %
2461 //
2462   Double_t x = px[0] + 0.47;              // rapidity shift
2463   return YUpsilonPPb5030ShFdummy(x, 4) * YUpsilonPP5030dummy(x);
2464 }
2465
2466 Double_t AliGenMUONlib::YUpsilonPbP5030(const Double_t *px, const Double_t */*dummy*/)
2467 {
2468 // Upsilon y
2469 // Pbp 5.03 TeV, minimum bias 0-100 %
2470 //
2471   Double_t x = -px[0] + 0.47;              // rapidity shift
2472   return YUpsilonPPb5030ShFdummy(x, 0) * YUpsilonPP5030dummy(x);
2473 }
2474
2475 Double_t AliGenMUONlib::YUpsilonPbP5030c1(const Double_t *px, const Double_t */*dummy*/)
2476 {
2477 // Upsilon y
2478 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
2479 //
2480   Double_t x = -px[0] + 0.47;              // rapidity shift
2481   return YUpsilonPPb5030ShFdummy(x, 1) * YUpsilonPP5030dummy(x);
2482 }
2483
2484 Double_t AliGenMUONlib::YUpsilonPbP5030c2(const Double_t *px, const Double_t */*dummy*/)
2485 {
2486 // Upsilon y
2487 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
2488 //
2489   Double_t x = -px[0] + 0.47;              // rapidity shift
2490   return YUpsilonPPb5030ShFdummy(x, 2) * YUpsilonPP5030dummy(x);
2491 }
2492
2493 Double_t AliGenMUONlib::YUpsilonPbP5030c3(const Double_t *px, const Double_t */*dummy*/)
2494 {
2495 // Upsilon y
2496 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
2497 //
2498   Double_t x = -px[0] + 0.47;              // rapidity shift
2499   return YUpsilonPPb5030ShFdummy(x, 3) * YUpsilonPP5030dummy(x);
2500 }
2501
2502 Double_t AliGenMUONlib::YUpsilonPbP5030c4(const Double_t *px, const Double_t */*dummy*/)
2503 {
2504 // Upsilon y
2505 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
2506 //
2507   Double_t x = -px[0] + 0.47;              // rapidity shift
2508   return YUpsilonPPb5030ShFdummy(x, 4) * YUpsilonPP5030dummy(x);
2509 }
2510
2511 Double_t AliGenMUONlib::YUpsilonPP8800dummy(Double_t px)
2512 {
2513     return AliGenMUONlib::YUpsilonPP8800(&px, (Double_t*) 0);
2514 }
2515
2516 Double_t AliGenMUONlib::YUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
2517 {
2518 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2519 //
2520 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
2521 //
2522     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2523     const Double_t c[7] = {8.6581e-01, 4.6111e-02, 7.6911e-03, 8.7313e-04,
2524                            -1.4700e-04,-5.0975e-05,-3.5718e-06}; 
2525     Double_t y;
2526     Int_t j;
2527     y = c[j = 6];
2528     while (j > 0) y = y * x + c[--j];
2529     if(y<0) y=0;
2530     //
2531     return 1 +(y-1)*f[n];
2532 }
2533
2534 Double_t AliGenMUONlib::YUpsilonPPb8800(const Double_t *px, const Double_t */*dummy*/)
2535 {
2536 // Upsilon y
2537 // pPb 8.8 TeV, minimum bias 0-100 %
2538 //
2539   Double_t x = px[0] + 0.47;              // rapidity shift
2540   return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2541 }
2542
2543 Double_t AliGenMUONlib::YUpsilonPPb8800c1(const Double_t *px, const Double_t */*dummy*/)
2544 {
2545 // Upsilon y
2546 // pPb 8.8 TeV, 1st centrality bin 0-20 %
2547 //
2548   Double_t x = px[0] + 0.47;              // rapidity shift
2549   return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2550 }
2551
2552 Double_t AliGenMUONlib::YUpsilonPPb8800c2(const Double_t *px, const Double_t */*dummy*/)
2553 {
2554 // Upsilon y
2555 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
2556 //
2557   Double_t x = px[0] + 0.47;              // rapidity shift
2558   return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2559 }
2560
2561 Double_t AliGenMUONlib::YUpsilonPPb8800c3(const Double_t *px, const Double_t */*dummy*/)
2562 {
2563 // Upsilon y
2564 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
2565 //
2566   Double_t x = px[0] + 0.47;              // rapidity shift
2567   return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2568 }
2569
2570 Double_t AliGenMUONlib::YUpsilonPPb8800c4(const Double_t *px, const Double_t */*dummy*/)
2571 {
2572 // Upsilon y
2573 // pPb 8.8 TeV, 4th centrality bin 60-100 %
2574 //
2575   Double_t x = px[0] + 0.47;              // rapidity shift
2576   return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2577 }
2578
2579 Double_t AliGenMUONlib::YUpsilonPbP8800(const Double_t *px, const Double_t */*dummy*/)
2580 {
2581 // Upsilon y
2582 // Pbp 8.8 TeV, minimum bias 0-100 %
2583 //
2584   Double_t x = -px[0] + 0.47;              // rapidity shift
2585   return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2586 }
2587
2588 Double_t AliGenMUONlib::YUpsilonPbP8800c1(const Double_t *px, const Double_t */*dummy*/)
2589 {
2590 // Upsilon y
2591 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
2592 //
2593   Double_t x = -px[0] + 0.47;              // rapidity shift
2594   return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2595 }
2596
2597 Double_t AliGenMUONlib::YUpsilonPbP8800c2(const Double_t *px, const Double_t */*dummy*/)
2598 {
2599 // Upsilon y
2600 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
2601 //
2602   Double_t x = -px[0] + 0.47;              // rapidity shift
2603   return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2604 }
2605
2606 Double_t AliGenMUONlib::YUpsilonPbP8800c3(const Double_t *px, const Double_t */*dummy*/)
2607 {
2608 // Upsilon y
2609 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
2610 //
2611   Double_t x = -px[0] + 0.47;              // rapidity shift
2612   return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2613 }
2614
2615 Double_t AliGenMUONlib::YUpsilonPbP8800c4(const Double_t *px, const Double_t */*dummy*/)
2616 {
2617 // Upsilon y
2618 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
2619 //
2620   Double_t x = -px[0] + 0.47;              // rapidity shift
2621   return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2622 }
2623
2624 Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
2625 {
2626 // Upsilon y
2627   const Double_t ky0 = 3.;
2628   const Double_t kb=1.;
2629   Double_t yu;
2630   Double_t y=TMath::Abs(*py);
2631   //
2632   if (y < ky0)
2633     yu=kb;
2634   else
2635     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2636   return yu;
2637 }
2638
2639
2640 Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2641 {
2642
2643 //
2644 // Upsilon y
2645 //
2646 //
2647 // R. Vogt 2002
2648 // PbPb 5.5 TeV
2649 // MRST HO
2650 // mc = 1.4 GeV, pt-kick 1 GeV
2651 //
2652
2653     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
2654                      -2.99753e-09, 1.28895e-05};
2655     Double_t x = TMath::Abs(px[0]);
2656     if (x > 5.55) return 0.;
2657     Int_t j;
2658     Double_t y = c[j = 6];
2659     while (j > 0) y  = y * x +c[--j];
2660     return y;
2661 }
2662
2663 Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
2664 {
2665     // Upsilon y
2666     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
2667     
2668 }
2669
2670 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
2671 {
2672     // Upsilon y
2673     return AliGenMUONlib::YUpsilonPP(px, dummy);
2674     
2675 }
2676
2677 Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
2678 {
2679     // Upsilon y
2680     return 1.;
2681     
2682 }
2683
2684 Double_t AliGenMUONlib::YUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2685 {
2686 // Upsilon y
2687 //
2688 // pp 10 TeV
2689 // scaled from YUpsilonPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
2690 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
2691 //
2692     Double_t c[4] = {1., -2.17877e-02, -6.52830e-04, 1.40578e-05};
2693     Double_t x = TMath::Abs(px[0]);
2694     if (x > 6.1) return 0.;
2695     Int_t j;
2696     Double_t y = c[j = 3];
2697     while (j > 0) y  = y * x*x +c[--j];
2698     return y;
2699 }
2700
2701 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2702 {
2703 // Upsilon y
2704 //
2705 // pp 8.8 TeV
2706 // rescaling of YUpsilonPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD 
2707 //
2708     Double_t c[4] = {1., -2.37621e-02, -6.29610e-04, 1.47976e-05};
2709     Double_t x = TMath::Abs(px[0]);
2710     if (x > 6.1) return 0.;
2711     Int_t j;
2712     Double_t y = c[j = 3];
2713     while (j > 0) y  = y * x*x +c[--j];
2714     return y;
2715 }
2716
2717 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9dummy(Double_t px)
2718 {
2719     return AliGenMUONlib::YUpsilonCDFscaledPP9(&px, (Double_t*) 0);
2720 }
2721
2722 Double_t AliGenMUONlib::YUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2723 {
2724 // Upsilon y
2725 //
2726 // pp 7 TeV
2727 // scaled from YUpsilonPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
2728 //
2729     Double_t c[4] = {1., -2.61009e-02, -6.83937e-04, 1.78451e-05};
2730     Double_t x = TMath::Abs(px[0]);
2731     if (x > 6.0) return 0.;
2732     Int_t j;
2733     Double_t y = c[j = 3];
2734     while (j > 0) y  = y * x*x +c[--j];
2735     return y;
2736 }
2737
2738 Double_t AliGenMUONlib::YUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2739 {
2740 // Upsilon y
2741 //
2742 // pp 3.94 TeV
2743 // rescaling of YUpsilonPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
2744 //
2745     Double_t c[4] = {1., -3.91924e-02, -4.26184e-04, 2.10914e-05};
2746     Double_t x = TMath::Abs(px[0]);
2747     if (x > 5.7) return 0.;
2748     Int_t j;
2749     Double_t y = c[j = 3];
2750     while (j > 0) y  = y * x*x +c[--j];
2751
2752     return y;
2753 }
2754
2755 Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2756 {
2757
2758 //
2759 // Upsilon y
2760 //
2761 //
2762 // R. Vogt 2002
2763 // p p  14. TeV
2764 // MRST HO
2765 // mc = 1.4 GeV, pt-kick 1 GeV
2766 //
2767     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, 
2768                      -6.20539e-10, 1.29943e-05};
2769     Double_t x = TMath::Abs(px[0]);
2770     if (x > 6.2) return 0.;
2771     Int_t j;
2772     Double_t y = c[j = 6];
2773     while (j > 0) y  = y * x +c[--j];
2774     return y;
2775 }
2776
2777 Double_t AliGenMUONlib::YUpsilonCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
2778 {
2779 // Upsilon y
2780 //
2781 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2782 //
2783     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2784                      -4.67538e-05,-2.11683e-06}; 
2785     Double_t y;
2786     Double_t x = px[0] + 0.47;              // rapidity shift
2787     Int_t j;
2788     y = c[j = 6];
2789     while (j > 0) y = y * x + c[--j];
2790     if(y<0) y=0;
2791
2792     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2793 }
2794
2795 Double_t AliGenMUONlib::YUpsilonCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
2796 {
2797 // Upsilon y
2798 //
2799 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2800 //
2801     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2802                      -4.67538e-05,-2.11683e-06}; 
2803     Double_t y;
2804     Double_t x = -px[0] + 0.47;              // rapidity shift
2805     Int_t j;
2806     y = c[j = 6];
2807     while (j > 0) y = y * x + c[--j];
2808     if(y<0) y=0;
2809
2810     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2811 }
2812
2813 Double_t AliGenMUONlib::YUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2814 {
2815 // Upsilon y
2816 //
2817 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2818 //
2819     Double_t c[4] = {8.27837e-01, 1.70115e-02, -1.26046e-03, 1.52091e-05}; 
2820     Double_t x = px[0]*px[0];
2821     Double_t y;
2822     Int_t j;
2823     y = c[j = 3];
2824     while (j > 0) y  = y * x + c[--j];
2825     if(y<0) y=0;
2826
2827     return y * AliGenMUONlib::YUpsilonCDFscaledPP4(px,dummy);
2828 }
2829
2830
2831 //                 particle composition
2832 //
2833 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
2834 {
2835 // y composition
2836     return 553;
2837 }
2838 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
2839 {
2840 // y composition
2841     return 100553;
2842 }
2843 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
2844 {
2845 // y composition
2846     return 200553;
2847 }
2848 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
2849 {
2850 // y composition
2851   // Using the LHCb pp data at 7 TeV: CERN-PH-EP-2012-051
2852   // (L. Manceau, S. Grigoryan)
2853   Int_t ip;
2854   Float_t r = gRandom->Rndm();  
2855   if (r < 0.687) {
2856     //  if (r < 0.712) {
2857     ip = 553;
2858   } else if (r < 0.903) {
2859     //  } else if (r < 0.896) {
2860     ip = 100553;
2861   } else {
2862     ip = 200553;
2863   }
2864   return ip;
2865 }
2866
2867
2868 //
2869 //                        Phi
2870 //
2871 //
2872 //    pt-distribution (by scaling of pion distribution)
2873 //____________________________________________________________
2874 Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
2875 {
2876 // Phi pT
2877   return PtScal(*px,7);
2878 }
2879 //    y-distribution
2880 Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
2881 {
2882 // Phi y
2883     Double_t *dum=0;
2884     return YJpsi(px,dum);
2885 }
2886 //                 particle composition
2887 //
2888 Int_t AliGenMUONlib::IpPhi(TRandom *)
2889 {
2890 // Phi composition
2891     return 333;
2892 }
2893
2894 //
2895 //                        omega
2896 //
2897 //
2898 //    pt-distribution (by scaling of pion distribution)
2899 //____________________________________________________________
2900 Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
2901 {
2902 // Omega pT
2903   return PtScal(*px,5);
2904 }
2905 //    y-distribution
2906 Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
2907 {
2908 // Omega y
2909     Double_t *dum=0;
2910     return YJpsi(px,dum);
2911 }
2912 //                 particle composition
2913 //
2914 Int_t AliGenMUONlib::IpOmega(TRandom *)
2915 {
2916 // Omega composition
2917     return 223;
2918 }
2919
2920
2921 //
2922 //                        Eta
2923 //
2924 //
2925 //    pt-distribution (by scaling of pion distribution)
2926 //____________________________________________________________
2927 Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
2928 {
2929 // Eta pT
2930   return PtScal(*px,3);
2931 }
2932 //    y-distribution
2933 Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
2934 {
2935 // Eta y
2936     Double_t *dum=0;
2937     return YJpsi(px,dum);
2938 }
2939 //                 particle composition
2940 //
2941 Int_t AliGenMUONlib::IpEta(TRandom *)
2942 {
2943 // Eta composition
2944     return 221;
2945 }
2946
2947 //
2948 //                        Charm
2949 //
2950 //
2951 //                    pt-distribution
2952 //____________________________________________________________
2953 Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
2954 {
2955 // Charm pT
2956   const Double_t kpt0 = 2.25;
2957   const Double_t kxn  = 3.17;
2958   Double_t x=*px;
2959   //
2960   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2961   return x/TMath::Power(pass1,kxn);
2962 }
2963
2964 Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
2965 {
2966 // Charm pT
2967   const Double_t kpt0 = 2.12;
2968   const Double_t kxn  = 2.78;
2969   Double_t x=*px;
2970   //
2971   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2972   return x/TMath::Power(pass1,kxn);
2973 }
2974 Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2975 {
2976 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
2977 // PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
2978 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
2979 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
2980 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
2981 // calculations for the following inputs: 
2982 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11 
2983 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV 
2984 // for j=0,1 & 2 respectively; 
2985 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
2986 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 
2987 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set 
2988 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
2989 // June 2008, Smbat.Grigoryan@cern.ch
2990
2991 // Charm pT
2992 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
2993 // for pp collisions at 14 TeV with one c-cbar pair per event.
2994 // Corresponding NLO total cross section is 5.68 mb
2995
2996
2997   const Double_t kpt0 = 2.2930;
2998   const Double_t kxn  = 3.1196;
2999   Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
3000   Double_t x=*px;
3001   //
3002   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3003   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3004 }
3005 Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3006 {
3007 // Charm pT
3008 // Corresponding NLO total cross section is 6.06 mb
3009   const Double_t kpt0 = 2.8669;
3010   const Double_t kxn  = 3.1044;
3011   Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
3012   Double_t x=*px;
3013   //
3014   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3015   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3016 }
3017 Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3018 {
3019 // Charm pT
3020 // Corresponding NLO total cross section is 6.06 mb
3021   const Double_t kpt0 = 1.8361;
3022   const Double_t kxn  = 3.2966;
3023   Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
3024   Double_t x=*px;
3025   //
3026   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3027   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3028 }
3029 Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3030 {
3031 // Charm pT
3032 // Corresponding NLO total cross section is 7.69 mb
3033   const Double_t kpt0 = 2.1280;
3034   const Double_t kxn  = 3.1397;
3035   Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
3036   Double_t x=*px;
3037   //
3038   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3039   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3040 }
3041 Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3042 {
3043 // Charm pT
3044 // Corresponding NLO total cross section is 4.81 mb
3045   const Double_t kpt0 = 2.4579;
3046   const Double_t kxn  = 3.1095;
3047   Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
3048   Double_t x=*px;
3049   //
3050   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3051   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3052 }
3053 Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3054 {
3055 // Charm pT
3056 // Corresponding NLO total cross section is 14.09 mb
3057   const Double_t kpt0 = 2.1272;
3058   const Double_t kxn  = 3.1904;
3059   Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
3060   Double_t x=*px;
3061   //
3062   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3063   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3064 }
3065 Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3066 {
3067 // Charm pT
3068 // Corresponding NLO total cross section is 1.52 mb
3069   const Double_t kpt0 = 2.8159;
3070   const Double_t kxn  = 3.0857;
3071   Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
3072   Double_t x=*px;
3073   //
3074   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3075   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3076 }
3077 Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3078 {
3079 // Charm pT
3080 // Corresponding NLO total cross section is 3.67 mb
3081   const Double_t kpt0 = 2.7297;
3082   const Double_t kxn  = 3.3019;
3083   Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
3084   Double_t x=*px;
3085   //
3086   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3087   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3088 }
3089 Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3090 {
3091 // Charm pT
3092 // Corresponding NLO total cross section is 3.38 mb
3093   const Double_t kpt0 = 2.3894;
3094   const Double_t kxn  = 3.1075;
3095   Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
3096   Double_t x=*px;
3097   //
3098   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3099   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3100 }
3101 Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3102 {
3103 // Charm pT
3104 // Corresponding NLO total cross section is 10.37 mb
3105   const Double_t kpt0 = 2.0187;
3106   const Double_t kxn  = 3.3011;
3107   Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
3108   Double_t x=*px;
3109   //
3110   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3111   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3112 }
3113 Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3114 {
3115 // Charm pT
3116 // Corresponding NLO total cross section is 7.22 mb
3117   const Double_t kpt0 = 2.1089;
3118   const Double_t kxn  = 3.1848;
3119   Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
3120   Double_t x=*px;
3121   //
3122   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3123   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3124 }
3125
3126 //                  y-distribution
3127 Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
3128 {
3129 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
3130 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3131 // shadowing + kt broadening 
3132
3133     Double_t x=px[0];
3134     Double_t c[2]={-2.42985e-03,-2.31001e-04};
3135     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
3136     Double_t ycharm;
3137     
3138     if (TMath::Abs(x)>8) {
3139       ycharm=0.;
3140     }
3141     else {
3142       ycharm=TMath::Power(y,3);
3143     }
3144     
3145     return ycharm;
3146 }
3147 Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3148 {
3149 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3150 // YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3151 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3152 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
3153 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
3154 // calculations for the following inputs: 
3155 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11 
3156 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV 
3157 // for j=0,1 & 2 respectively; 
3158 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
3159 // with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for 
3160 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3161 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3162 // June 2008, Smbat.Grigoryan@cern.ch
3163
3164 // Charm y
3165 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
3166 // for pp collisions at 14 TeV with one c-cbar pair per event.
3167 // Corresponding NLO total cross section is 5.68 mb
3168
3169     Double_t x=px[0];
3170     Double_t c[2]={7.0909e-03,6.1967e-05};
3171     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3172     Double_t ycharm;
3173     
3174     if (TMath::Abs(x)>9) {
3175       ycharm=0.;
3176     }
3177     else {
3178       ycharm=TMath::Power(y,3);
3179     }
3180     
3181     return ycharm;
3182 }
3183 Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3184 {
3185 // Charm y
3186 // Corresponding NLO total cross section is 6.06 mb
3187     Double_t x=px[0];
3188     Double_t c[2]={6.9707e-03,6.0971e-05};
3189     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3190     Double_t ycharm;
3191     
3192     if (TMath::Abs(x)>9) {
3193       ycharm=0.;
3194     }
3195     else {
3196       ycharm=TMath::Power(y,3);
3197     }
3198     
3199     return ycharm;
3200 }
3201 Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3202 {
3203 // Charm y
3204 // Corresponding NLO total cross section is 6.06 mb
3205     Double_t x=px[0];
3206     Double_t c[2]={7.1687e-03,6.5303e-05};
3207     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3208     Double_t ycharm;
3209     
3210     if (TMath::Abs(x)>9) {
3211       ycharm=0.;
3212     }
3213     else {
3214       ycharm=TMath::Power(y,3);
3215     }
3216     
3217     return ycharm;
3218 }
3219 Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3220 {
3221 // Charm y
3222 // Corresponding NLO total cross section is 7.69 mb
3223     Double_t x=px[0];
3224     Double_t c[2]={5.9090e-03,7.1854e-05};
3225     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3226     Double_t ycharm;
3227     
3228     if (TMath::Abs(x)>9) {
3229       ycharm=0.;
3230     }
3231     else {
3232       ycharm=TMath::Power(y,3);
3233     }
3234     
3235     return ycharm;
3236 }
3237 Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3238 {
3239 // Charm y
3240 // Corresponding NLO total cross section is 4.81 mb
3241     Double_t x=px[0];
3242     Double_t c[2]={8.0882e-03,5.5872e-05};
3243     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3244     Double_t ycharm;
3245     
3246     if (TMath::Abs(x)>9) {
3247       ycharm=0.;
3248     }
3249     else {
3250       ycharm=TMath::Power(y,3);
3251     }
3252     
3253     return ycharm;
3254 }
3255 Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3256 {
3257 // Charm y
3258 // Corresponding NLO total cross section is 14.09 mb
3259     Double_t x=px[0];
3260     Double_t c[2]={7.2520e-03,6.2691e-05};
3261     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3262     Double_t ycharm;
3263     
3264     if (TMath::Abs(x)>9) {
3265       ycharm=0.;
3266     }
3267     else {
3268       ycharm=TMath::Power(y,3);
3269     }
3270     
3271     return ycharm;
3272 }
3273 Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3274 {
3275 // Charm y
3276 // Corresponding NLO total cross section is 1.52 mb
3277     Double_t x=px[0];
3278     Double_t c[2]={1.1040e-04,1.4498e-04};
3279     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3280     Double_t ycharm;
3281     
3282     if (TMath::Abs(x)>9) {
3283       ycharm=0.;
3284     }
3285     else {
3286       ycharm=TMath::Power(y,3);
3287     }
3288     
3289     return ycharm;
3290 }
3291 Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3292 {
3293 // Charm y
3294 // Corresponding NLO total cross section is 3.67 mb
3295     Double_t x=px[0];
3296     Double_t c[2]={-3.1328e-03,1.8270e-04};
3297     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3298     Double_t ycharm;
3299     
3300     if (TMath::Abs(x)>9) {
3301       ycharm=0.;
3302     }
3303     else {
3304       ycharm=TMath::Power(y,3);
3305     }
3306     
3307     return ycharm;
3308 }
3309 Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3310 {
3311 // Charm y
3312 // Corresponding NLO total cross section is 3.38 mb
3313     Double_t x=px[0];
3314     Double_t c[2]={7.0865e-03,6.2532e-05};
3315     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3316     Double_t ycharm;
3317     
3318     if (TMath::Abs(x)>9) {
3319       ycharm=0.;
3320     }
3321     else {
3322       ycharm=TMath::Power(y,3);
3323     }
3324     
3325     return ycharm;
3326 }
3327 Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3328 {
3329 // Charm y
3330 // Corresponding NLO total cross section is 10.37 mb
3331     Double_t x=px[0];
3332     Double_t c[2]={7.7070e-03,5.3533e-05};
3333     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3334     Double_t ycharm;
3335     
3336     if (TMath::Abs(x)>9) {
3337       ycharm=0.;
3338     }
3339     else {
3340       ycharm=TMath::Power(y,3);
3341     }
3342     
3343     return ycharm;
3344 }
3345 Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3346 {
3347 // Charm y
3348 // Corresponding NLO total cross section is 7.22 mb
3349     Double_t x=px[0];
3350     Double_t c[2]={7.9195e-03,5.3823e-05};
3351     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3352     Double_t ycharm;
3353     
3354     if (TMath::Abs(x)>9) {
3355       ycharm=0.;
3356     }
3357     else {
3358       ycharm=TMath::Power(y,3);
3359     }
3360     
3361     return ycharm;
3362 }
3363
3364
3365 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
3366 {  
3367 // Charm composition
3368     Float_t random;
3369     Int_t ip;
3370 //    411,421,431,4122
3371     random = ran->Rndm();
3372 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
3373 //  >>>>> cf. tab 4 p 11
3374   
3375     if (random < 0.30) {                       
3376         ip=421;
3377     } else if (random < 0.60) {
3378         ip=-421;
3379     } else if (random < 0.70) {
3380         ip=411;
3381     } else if (random < 0.80) {
3382         ip=-411;
3383     } else if (random < 0.86) {
3384         ip=431;
3385     } else if (random < 0.92) {
3386         ip=-431;        
3387     } else if (random < 0.96) {
3388         ip=4122;
3389     } else {
3390         ip=-4122;
3391     }
3392     
3393     return ip;
3394 }
3395
3396 //
3397 //                        Beauty
3398 //
3399 //
3400 //                    pt-distribution
3401 //____________________________________________________________
3402 Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
3403 {
3404 // Beauty pT
3405   const Double_t kpt0 = 6.53;
3406   const Double_t kxn  = 3.59;
3407   Double_t x=*px;
3408   //
3409   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3410   return x/TMath::Power(pass1,kxn);
3411 }
3412
3413 Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
3414 {
3415 // Beauty pT
3416   const Double_t kpt0 = 6.14;
3417   const Double_t kxn  = 2.93;
3418   Double_t x=*px;
3419   //
3420   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3421   return x/TMath::Power(pass1,kxn);
3422 }
3423 Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3424 {
3425 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3426 // PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
3427 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3428 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
3429 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
3430 // calculations for the following inputs: 
3431 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004 
3432 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV 
3433 // for j=0,1 & 2 respectively; 
3434 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
3435 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for 
3436 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3437 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3438 // June 2008, Smbat.Grigoryan@cern.ch
3439
3440 // Beauty pT
3441 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3442 // for pp collisions at 14 TeV with one b-bbar pair per event.
3443 // Corresponding NLO total cross section is 0.494 mb
3444
3445   const Double_t kpt0 = 8.0575;
3446   const Double_t kxn  = 3.1921;
3447   Double_t x=*px;
3448   //
3449   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3450   return x/TMath::Power(pass1,kxn);
3451 }
3452 Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3453 {
3454 // Beauty pT
3455 // Corresponding NLO total cross section is 0.445 mb
3456   const Double_t kpt0 = 8.6239;
3457   const Double_t kxn  = 3.2911;
3458   Double_t x=*px;
3459   //
3460   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3461   return x/TMath::Power(pass1,kxn);
3462 }
3463 Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3464 {
3465 // Beauty pT
3466 // Corresponding NLO total cross section is 0.445 mb
3467   const Double_t kpt0 = 7.3367;
3468   const Double_t kxn  = 3.0692;
3469   Double_t x=*px;
3470   //
3471   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3472   return x/TMath::Power(pass1,kxn);
3473 }
3474 Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3475 {
3476 // Beauty pT
3477 // Corresponding NLO total cross section is 0.518 mb
3478   const Double_t kpt0 = 7.6409;
3479   const Double_t kxn  = 3.1364;
3480   Double_t x=*px;
3481   //
3482   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3483   return x/TMath::Power(pass1,kxn);
3484 }
3485 Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3486 {
3487 // Beauty pT
3488 // Corresponding NLO total cross section is 0.384 mb
3489   const Double_t kpt0 = 8.4948;
3490   const Double_t kxn  = 3.2546;
3491   Double_t x=*px;
3492   //
3493   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3494   return x/TMath::Power(pass1,kxn);
3495 }
3496 Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3497 {
3498 // Beauty pT
3499 // Corresponding NLO total cross section is 0.648 mb
3500   const Double_t kpt0 = 7.6631;
3501   const Double_t kxn  = 3.1621;
3502   Double_t x=*px;
3503   //
3504   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3505   return x/TMath::Power(pass1,kxn);
3506 }
3507 Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3508 {
3509 // Beauty pT
3510 // Corresponding NLO total cross section is 0.294 mb
3511   const Double_t kpt0 = 8.7245;
3512   const Double_t kxn  = 3.2213;
3513   Double_t x=*px;
3514   //
3515   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3516   return x/TMath::Power(pass1,kxn);
3517 }
3518 Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3519 {
3520 // Beauty pT
3521 // Corresponding NLO total cross section is 0.475 mb
3522   const Double_t kpt0 = 8.5296;
3523   const Double_t kxn  = 3.2187;
3524   Double_t x=*px;
3525   //
3526   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3527   return x/TMath::Power(pass1,kxn);
3528 }
3529 Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3530 {
3531 // Beauty pT
3532 // Corresponding NLO total cross section is 0.324 mb
3533   const Double_t kpt0 = 7.9440;
3534   const Double_t kxn  = 3.1614;
3535   Double_t x=*px;
3536   //
3537   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3538   return x/TMath::Power(pass1,kxn);
3539 }
3540 Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3541 {
3542 // Beauty pT
3543 // Corresponding NLO total cross section is 0.536 mb
3544   const Double_t kpt0 = 8.2408;
3545   const Double_t kxn  = 3.3029;
3546   Double_t x=*px;
3547   //
3548   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3549   return x/TMath::Power(pass1,kxn);
3550 }
3551 Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3552 {
3553 // Beauty pT
3554 // Corresponding NLO total cross section is 0.420 mb
3555   const Double_t kpt0 = 7.8041;
3556   const Double_t kxn  = 3.2094;
3557   Double_t x=*px;
3558   //
3559   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3560   return x/TMath::Power(pass1,kxn);
3561 }
3562
3563 //                     y-distribution
3564 Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
3565 {
3566 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
3567 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3568 // shadowing + kt broadening 
3569
3570     Double_t x=px[0];
3571     Double_t c[2]={-1.27590e-02,-2.42731e-04};
3572     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
3573     Double_t ybeauty;
3574     
3575     if (TMath::Abs(x)>6) {
3576       ybeauty=0.;
3577     }
3578     else {
3579       ybeauty=TMath::Power(y,3);
3580     }
3581     
3582     return ybeauty;
3583 }
3584 Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3585 {
3586 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3587 // YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3588 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3589 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
3590 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
3591 // calculations for the following inputs: 
3592 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004 
3593 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV 
3594 // for j=0,1 & 2 respectively; 
3595 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
3596 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 
3597 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set 
3598 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3599 // June 2008, Smbat.Grigoryan@cern.ch
3600
3601 // Beauty y
3602 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3603 // for pp collisions at 14 TeV with one b-bbar pair per event.
3604 // Corresponding NLO total cross section is 0.494 mb
3605
3606
3607     Double_t x=px[0];
3608     Double_t c[2]={1.2350e-02,9.2667e-05};
3609     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3610     Double_t ybeauty;
3611     
3612     if (TMath::Abs(x)>7.6) {
3613       ybeauty=0.;
3614     }
3615     else {
3616       ybeauty=TMath::Power(y,3);
3617     }
3618     
3619     return ybeauty;
3620 }
3621 Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3622 {
3623 // Beauty y
3624 // Corresponding NLO total cross section is 0.445 mb
3625     Double_t x=px[0];
3626     Double_t c[2]={1.2292e-02,9.1847e-05};
3627     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3628     Double_t ybeauty;
3629     
3630     if (TMath::Abs(x)>7.6) {
3631       ybeauty=0.;
3632     }
3633     else {
3634       ybeauty=TMath::Power(y,3);
3635     }
3636     
3637     return ybeauty;
3638 }
3639 Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3640 {
3641 // Beauty y
3642 // Corresponding NLO total cross section is 0.445 mb
3643     Double_t x=px[0];
3644     Double_t c[2]={1.2436e-02,9.3709e-05};
3645     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3646     Double_t ybeauty;
3647     
3648     if (TMath::Abs(x)>7.6) {
3649       ybeauty=0.;
3650     }
3651     else {
3652       ybeauty=TMath::Power(y,3);
3653     }
3654     
3655     return ybeauty;
3656 }
3657 Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3658 {
3659 // Beauty y
3660 // Corresponding NLO total cross section is 0.518 mb
3661     Double_t x=px[0];
3662     Double_t c[2]={1.1714e-02,1.0068e-04};
3663     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3664     Double_t ybeauty;
3665     
3666     if (TMath::Abs(x)>7.6) {
3667       ybeauty=0.;
3668     }
3669     else {
3670       ybeauty=TMath::Power(y,3);
3671     }
3672     
3673     return ybeauty;
3674 }
3675 Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3676 {
3677 // Beauty y
3678 // Corresponding NLO total cross section is 0.384 mb
3679     Double_t x=px[0];
3680     Double_t c[2]={1.2944e-02,8.5500e-05};
3681     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3682     Double_t ybeauty;
3683     
3684     if (TMath::Abs(x)>7.6) {
3685       ybeauty=0.;
3686     }
3687     else {
3688       ybeauty=TMath::Power(y,3);
3689     }
3690     
3691     return ybeauty;
3692 }
3693 Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3694 {
3695 // Beauty y
3696 // Corresponding NLO total cross section is 0.648 mb
3697     Double_t x=px[0];
3698     Double_t c[2]={1.2455e-02,9.2713e-05};
3699     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3700     Double_t ybeauty;
3701     
3702     if (TMath::Abs(x)>7.6) {
3703       ybeauty=0.;
3704     }
3705     else {
3706       ybeauty=TMath::Power(y,3);
3707     }
3708     
3709     return ybeauty;
3710 }
3711 Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3712 {
3713 // Beauty y
3714 // Corresponding NLO total cross section is 0.294 mb
3715     Double_t x=px[0];
3716     Double_t c[2]={1.0897e-02,1.1878e-04};
3717     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3718     Double_t ybeauty;
3719     
3720     if (TMath::Abs(x)>7.6) {
3721       ybeauty=0.;
3722     }
3723     else {
3724       ybeauty=TMath::Power(y,3);
3725     }
3726     
3727     return ybeauty;
3728 }
3729 Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3730 {
3731 // Beauty y
3732 // Corresponding NLO total cross section is 0.475 mb
3733     Double_t x=px[0];
3734     Double_t c[2]={1.0912e-02,1.1858e-04};
3735     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3736     Double_t ybeauty;
3737     
3738     if (TMath::Abs(x)>7.6) {
3739       ybeauty=0.;
3740     }
3741     else {
3742       ybeauty=TMath::Power(y,3);
3743     }
3744     
3745     return ybeauty;
3746 }
3747 Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3748 {
3749 // Beauty y
3750 // Corresponding NLO total cross section is 0.324 mb
3751     Double_t x=px[0];
3752     Double_t c[2]={1.2378e-02,9.2490e-05};
3753     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3754     Double_t ybeauty;
3755     
3756     if (TMath::Abs(x)>7.6) {
3757       ybeauty=0.;
3758     }
3759     else {
3760       ybeauty=TMath::Power(y,3);
3761     }
3762     
3763     return ybeauty;
3764 }
3765 Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3766 {
3767 // Beauty y
3768 // Corresponding NLO total cross section is 0.536 mb
3769     Double_t x=px[0];
3770     Double_t c[2]={1.2886e-02,8.2912e-05};
3771     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3772     Double_t ybeauty;
3773     
3774     if (TMath::Abs(x)>7.6) {
3775       ybeauty=0.;
3776     }
3777     else {
3778       ybeauty=TMath::Power(y,3);
3779     }
3780     
3781     return ybeauty;
3782 }
3783 Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3784 {
3785 // Beauty y
3786 // Corresponding NLO total cross section is 0.420 mb
3787     Double_t x=px[0];
3788     Double_t c[2]={1.3106e-02,8.0115e-05};
3789     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3790     Double_t ybeauty;
3791     
3792     if (TMath::Abs(x)>7.6) {
3793       ybeauty=0.;
3794     }
3795     else {
3796       ybeauty=TMath::Power(y,3);
3797     }
3798     
3799     return ybeauty;
3800 }
3801
3802 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
3803 {  
3804 // Beauty Composition
3805     Float_t random;
3806     Int_t ip;
3807     random = ran->Rndm(); 
3808     
3809 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
3810 //  >>>>> cf. tab 4 p 11
3811     
3812  if (random < 0.20) {                       
3813         ip=511;
3814     } else if (random < 0.40) {
3815         ip=-511;
3816     } else if (random < 0.605) {
3817         ip=521;
3818     } else if (random < 0.81) {
3819         ip=-521;
3820     } else if (random < 0.87) {
3821         ip=531;
3822     } else if (random < 0.93) {
3823         ip=-531;