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