Setter for calling Generate n-times
[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::PtUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
1616 {
1617 // Upsilon pT
1618 // pp 2.76 TeV
1619 //
1620   return PtUpsilonPPdummy(*px,2760);
1621 }
1622
1623 Double_t AliGenMUONlib::PtUpsilonPP4400(const Double_t *px, const Double_t */*dummy*/)
1624 {
1625 // Upsilon pT
1626 // pp 4.4 TeV
1627 //
1628   return PtUpsilonPPdummy(*px,4400);
1629 }
1630
1631 Double_t AliGenMUONlib::PtUpsilonPP5030(const Double_t *px, const Double_t */*dummy*/)
1632 {
1633 // Upsilon pT
1634 // pp 5.03 TeV
1635 //
1636   return PtUpsilonPPdummy(*px,5030);
1637 }
1638
1639 Double_t AliGenMUONlib::PtUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
1640 {
1641 // Upsilon pT
1642 // pp 8.8 TeV
1643 //
1644   return PtUpsilonPPdummy(*px,8800);
1645 }
1646
1647 Double_t AliGenMUONlib::PtUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
1648 {
1649 // Usilon shadowing factor vs pT for PbPb min. bias and 11 centr. bins (in 2.5<y<4)
1650 //
1651 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
1652 //
1653   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
1654                            0.428, 0.317, 0.231, 0.156};
1655   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
1656                            0.106, 0.041, 0.013, 0.002};
1657   const Double_t c1[7] = {1.9089e+00, 1.2969e-03, 8.9786e-05,-5.3062e-06,
1658                           -1.0046e-06,6.1446e-08, 1.0885e-09};
1659   const Double_t c2[7] = {8.8423e-01,-8.7488e-05, 5.9857e-04,-5.7959e-05, 
1660                           2.0059e-06,-2.7343e-08, 6.6053e-10};
1661   Double_t y1, y2;
1662   Int_t j;
1663   y1 = c1[j = 6]; y2 = c2[6];
1664   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
1665   
1666   y1 /= 1.+c1[6]*TMath::Power(x,6);
1667   y2 /= 1.+c2[6]*TMath::Power(x,6);
1668   //  
1669   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
1670   if(y1<0) y1=0;
1671   return y1;
1672 }
1673
1674 Double_t AliGenMUONlib::PtUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
1675 {
1676 // Upsilon pT
1677 // PbPb 2.76 TeV, minimum bias 0-100 %
1678 //
1679   return PtUpsilonPbPb2760ShFdummy(*px, 0) * PtUpsilonPP2760(px, dummy);
1680 }
1681
1682 Double_t AliGenMUONlib::PtUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
1683 {
1684 // Upsilon pT
1685 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
1686 //
1687   return PtUpsilonPbPb2760ShFdummy(*px, 1) * PtUpsilonPP2760(px, dummy);
1688 }
1689
1690 Double_t AliGenMUONlib::PtUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
1691 {
1692 // Upsilon pT
1693 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
1694 //
1695   return PtUpsilonPbPb2760ShFdummy(*px, 2) * PtUpsilonPP2760(px, dummy);
1696 }
1697
1698 Double_t AliGenMUONlib::PtUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
1699 {
1700 // Upsilon pT
1701 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
1702 //
1703   return PtUpsilonPbPb2760ShFdummy(*px, 3) * PtUpsilonPP2760(px, dummy);
1704 }
1705
1706 Double_t AliGenMUONlib::PtUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
1707 {
1708 // Upsilon pT
1709 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
1710 //
1711   return PtUpsilonPbPb2760ShFdummy(*px, 4) * PtUpsilonPP2760(px, dummy);
1712 }
1713
1714 Double_t AliGenMUONlib::PtUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
1715 {
1716 // Upsilon pT
1717 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
1718 //
1719   return PtUpsilonPbPb2760ShFdummy(*px, 5) * PtUpsilonPP2760(px, dummy);
1720 }
1721
1722 Double_t AliGenMUONlib::PtUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
1723 {
1724 // Upsilon pT
1725 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
1726 //
1727   return PtUpsilonPbPb2760ShFdummy(*px, 6) * PtUpsilonPP2760(px, dummy);
1728 }
1729
1730 Double_t AliGenMUONlib::PtUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
1731 {
1732 // Upsilon pT
1733 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
1734 //
1735   return PtUpsilonPbPb2760ShFdummy(*px, 7) * PtUpsilonPP2760(px, dummy);
1736 }
1737
1738 Double_t AliGenMUONlib::PtUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
1739 {
1740 // Upsilon pT
1741 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
1742 //
1743   return PtUpsilonPbPb2760ShFdummy(*px, 8) * PtUpsilonPP2760(px, dummy);
1744 }
1745
1746 Double_t AliGenMUONlib::PtUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
1747 {
1748 // Upsilon pT
1749 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
1750 //
1751   return PtUpsilonPbPb2760ShFdummy(*px, 9) * PtUpsilonPP2760(px, dummy);
1752 }
1753
1754 Double_t AliGenMUONlib::PtUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
1755 {
1756 // Upsilon pT
1757 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
1758 //
1759   return PtUpsilonPbPb2760ShFdummy(*px, 10) * PtUpsilonPP2760(px, dummy);
1760 }
1761
1762 Double_t AliGenMUONlib::PtUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
1763 {
1764 // Upsilon pT
1765 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
1766 //
1767   return PtUpsilonPbPb2760ShFdummy(*px, 11) * PtUpsilonPP2760(px, dummy);
1768 }
1769
1770 Double_t AliGenMUONlib::PtUpsilonPPb5030ShFdummy(Double_t x, Int_t n)
1771 {
1772 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1773 //
1774 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
1775 //
1776   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1777   const Double_t c[5] = {8.069e-01, 1.407e-04, 4.372e-04,-2.797e-05, 4.405e-06};
1778   Double_t y;
1779   Int_t j;
1780   y = c[j = 4];
1781   while (j > 0) y  = y * x + c[--j];
1782   y /= 1 + c[4]*TMath::Power(x,4);
1783   //  
1784   return 1 + (y-1)*f[n];
1785 }
1786
1787 Double_t AliGenMUONlib::PtUpsilonPPb5030(const Double_t *px, const Double_t *dummy)
1788 {
1789 // Upsilon pT
1790 // pPb 5.03 TeV, minimum bias 0-100 %
1791 //
1792   return PtUpsilonPPb5030ShFdummy(*px, 0) * PtUpsilonPP5030(px, dummy);
1793 }
1794
1795 Double_t AliGenMUONlib::PtUpsilonPPb5030c1(const Double_t *px, const Double_t *dummy)
1796 {
1797 // Upsilon pT
1798 // pPb 5.03 TeV, 1st centrality bin 0-20 %
1799 //
1800   return PtUpsilonPPb5030ShFdummy(*px, 1) * PtUpsilonPP5030(px, dummy);
1801 }
1802
1803 Double_t AliGenMUONlib::PtUpsilonPPb5030c2(const Double_t *px, const Double_t *dummy)
1804 {
1805 // Upsilon pT
1806 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
1807 //
1808   return PtUpsilonPPb5030ShFdummy(*px, 2) * PtUpsilonPP5030(px, dummy);
1809 }
1810
1811 Double_t AliGenMUONlib::PtUpsilonPPb5030c3(const Double_t *px, const Double_t *dummy)
1812 {
1813 // Upsilon pT
1814 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
1815 //
1816   return PtUpsilonPPb5030ShFdummy(*px, 3) * PtUpsilonPP5030(px, dummy);
1817 }
1818
1819 Double_t AliGenMUONlib::PtUpsilonPPb5030c4(const Double_t *px, const Double_t *dummy)
1820 {
1821 // Upsilon pT
1822 // pPb 5.03 TeV, 4th centrality bin 60-100 %
1823 //
1824   return PtUpsilonPPb5030ShFdummy(*px, 4) * PtUpsilonPP5030(px, dummy);
1825 }
1826
1827 Double_t AliGenMUONlib::PtUpsilonPbP5030ShFdummy(Double_t x, Int_t n)
1828 {
1829 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1830 //
1831 // Pbp 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
1832 //
1833   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1834   const Double_t c[5] = {1.122, 2.565e-03,-3.025e-04, 4.113e-06, 6.140e-07};
1835   Double_t y;
1836   Int_t j;
1837   y = c[j = 4];
1838   while (j > 0) y  = y * x + c[--j];
1839   y /= 1 + c[4]*TMath::Power(x,4);
1840   //  
1841   return 1 + (y-1)*f[n];
1842 }
1843
1844 Double_t AliGenMUONlib::PtUpsilonPbP5030(const Double_t *px, const Double_t *dummy)
1845 {
1846 // Upsilon pT
1847 // Pbp 5.03 TeV, minimum bias 0-100 %
1848 //
1849   return PtUpsilonPbP5030ShFdummy(*px, 0) * PtUpsilonPP5030(px, dummy);
1850 }
1851
1852 Double_t AliGenMUONlib::PtUpsilonPbP5030c1(const Double_t *px, const Double_t *dummy)
1853 {
1854 // Upsilon pT
1855 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
1856 //
1857   return PtUpsilonPbP5030ShFdummy(*px, 1) * PtUpsilonPP5030(px, dummy);
1858 }
1859
1860 Double_t AliGenMUONlib::PtUpsilonPbP5030c2(const Double_t *px, const Double_t *dummy)
1861 {
1862 // Upsilon pT
1863 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
1864 //
1865   return PtUpsilonPbP5030ShFdummy(*px, 2) * PtUpsilonPP5030(px, dummy);
1866 }
1867
1868 Double_t AliGenMUONlib::PtUpsilonPbP5030c3(const Double_t *px, const Double_t *dummy)
1869 {
1870 // Upsilon pT
1871 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
1872 //
1873   return PtUpsilonPbP5030ShFdummy(*px, 3) * PtUpsilonPP5030(px, dummy);
1874 }
1875
1876 Double_t AliGenMUONlib::PtUpsilonPbP5030c4(const Double_t *px, const Double_t *dummy)
1877 {
1878 // Upsilon pT
1879 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
1880 //
1881   return PtUpsilonPbP5030ShFdummy(*px, 4) * PtUpsilonPP5030(px, dummy);
1882 }
1883
1884 Double_t AliGenMUONlib::PtUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
1885 {
1886 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1887 //
1888 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1889 //
1890   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1891   const Double_t c[5] = {7.6561e-01, 1.1360e-04, 4.9596e-04,-3.0287e-05, 3.7555e-06};
1892   Double_t y;
1893   Int_t j;
1894   y = c[j = 4];
1895   while (j > 0) y  = y * x + c[--j];
1896   y /= 1 + c[4]*TMath::Power(x,4);
1897   //  
1898   return 1 + (y-1)*f[n];
1899 }
1900
1901 Double_t AliGenMUONlib::PtUpsilonPPb8800(const Double_t *px, const Double_t *dummy)
1902 {
1903 // Upsilon pT
1904 // pPb 8.8 TeV, minimum bias 0-100 %
1905 //
1906   return PtUpsilonPPb8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1907 }
1908
1909 Double_t AliGenMUONlib::PtUpsilonPPb8800c1(const Double_t *px, const Double_t *dummy)
1910 {
1911 // Upsilon pT
1912 // pPb 8.8 TeV, 1st centrality bin 0-20 %
1913 //
1914   return PtUpsilonPPb8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1915 }
1916
1917 Double_t AliGenMUONlib::PtUpsilonPPb8800c2(const Double_t *px, const Double_t *dummy)
1918 {
1919 // Upsilon pT
1920 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
1921 //
1922   return PtUpsilonPPb8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1923 }
1924
1925 Double_t AliGenMUONlib::PtUpsilonPPb8800c3(const Double_t *px, const Double_t *dummy)
1926 {
1927 // Upsilon pT
1928 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
1929 //
1930   return PtUpsilonPPb8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1931 }
1932
1933 Double_t AliGenMUONlib::PtUpsilonPPb8800c4(const Double_t *px, const Double_t *dummy)
1934 {
1935 // Upsilon pT
1936 // pPb 8.8 TeV, 4th centrality bin 60-100 %
1937 //
1938   return PtUpsilonPPb8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1939 }
1940
1941 Double_t AliGenMUONlib::PtUpsilonPbP8800ShFdummy(Double_t x, Int_t n)
1942 {
1943 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1944 //
1945 // Pbp 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1946 //
1947   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1948   const Double_t c[5] = {1.0975, 3.1905e-03,-2.0477e-04, 8.5270e-06, 2.5343e-06};
1949   Double_t y;
1950   Int_t j;
1951   y = c[j = 4];
1952   while (j > 0) y  = y * x + c[--j];
1953   y /= 1 + c[4]*TMath::Power(x,4);
1954   //  
1955   return 1 + (y-1)*f[n];
1956 }
1957
1958 Double_t AliGenMUONlib::PtUpsilonPbP8800(const Double_t *px, const Double_t *dummy)
1959 {
1960 // Upsilon pT
1961 // Pbp 8.8 TeV, minimum bias 0-100 %
1962 //
1963   return PtUpsilonPbP8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1964 }
1965
1966 Double_t AliGenMUONlib::PtUpsilonPbP8800c1(const Double_t *px, const Double_t *dummy)
1967 {
1968 // Upsilon pT
1969 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
1970 //
1971   return PtUpsilonPbP8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1972 }
1973
1974 Double_t AliGenMUONlib::PtUpsilonPbP8800c2(const Double_t *px, const Double_t *dummy)
1975 {
1976 // Upsilon pT
1977 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
1978 //
1979   return PtUpsilonPbP8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1980 }
1981
1982 Double_t AliGenMUONlib::PtUpsilonPbP8800c3(const Double_t *px, const Double_t *dummy)
1983 {
1984 // Upsilon pT
1985 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
1986 //
1987   return PtUpsilonPbP8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1988 }
1989
1990 Double_t AliGenMUONlib::PtUpsilonPbP8800c4(const Double_t *px, const Double_t *dummy)
1991 {
1992 // Upsilon pT
1993 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
1994 //
1995   return PtUpsilonPbP8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1996 }
1997
1998 Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
1999 {
2000 // Upsilon pT
2001   const Double_t kpt0 = 5.3;
2002   const Double_t kxn  = 2.5;
2003   Double_t x=*px;
2004   //
2005   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2006   return x/TMath::Power(pass1,kxn);
2007 }
2008
2009 Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
2010 {
2011 // Upsilon pT
2012   const Double_t kpt0 = 7.753;
2013   const Double_t kxn  = 3.042;
2014   Double_t x=*px;
2015   //
2016   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2017   return x/TMath::Power(pass1,kxn);
2018 }
2019
2020 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
2021 {
2022 // Upsilon pT
2023 //
2024 // pp 14 TeV
2025 //
2026 // scaled from CDF data at 2 TeV
2027
2028   const Double_t kpt0 = 8.610;
2029   const Double_t kxn  = 3.051;
2030   Double_t x=*px;
2031   //
2032   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2033   return x/TMath::Power(pass1,kxn);
2034 }
2035
2036 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2037 {
2038 // Upsilon pT
2039 //
2040 // pp 10 TeV
2041 //
2042 // scaled from CDF data at 2 TeV
2043
2044   const Double_t kpt0 = 8.235;
2045   const Double_t kxn  = 3.051;
2046   Double_t x=*px;
2047   //
2048   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2049   return x/TMath::Power(pass1,kxn);
2050 }
2051
2052 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2053 {
2054 // Upsilon pT
2055 //
2056 // pp 8.8 TeV
2057 // scaled from CDF data at 2 TeV
2058 //
2059   const Double_t kpt0 = 8.048;
2060   const Double_t kxn  = 3.051;
2061   Double_t x=*px;
2062   //
2063   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2064   return x/TMath::Power(pass1,kxn);
2065 }
2066
2067 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2068 {
2069 // Upsilon pT
2070 //
2071 // pp 7 TeV
2072 //
2073 // scaled from CDF data at 2 TeV
2074
2075   const Double_t kpt0 = 7.817;
2076   const Double_t kxn  = 3.051;
2077   Double_t x=*px;
2078   //
2079   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2080   return x/TMath::Power(pass1,kxn);
2081 }
2082
2083 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2084 {
2085 // Upsilon pT
2086 //
2087 // pp 3.94 TeV
2088 // scaled from CDF data at 2 TeV
2089 //
2090   const Double_t kpt0 = 7.189;
2091   const Double_t kxn  = 3.051;
2092   Double_t x=*px;
2093   //
2094   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2095   return x/TMath::Power(pass1,kxn);
2096 }
2097
2098 Double_t AliGenMUONlib::PtUpsilonCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
2099 {
2100 // Upsilon pT
2101 //
2102 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2103 //
2104   Double_t c[5] = {7.64952e-01, 1.12501e-04, 4.96038e-04, -3.03198e-05, 3.74035e-06};
2105   Double_t x=*px;
2106   Double_t y;
2107   Int_t j;
2108   y = c[j = 4];
2109   while (j > 0) y  = y * x + c[--j];
2110   //  
2111   Double_t d = 1.+c[4]*TMath::Power(x,4);
2112   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
2113 }
2114
2115 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
2116 {
2117 // Upsilon pT
2118 //
2119 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2120 //
2121   Double_t c[5] = {1.09881e+00, 3.08329e-03, -2.00356e-04, 8.28991e-06, 2.52576e-06};
2122   Double_t x=*px;
2123   Double_t y;
2124   Int_t j;
2125   y = c[j = 4];
2126   while (j > 0) y  = y * x + c[--j];
2127   //  
2128   Double_t d = 1.+c[4]*TMath::Power(x,4);
2129   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
2130 }
2131
2132 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2133 {
2134 // Upsilon pT
2135 //
2136 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2137 //
2138   Double_t c[5] = {8.65872e-01, 2.05465e-03, 2.56063e-04, -1.65598e-05, 2.29209e-06};
2139   Double_t x=*px;
2140   Double_t y;
2141   Int_t j;
2142   y = c[j = 4];
2143   while (j > 0) y  = y * x + c[--j];
2144   //  
2145   Double_t d = 1.+c[4]*TMath::Power(x,4);
2146   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP4(px,dummy);
2147 }
2148
2149 Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
2150 {
2151   return 1.;
2152 }
2153
2154 Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2155 {
2156 //
2157 // Upsilon pT
2158 //
2159 //
2160 // R. Vogt 2002
2161 // PbPb 5.5 TeV
2162 // MRST HO
2163 // mc = 1.4 GeV, pt-kick 1 GeV
2164 //
2165     Float_t x = px[0];
2166     Double_t c[8] = {
2167         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,       
2168         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
2169     };
2170     Double_t y;
2171     if (x < 10.) {
2172         Int_t j;
2173         y = c[j = 7];
2174         while (j > 0) y  = y * x +c[--j];
2175         y = x * TMath::Exp(y);
2176     } else {
2177         y = 0.;
2178     }
2179     return y;
2180 }
2181
2182 Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2183 {
2184 //
2185 // Upsilon pT
2186 //
2187 //
2188 // R. Vogt 2002
2189 // pp 14 TeV
2190 // MRST HO
2191 // mc = 1.4 GeV, pt-kick 1 GeV
2192 //
2193     Float_t x = px[0];
2194     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,   
2195                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
2196     
2197     Double_t y;
2198     if (x < 10.) {
2199         Int_t j;
2200         y = c[j = 7];
2201         while (j > 0) y  = y * x +c[--j];
2202         y = x * TMath::Exp(y);
2203     } else {
2204         y = 0.;
2205     }
2206     return y;
2207 }
2208
2209 //
2210 //                    y-distribution
2211 //
2212 //____________________________________________________________
2213 Double_t AliGenMUONlib::YUpsilonPPdummy(Double_t x, Double_t energy)
2214 {
2215 // Upsilon y
2216 // pp
2217 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
2218 //
2219     x = x/TMath::Log(energy/9.46);
2220     x = x*x;
2221     Double_t y = TMath::Exp(-x/0.4/0.4/2);
2222     if(x > 1) y=0;
2223     return y;
2224 }
2225
2226 Double_t AliGenMUONlib::YUpsilonPPpoly(Double_t x, Double_t energy)
2227 {
2228 // Upsilon y
2229 // pp
2230 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
2231 //
2232     x = x/TMath::Log(energy/9.46);
2233     x = x*x;
2234     Double_t y = 1 - 6.9*x*x;
2235     if(y < 0) y=0;
2236     return y;
2237 }
2238
2239 Double_t AliGenMUONlib::YUpsilonPP7000(const Double_t *px, const Double_t */*dummy*/)
2240 {
2241 // Upsilon y
2242 // pp 7 TeV
2243 //
2244   return YUpsilonPPdummy(*px, 7000);
2245 }
2246
2247 Double_t AliGenMUONlib::YUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
2248 {
2249 // Upsilon y
2250 // pp 2.76 TeV
2251 //
2252   return YUpsilonPPdummy(*px, 2760);
2253 }
2254
2255 Double_t AliGenMUONlib::YUpsilonPP4400(const Double_t *px, const Double_t */*dummy*/)
2256 {
2257 // Upsilon y
2258 // pp 4.4 TeV
2259 //
2260   return YUpsilonPPdummy(*px, 4400);
2261 }
2262
2263 Double_t AliGenMUONlib::YUpsilonPP5030(const Double_t *px, const Double_t */*dummy*/)
2264 {
2265 // Upsilon y
2266 // pp 5.03 TeV
2267 //
2268   return YUpsilonPPdummy(*px, 5030);
2269 }
2270
2271 Double_t AliGenMUONlib::YUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
2272 {
2273 // Upsilon y
2274 // pp 8.8 TeV
2275 //
2276   return YUpsilonPPdummy(*px, 8800);
2277 }
2278
2279 Double_t AliGenMUONlib::YUpsilonPPpoly7000(const Double_t *px, const Double_t */*dummy*/)
2280 {
2281 // Upsilon y
2282 // pp 7 TeV
2283 //
2284   return YUpsilonPPpoly(*px, 7000);
2285 }
2286
2287 Double_t AliGenMUONlib::YUpsilonPPpoly2760(const Double_t *px, const Double_t */*dummy*/)
2288 {
2289 // Upsilon y
2290 // pp 2.76 TeV
2291 //
2292   return YUpsilonPPpoly(*px, 2760);
2293 }
2294
2295 Double_t AliGenMUONlib::YUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
2296 {
2297 // Upsilon shadowing factor vs y for PbPb min. bias and 11 centr. bins
2298 //
2299 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
2300 //
2301   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
2302                            0.428, 0.317, 0.231, 0.156};
2303   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
2304                            0.106, 0.041, 0.013, 0.002};
2305   const Double_t c1[5] = {1.8547e+00, 1.6717e-02,-2.1285e-04,-9.7662e-05, 2.5768e-06};
2306   const Double_t c2[5] = {8.6029e-01, 1.1742e-02,-2.7944e-04,-6.7973e-05, 1.8838e-06}; 
2307
2308   x = x*x;
2309   Double_t y1, y2;
2310   Int_t j;
2311   y1 = c1[j = 4]; y2 = c2[4];
2312   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
2313   
2314   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
2315   if(y1<0) y1=0;
2316   return y1;
2317 }
2318
2319 Double_t AliGenMUONlib::YUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
2320 {
2321 // Upsilon y
2322 // PbPb 2.76 TeV, minimum bias 0-100 %
2323 //
2324   return YUpsilonPbPb2760ShFdummy(*px, 0) * YUpsilonPP2760(px, dummy);
2325 }
2326
2327 Double_t AliGenMUONlib::YUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
2328 {
2329 // Upsilon y
2330 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
2331 //
2332   return YUpsilonPbPb2760ShFdummy(*px, 1) * YUpsilonPP2760(px, dummy);
2333 }
2334
2335 Double_t AliGenMUONlib::YUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
2336 {
2337 // Upsilon y
2338 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
2339 //
2340   return YUpsilonPbPb2760ShFdummy(*px, 2) * YUpsilonPP2760(px, dummy);
2341 }
2342
2343 Double_t AliGenMUONlib::YUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
2344 {
2345 // Upsilon y
2346 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
2347 //
2348   return YUpsilonPbPb2760ShFdummy(*px, 3) * YUpsilonPP2760(px, dummy);
2349 }
2350
2351 Double_t AliGenMUONlib::YUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
2352 {
2353 // Upsilon y
2354 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
2355 //
2356   return YUpsilonPbPb2760ShFdummy(*px, 4) * YUpsilonPP2760(px, dummy);
2357 }
2358
2359 Double_t AliGenMUONlib::YUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
2360 {
2361 // Upsilon y
2362 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
2363 //
2364   return YUpsilonPbPb2760ShFdummy(*px, 5) * YUpsilonPP2760(px, dummy);
2365 }
2366
2367 Double_t AliGenMUONlib::YUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
2368 {
2369 // Upsilon y
2370 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
2371 //
2372   return YUpsilonPbPb2760ShFdummy(*px, 6) * YUpsilonPP2760(px, dummy);
2373 }
2374
2375 Double_t AliGenMUONlib::YUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
2376 {
2377 // Upsilon y
2378 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
2379 //
2380   return YUpsilonPbPb2760ShFdummy(*px, 7) * YUpsilonPP2760(px, dummy);
2381 }
2382
2383 Double_t AliGenMUONlib::YUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
2384 {
2385 // Upsilon y
2386 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
2387 //
2388   return YUpsilonPbPb2760ShFdummy(*px, 8) * YUpsilonPP2760(px, dummy);
2389 }
2390
2391 Double_t AliGenMUONlib::YUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
2392 {
2393 // Upsilon y
2394 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
2395 //
2396   return YUpsilonPbPb2760ShFdummy(*px, 9) * YUpsilonPP2760(px, dummy);
2397 }
2398
2399 Double_t AliGenMUONlib::YUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
2400 {
2401 // Upsilon y
2402 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
2403 //
2404   return YUpsilonPbPb2760ShFdummy(*px, 10) * YUpsilonPP2760(px, dummy);
2405 }
2406
2407 Double_t AliGenMUONlib::YUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
2408 {
2409 // Upsilon y
2410 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
2411 //
2412   return YUpsilonPbPb2760ShFdummy(*px, 11) * YUpsilonPP2760(px, dummy);
2413 }
2414
2415 Double_t AliGenMUONlib::YUpsilonPP5030dummy(Double_t px)
2416 {
2417     return AliGenMUONlib::YUpsilonPP5030(&px, (Double_t*) 0);
2418 }
2419
2420 Double_t AliGenMUONlib::YUpsilonPPb5030ShFdummy(Double_t x, Int_t n)
2421 {
2422 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2423 //
2424 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
2425 //
2426     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2427     const Double_t c[7] = {8.885e-01, 4.620e-02, 1.158e-02, 4.959e-04,-4.422e-04,-5.345e-05, 0.};
2428     Double_t y;
2429     Int_t j;
2430     y = c[j = 6];
2431     while (j > 0) y = y * x + c[--j];
2432     if(y<0) y=0;
2433     //
2434     return 1 +(y-1)*f[n];
2435 }
2436
2437 Double_t AliGenMUONlib::YUpsilonPPb5030(const Double_t *px, const Double_t */*dummy*/)
2438 {
2439 // Upsilon y
2440 // pPb 5.03 TeV, minimum bias 0-100 %
2441 //
2442   Double_t x = px[0] + 0.47;              // rapidity shift
2443   return YUpsilonPPb5030ShFdummy(x, 0) * YUpsilonPP5030dummy(x);
2444 }
2445
2446 Double_t AliGenMUONlib::YUpsilonPPb5030c1(const Double_t *px, const Double_t */*dummy*/)
2447 {
2448 // Upsilon y
2449 // pPb 5.03 TeV, 1st centrality bin 0-20 %
2450 //
2451   Double_t x = px[0] + 0.47;              // rapidity shift
2452   return YUpsilonPPb5030ShFdummy(x, 1) * YUpsilonPP5030dummy(x);
2453 }
2454
2455 Double_t AliGenMUONlib::YUpsilonPPb5030c2(const Double_t *px, const Double_t */*dummy*/)
2456 {
2457 // Upsilon y
2458 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
2459 //
2460   Double_t x = px[0] + 0.47;              // rapidity shift
2461   return YUpsilonPPb5030ShFdummy(x, 2) * YUpsilonPP5030dummy(x);
2462 }
2463
2464 Double_t AliGenMUONlib::YUpsilonPPb5030c3(const Double_t *px, const Double_t */*dummy*/)
2465 {
2466 // Upsilon y
2467 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
2468 //
2469   Double_t x = px[0] + 0.47;              // rapidity shift
2470   return YUpsilonPPb5030ShFdummy(x, 3) * YUpsilonPP5030dummy(x);
2471 }
2472
2473 Double_t AliGenMUONlib::YUpsilonPPb5030c4(const Double_t *px, const Double_t */*dummy*/)
2474 {
2475 // Upsilon y
2476 // pPb 5.03 TeV, 4th centrality bin 60-100 %
2477 //
2478   Double_t x = px[0] + 0.47;              // rapidity shift
2479   return YUpsilonPPb5030ShFdummy(x, 4) * YUpsilonPP5030dummy(x);
2480 }
2481
2482 Double_t AliGenMUONlib::YUpsilonPbP5030(const Double_t *px, const Double_t */*dummy*/)
2483 {
2484 // Upsilon y
2485 // Pbp 5.03 TeV, minimum bias 0-100 %
2486 //
2487   Double_t x = -px[0] + 0.47;              // rapidity shift
2488   return YUpsilonPPb5030ShFdummy(x, 0) * YUpsilonPP5030dummy(x);
2489 }
2490
2491 Double_t AliGenMUONlib::YUpsilonPbP5030c1(const Double_t *px, const Double_t */*dummy*/)
2492 {
2493 // Upsilon y
2494 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
2495 //
2496   Double_t x = -px[0] + 0.47;              // rapidity shift
2497   return YUpsilonPPb5030ShFdummy(x, 1) * YUpsilonPP5030dummy(x);
2498 }
2499
2500 Double_t AliGenMUONlib::YUpsilonPbP5030c2(const Double_t *px, const Double_t */*dummy*/)
2501 {
2502 // Upsilon y
2503 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
2504 //
2505   Double_t x = -px[0] + 0.47;              // rapidity shift
2506   return YUpsilonPPb5030ShFdummy(x, 2) * YUpsilonPP5030dummy(x);
2507 }
2508
2509 Double_t AliGenMUONlib::YUpsilonPbP5030c3(const Double_t *px, const Double_t */*dummy*/)
2510 {
2511 // Upsilon y
2512 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
2513 //
2514   Double_t x = -px[0] + 0.47;              // rapidity shift
2515   return YUpsilonPPb5030ShFdummy(x, 3) * YUpsilonPP5030dummy(x);
2516 }
2517
2518 Double_t AliGenMUONlib::YUpsilonPbP5030c4(const Double_t *px, const Double_t */*dummy*/)
2519 {
2520 // Upsilon y
2521 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
2522 //
2523   Double_t x = -px[0] + 0.47;              // rapidity shift
2524   return YUpsilonPPb5030ShFdummy(x, 4) * YUpsilonPP5030dummy(x);
2525 }
2526
2527 Double_t AliGenMUONlib::YUpsilonPP8800dummy(Double_t px)
2528 {
2529     return AliGenMUONlib::YUpsilonPP8800(&px, (Double_t*) 0);
2530 }
2531
2532 Double_t AliGenMUONlib::YUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
2533 {
2534 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2535 //
2536 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
2537 //
2538     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2539     const Double_t c[7] = {8.6581e-01, 4.6111e-02, 7.6911e-03, 8.7313e-04,
2540                            -1.4700e-04,-5.0975e-05,-3.5718e-06}; 
2541     Double_t y;
2542     Int_t j;
2543     y = c[j = 6];
2544     while (j > 0) y = y * x + c[--j];
2545     if(y<0) y=0;
2546     //
2547     return 1 +(y-1)*f[n];
2548 }
2549
2550 Double_t AliGenMUONlib::YUpsilonPPb8800(const Double_t *px, const Double_t */*dummy*/)
2551 {
2552 // Upsilon y
2553 // pPb 8.8 TeV, minimum bias 0-100 %
2554 //
2555   Double_t x = px[0] + 0.47;              // rapidity shift
2556   return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2557 }
2558
2559 Double_t AliGenMUONlib::YUpsilonPPb8800c1(const Double_t *px, const Double_t */*dummy*/)
2560 {
2561 // Upsilon y
2562 // pPb 8.8 TeV, 1st centrality bin 0-20 %
2563 //
2564   Double_t x = px[0] + 0.47;              // rapidity shift
2565   return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2566 }
2567
2568 Double_t AliGenMUONlib::YUpsilonPPb8800c2(const Double_t *px, const Double_t */*dummy*/)
2569 {
2570 // Upsilon y
2571 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
2572 //
2573   Double_t x = px[0] + 0.47;              // rapidity shift
2574   return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2575 }
2576
2577 Double_t AliGenMUONlib::YUpsilonPPb8800c3(const Double_t *px, const Double_t */*dummy*/)
2578 {
2579 // Upsilon y
2580 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
2581 //
2582   Double_t x = px[0] + 0.47;              // rapidity shift
2583   return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2584 }
2585
2586 Double_t AliGenMUONlib::YUpsilonPPb8800c4(const Double_t *px, const Double_t */*dummy*/)
2587 {
2588 // Upsilon y
2589 // pPb 8.8 TeV, 4th centrality bin 60-100 %
2590 //
2591   Double_t x = px[0] + 0.47;              // rapidity shift
2592   return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2593 }
2594
2595 Double_t AliGenMUONlib::YUpsilonPbP8800(const Double_t *px, const Double_t */*dummy*/)
2596 {
2597 // Upsilon y
2598 // Pbp 8.8 TeV, minimum bias 0-100 %
2599 //
2600   Double_t x = -px[0] + 0.47;              // rapidity shift
2601   return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2602 }
2603
2604 Double_t AliGenMUONlib::YUpsilonPbP8800c1(const Double_t *px, const Double_t */*dummy*/)
2605 {
2606 // Upsilon y
2607 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
2608 //
2609   Double_t x = -px[0] + 0.47;              // rapidity shift
2610   return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2611 }
2612
2613 Double_t AliGenMUONlib::YUpsilonPbP8800c2(const Double_t *px, const Double_t */*dummy*/)
2614 {
2615 // Upsilon y
2616 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
2617 //
2618   Double_t x = -px[0] + 0.47;              // rapidity shift
2619   return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2620 }
2621
2622 Double_t AliGenMUONlib::YUpsilonPbP8800c3(const Double_t *px, const Double_t */*dummy*/)
2623 {
2624 // Upsilon y
2625 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
2626 //
2627   Double_t x = -px[0] + 0.47;              // rapidity shift
2628   return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2629 }
2630
2631 Double_t AliGenMUONlib::YUpsilonPbP8800c4(const Double_t *px, const Double_t */*dummy*/)
2632 {
2633 // Upsilon y
2634 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
2635 //
2636   Double_t x = -px[0] + 0.47;              // rapidity shift
2637   return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2638 }
2639
2640 Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
2641 {
2642 // Upsilon y
2643   const Double_t ky0 = 3.;
2644   const Double_t kb=1.;
2645   Double_t yu;
2646   Double_t y=TMath::Abs(*py);
2647   //
2648   if (y < ky0)
2649     yu=kb;
2650   else
2651     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2652   return yu;
2653 }
2654
2655
2656 Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2657 {
2658
2659 //
2660 // Upsilon y
2661 //
2662 //
2663 // R. Vogt 2002
2664 // PbPb 5.5 TeV
2665 // MRST HO
2666 // mc = 1.4 GeV, pt-kick 1 GeV
2667 //
2668
2669     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
2670                      -2.99753e-09, 1.28895e-05};
2671     Double_t x = TMath::Abs(px[0]);
2672     if (x > 5.55) return 0.;
2673     Int_t j;
2674     Double_t y = c[j = 6];
2675     while (j > 0) y  = y * x +c[--j];
2676     return y;
2677 }
2678
2679 Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
2680 {
2681     // Upsilon y
2682     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
2683     
2684 }
2685
2686 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
2687 {
2688     // Upsilon y
2689     return AliGenMUONlib::YUpsilonPP(px, dummy);
2690     
2691 }
2692
2693 Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
2694 {
2695     // Upsilon y
2696     return 1.;
2697     
2698 }
2699
2700 Double_t AliGenMUONlib::YUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2701 {
2702 // Upsilon y
2703 //
2704 // pp 10 TeV
2705 // scaled from YUpsilonPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
2706 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
2707 //
2708     Double_t c[4] = {1., -2.17877e-02, -6.52830e-04, 1.40578e-05};
2709     Double_t x = TMath::Abs(px[0]);
2710     if (x > 6.1) return 0.;
2711     Int_t j;
2712     Double_t y = c[j = 3];
2713     while (j > 0) y  = y * x*x +c[--j];
2714     return y;
2715 }
2716
2717 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2718 {
2719 // Upsilon y
2720 //
2721 // pp 8.8 TeV
2722 // rescaling of YUpsilonPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD 
2723 //
2724     Double_t c[4] = {1., -2.37621e-02, -6.29610e-04, 1.47976e-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::YUpsilonCDFscaledPP9dummy(Double_t px)
2734 {
2735     return AliGenMUONlib::YUpsilonCDFscaledPP9(&px, (Double_t*) 0);
2736 }
2737
2738 Double_t AliGenMUONlib::YUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2739 {
2740 // Upsilon y
2741 //
2742 // pp 7 TeV
2743 // scaled from YUpsilonPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
2744 //
2745     Double_t c[4] = {1., -2.61009e-02, -6.83937e-04, 1.78451e-05};
2746     Double_t x = TMath::Abs(px[0]);
2747     if (x > 6.0) return 0.;
2748     Int_t j;
2749     Double_t y = c[j = 3];
2750     while (j > 0) y  = y * x*x +c[--j];
2751     return y;
2752 }
2753
2754 Double_t AliGenMUONlib::YUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2755 {
2756 // Upsilon y
2757 //
2758 // pp 3.94 TeV
2759 // rescaling of YUpsilonPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
2760 //
2761     Double_t c[4] = {1., -3.91924e-02, -4.26184e-04, 2.10914e-05};
2762     Double_t x = TMath::Abs(px[0]);
2763     if (x > 5.7) return 0.;
2764     Int_t j;
2765     Double_t y = c[j = 3];
2766     while (j > 0) y  = y * x*x +c[--j];
2767
2768     return y;
2769 }
2770
2771 Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2772 {
2773
2774 //
2775 // Upsilon y
2776 //
2777 //
2778 // R. Vogt 2002
2779 // p p  14. TeV
2780 // MRST HO
2781 // mc = 1.4 GeV, pt-kick 1 GeV
2782 //
2783     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, 
2784                      -6.20539e-10, 1.29943e-05};
2785     Double_t x = TMath::Abs(px[0]);
2786     if (x > 6.2) return 0.;
2787     Int_t j;
2788     Double_t y = c[j = 6];
2789     while (j > 0) y  = y * x +c[--j];
2790     return y;
2791 }
2792
2793 Double_t AliGenMUONlib::YUpsilonCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
2794 {
2795 // Upsilon y
2796 //
2797 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2798 //
2799     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2800                      -4.67538e-05,-2.11683e-06}; 
2801     Double_t y;
2802     Double_t x = px[0] + 0.47;              // rapidity shift
2803     Int_t j;
2804     y = c[j = 6];
2805     while (j > 0) y = y * x + c[--j];
2806     if(y<0) y=0;
2807
2808     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2809 }
2810
2811 Double_t AliGenMUONlib::YUpsilonCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
2812 {
2813 // Upsilon y
2814 //
2815 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2816 //
2817     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2818                      -4.67538e-05,-2.11683e-06}; 
2819     Double_t y;
2820     Double_t x = -px[0] + 0.47;              // rapidity shift
2821     Int_t j;
2822     y = c[j = 6];
2823     while (j > 0) y = y * x + c[--j];
2824     if(y<0) y=0;
2825
2826     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2827 }
2828
2829 Double_t AliGenMUONlib::YUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2830 {
2831 // Upsilon y
2832 //
2833 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2834 //
2835     Double_t c[4] = {8.27837e-01, 1.70115e-02, -1.26046e-03, 1.52091e-05}; 
2836     Double_t x = px[0]*px[0];
2837     Double_t y;
2838     Int_t j;
2839     y = c[j = 3];
2840     while (j > 0) y  = y * x + c[--j];
2841     if(y<0) y=0;
2842
2843     return y * AliGenMUONlib::YUpsilonCDFscaledPP4(px,dummy);
2844 }
2845
2846
2847 //                 particle composition
2848 //
2849 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
2850 {
2851 // y composition
2852     return 553;
2853 }
2854 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
2855 {
2856 // y composition
2857     return 100553;
2858 }
2859 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
2860 {
2861 // y composition
2862     return 200553;
2863 }
2864 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
2865 {
2866 // y composition
2867   // Using the LHCb pp data at 7 TeV: CERN-PH-EP-2012-051
2868   // (L. Manceau, S. Grigoryan)
2869   Int_t ip;
2870   Float_t r = gRandom->Rndm();  
2871   if (r < 0.687) {
2872     //  if (r < 0.712) {
2873     ip = 553;
2874   } else if (r < 0.903) {
2875     //  } else if (r < 0.896) {
2876     ip = 100553;
2877   } else {
2878     ip = 200553;
2879   }
2880   return ip;
2881 }
2882
2883
2884 //
2885 //                        Phi
2886 //
2887 //
2888 //    pt-distribution (by scaling of pion distribution)
2889 //____________________________________________________________
2890 Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
2891 {
2892 // Phi pT
2893   return PtScal(*px,7);
2894 }
2895 //    y-distribution
2896 Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
2897 {
2898 // Phi y
2899     Double_t *dum=0;
2900     return YJpsi(px,dum);
2901 }
2902 //                 particle composition
2903 //
2904 Int_t AliGenMUONlib::IpPhi(TRandom *)
2905 {
2906 // Phi composition
2907     return 333;
2908 }
2909
2910 //
2911 //                        omega
2912 //
2913 //
2914 //    pt-distribution (by scaling of pion distribution)
2915 //____________________________________________________________
2916 Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
2917 {
2918 // Omega pT
2919   return PtScal(*px,5);
2920 }
2921 //    y-distribution
2922 Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
2923 {
2924 // Omega y
2925     Double_t *dum=0;
2926     return YJpsi(px,dum);
2927 }
2928 //                 particle composition
2929 //
2930 Int_t AliGenMUONlib::IpOmega(TRandom *)
2931 {
2932 // Omega composition
2933     return 223;
2934 }
2935
2936
2937 //
2938 //                        Eta
2939 //
2940 //
2941 //    pt-distribution (by scaling of pion distribution)
2942 //____________________________________________________________
2943 Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
2944 {
2945 // Eta pT
2946   return PtScal(*px,3);
2947 }
2948 //    y-distribution
2949 Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
2950 {
2951 // Eta y
2952     Double_t *dum=0;
2953     return YJpsi(px,dum);
2954 }
2955 //                 particle composition
2956 //
2957 Int_t AliGenMUONlib::IpEta(TRandom *)
2958 {
2959 // Eta composition
2960     return 221;
2961 }
2962
2963 //
2964 //                        Charm
2965 //
2966 //
2967 //                    pt-distribution
2968 //____________________________________________________________
2969 Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
2970 {
2971 // Charm pT
2972   const Double_t kpt0 = 2.25;
2973   const Double_t kxn  = 3.17;
2974   Double_t x=*px;
2975   //
2976   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2977   return x/TMath::Power(pass1,kxn);
2978 }
2979
2980 Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
2981 {
2982 // Charm pT
2983   const Double_t kpt0 = 2.12;
2984   const Double_t kxn  = 2.78;
2985   Double_t x=*px;
2986   //
2987   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2988   return x/TMath::Power(pass1,kxn);
2989 }
2990 Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2991 {
2992 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
2993 // PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
2994 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
2995 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
2996 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
2997 // calculations for the following inputs: 
2998 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11 
2999 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV 
3000 // for j=0,1 & 2 respectively; 
3001 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
3002 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 
3003 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set 
3004 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3005 // June 2008, Smbat.Grigoryan@cern.ch
3006
3007 // Charm pT
3008 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
3009 // for pp collisions at 14 TeV with one c-cbar pair per event.
3010 // Corresponding NLO total cross section is 5.68 mb
3011
3012
3013   const Double_t kpt0 = 2.2930;
3014   const Double_t kxn  = 3.1196;
3015   Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
3016   Double_t x=*px;
3017   //
3018   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3019   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3020 }
3021 Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3022 {
3023 // Charm pT
3024 // Corresponding NLO total cross section is 6.06 mb
3025   const Double_t kpt0 = 2.8669;
3026   const Double_t kxn  = 3.1044;
3027   Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
3028   Double_t x=*px;
3029   //
3030   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3031   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3032 }
3033 Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3034 {
3035 // Charm pT
3036 // Corresponding NLO total cross section is 6.06 mb
3037   const Double_t kpt0 = 1.8361;
3038   const Double_t kxn  = 3.2966;
3039   Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
3040   Double_t x=*px;
3041   //
3042   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3043   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3044 }
3045 Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3046 {
3047 // Charm pT
3048 // Corresponding NLO total cross section is 7.69 mb
3049   const Double_t kpt0 = 2.1280;
3050   const Double_t kxn  = 3.1397;
3051   Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
3052   Double_t x=*px;
3053   //
3054   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3055   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3056 }
3057 Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3058 {
3059 // Charm pT
3060 // Corresponding NLO total cross section is 4.81 mb
3061   const Double_t kpt0 = 2.4579;
3062   const Double_t kxn  = 3.1095;
3063   Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
3064   Double_t x=*px;
3065   //
3066   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3067   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3068 }
3069 Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3070 {
3071 // Charm pT
3072 // Corresponding NLO total cross section is 14.09 mb
3073   const Double_t kpt0 = 2.1272;
3074   const Double_t kxn  = 3.1904;
3075   Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
3076   Double_t x=*px;
3077   //
3078   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3079   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3080 }
3081 Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3082 {
3083 // Charm pT
3084 // Corresponding NLO total cross section is 1.52 mb
3085   const Double_t kpt0 = 2.8159;
3086   const Double_t kxn  = 3.0857;
3087   Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
3088   Double_t x=*px;
3089   //
3090   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3091   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3092 }
3093 Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3094 {
3095 // Charm pT
3096 // Corresponding NLO total cross section is 3.67 mb
3097   const Double_t kpt0 = 2.7297;
3098   const Double_t kxn  = 3.3019;
3099   Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
3100   Double_t x=*px;
3101   //
3102   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3103   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3104 }
3105 Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3106 {
3107 // Charm pT
3108 // Corresponding NLO total cross section is 3.38 mb
3109   const Double_t kpt0 = 2.3894;
3110   const Double_t kxn  = 3.1075;
3111   Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
3112   Double_t x=*px;
3113   //
3114   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3115   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3116 }
3117 Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3118 {
3119 // Charm pT
3120 // Corresponding NLO total cross section is 10.37 mb
3121   const Double_t kpt0 = 2.0187;
3122   const Double_t kxn  = 3.3011;
3123   Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
3124   Double_t x=*px;
3125   //
3126   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3127   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3128 }
3129 Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3130 {
3131 // Charm pT
3132 // Corresponding NLO total cross section is 7.22 mb
3133   const Double_t kpt0 = 2.1089;
3134   const Double_t kxn  = 3.1848;
3135   Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
3136   Double_t x=*px;
3137   //
3138   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3139   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3140 }
3141
3142 //                  y-distribution
3143 Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
3144 {
3145 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
3146 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3147 // shadowing + kt broadening 
3148
3149     Double_t x=px[0];
3150     Double_t c[2]={-2.42985e-03,-2.31001e-04};
3151     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
3152     Double_t ycharm;
3153     
3154     if (TMath::Abs(x)>8) {
3155       ycharm=0.;
3156     }
3157     else {
3158       ycharm=TMath::Power(y,3);
3159     }
3160     
3161     return ycharm;
3162 }
3163 Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3164 {
3165 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3166 // YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3167 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3168 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
3169 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
3170 // calculations for the following inputs: 
3171 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11 
3172 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV 
3173 // for j=0,1 & 2 respectively; 
3174 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
3175 // with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for 
3176 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3177 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3178 // June 2008, Smbat.Grigoryan@cern.ch
3179
3180 // Charm y
3181 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
3182 // for pp collisions at 14 TeV with one c-cbar pair per event.
3183 // Corresponding NLO total cross section is 5.68 mb
3184
3185     Double_t x=px[0];
3186     Double_t c[2]={7.0909e-03,6.1967e-05};
3187     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3188     Double_t ycharm;
3189     
3190     if (TMath::Abs(x)>9) {
3191       ycharm=0.;
3192     }
3193     else {
3194       ycharm=TMath::Power(y,3);
3195     }
3196     
3197     return ycharm;
3198 }
3199 Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3200 {
3201 // Charm y
3202 // Corresponding NLO total cross section is 6.06 mb
3203     Double_t x=px[0];
3204     Double_t c[2]={6.9707e-03,6.0971e-05};
3205     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3206     Double_t ycharm;
3207     
3208     if (TMath::Abs(x)>9) {
3209       ycharm=0.;
3210     }
3211     else {
3212       ycharm=TMath::Power(y,3);
3213     }
3214     
3215     return ycharm;
3216 }
3217 Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3218 {
3219 // Charm y
3220 // Corresponding NLO total cross section is 6.06 mb
3221     Double_t x=px[0];
3222     Double_t c[2]={7.1687e-03,6.5303e-05};
3223     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3224     Double_t ycharm;
3225     
3226     if (TMath::Abs(x)>9) {
3227       ycharm=0.;
3228     }
3229     else {
3230       ycharm=TMath::Power(y,3);
3231     }
3232     
3233     return ycharm;
3234 }
3235 Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3236 {
3237 // Charm y
3238 // Corresponding NLO total cross section is 7.69 mb
3239     Double_t x=px[0];
3240     Double_t c[2]={5.9090e-03,7.1854e-05};
3241     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3242     Double_t ycharm;
3243     
3244     if (TMath::Abs(x)>9) {
3245       ycharm=0.;
3246     }
3247     else {
3248       ycharm=TMath::Power(y,3);
3249     }
3250     
3251     return ycharm;
3252 }
3253 Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3254 {
3255 // Charm y
3256 // Corresponding NLO total cross section is 4.81 mb
3257     Double_t x=px[0];
3258     Double_t c[2]={8.0882e-03,5.5872e-05};
3259     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3260     Double_t ycharm;
3261     
3262     if (TMath::Abs(x)>9) {
3263       ycharm=0.;
3264     }
3265     else {
3266       ycharm=TMath::Power(y,3);
3267     }
3268     
3269     return ycharm;
3270 }
3271 Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3272 {
3273 // Charm y
3274 // Corresponding NLO total cross section is 14.09 mb
3275     Double_t x=px[0];
3276     Double_t c[2]={7.2520e-03,6.2691e-05};
3277     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3278     Double_t ycharm;
3279     
3280     if (TMath::Abs(x)>9) {
3281       ycharm=0.;
3282     }
3283     else {
3284       ycharm=TMath::Power(y,3);
3285     }
3286     
3287     return ycharm;
3288 }
3289 Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3290 {
3291 // Charm y
3292 // Corresponding NLO total cross section is 1.52 mb
3293     Double_t x=px[0];
3294     Double_t c[2]={1.1040e-04,1.4498e-04};
3295     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3296     Double_t ycharm;
3297     
3298     if (TMath::Abs(x)>9) {
3299       ycharm=0.;
3300     }
3301     else {
3302       ycharm=TMath::Power(y,3);
3303     }
3304     
3305     return ycharm;
3306 }
3307 Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3308 {
3309 // Charm y
3310 // Corresponding NLO total cross section is 3.67 mb
3311     Double_t x=px[0];
3312     Double_t c[2]={-3.1328e-03,1.8270e-04};
3313     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3314     Double_t ycharm;
3315     
3316     if (TMath::Abs(x)>9) {
3317       ycharm=0.;
3318     }
3319     else {
3320       ycharm=TMath::Power(y,3);
3321     }
3322     
3323     return ycharm;
3324 }
3325 Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3326 {
3327 // Charm y
3328 // Corresponding NLO total cross section is 3.38 mb
3329     Double_t x=px[0];
3330     Double_t c[2]={7.0865e-03,6.2532e-05};
3331     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3332     Double_t ycharm;
3333     
3334     if (TMath::Abs(x)>9) {
3335       ycharm=0.;
3336     }
3337     else {
3338       ycharm=TMath::Power(y,3);
3339     }
3340     
3341     return ycharm;
3342 }
3343 Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3344 {
3345 // Charm y
3346 // Corresponding NLO total cross section is 10.37 mb
3347     Double_t x=px[0];
3348     Double_t c[2]={7.7070e-03,5.3533e-05};
3349     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3350     Double_t ycharm;
3351     
3352     if (TMath::Abs(x)>9) {
3353       ycharm=0.;
3354     }
3355     else {
3356       ycharm=TMath::Power(y,3);
3357     }
3358     
3359     return ycharm;
3360 }
3361 Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3362 {
3363 // Charm y
3364 // Corresponding NLO total cross section is 7.22 mb
3365     Double_t x=px[0];
3366     Double_t c[2]={7.9195e-03,5.3823e-05};
3367     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3368     Double_t ycharm;
3369     
3370     if (TMath::Abs(x)>9) {
3371       ycharm=0.;
3372     }
3373     else {
3374       ycharm=TMath::Power(y,3);
3375     }
3376     
3377     return ycharm;
3378 }
3379
3380
3381 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
3382 {  
3383 // Charm composition
3384     Float_t random;
3385     Int_t ip;
3386 //    411,421,431,4122
3387     random = ran->Rndm();
3388 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
3389 //  >>>>> cf. tab 4 p 11
3390   
3391     if (random < 0.30) {                       
3392         ip=421;
3393     } else if (random < 0.60) {
3394         ip=-421;
3395     } else if (random < 0.70) {
3396         ip=411;
3397     } else if (random < 0.80) {
3398         ip=-411;
3399     } else if (random < 0.86) {
3400         ip=431;
3401     } else if (random < 0.92) {
3402         ip=-431;        
3403     } else if (random < 0.96) {
3404         ip=4122;
3405     } else {
3406         ip=-4122;
3407     }
3408     
3409     return ip;
3410 }
3411
3412 //
3413 //                        Beauty
3414 //
3415 //
3416 //                    pt-distribution
3417 //____________________________________________________________
3418 Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
3419 {
3420 // Beauty pT
3421   const Double_t kpt0 = 6.53;
3422   const Double_t kxn  = 3.59;
3423   Double_t x=*px;
3424   //
3425   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3426   return x/TMath::Power(pass1,kxn);
3427 }
3428
3429 Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
3430 {
3431 // Beauty pT
3432   const Double_t kpt0 = 6.14;
3433   const Double_t kxn  = 2.93;
3434   Double_t x=*px;
3435   //
3436   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3437   return x/TMath::Power(pass1,kxn);
3438 }
3439 Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3440 {
3441 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3442 // PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
3443 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3444 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
3445 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
3446 // calculations for the following inputs: 
3447 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004 
3448 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV 
3449 // for j=0,1 & 2 respectively; 
3450 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
3451 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for 
3452 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3453 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3454 // June 2008, Smbat.Grigoryan@cern.ch
3455
3456 // Beauty pT
3457 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3458 // for pp collisions at 14 TeV with one b-bbar pair per event.
3459 // Corresponding NLO total cross section is 0.494 mb
3460
3461   const Double_t kpt0 = 8.0575;
3462   const Double_t kxn  = 3.1921;
3463   Double_t x=*px;
3464   //
3465   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3466   return x/TMath::Power(pass1,kxn);
3467 }
3468 Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3469 {
3470 // Beauty pT
3471 // Corresponding NLO total cross section is 0.445 mb
3472   const Double_t kpt0 = 8.6239;
3473   const Double_t kxn  = 3.2911;
3474   Double_t x=*px;
3475   //
3476   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3477   return x/TMath::Power(pass1,kxn);
3478 }
3479 Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3480 {
3481 // Beauty pT
3482 // Corresponding NLO total cross section is 0.445 mb
3483   const Double_t kpt0 = 7.3367;
3484   const Double_t kxn  = 3.0692;
3485   Double_t x=*px;
3486   //
3487   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3488   return x/TMath::Power(pass1,kxn);
3489 }
3490 Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3491 {
3492 // Beauty pT
3493 // Corresponding NLO total cross section is 0.518 mb
3494   const Double_t kpt0 = 7.6409;
3495   const Double_t kxn  = 3.1364;
3496   Double_t x=*px;
3497   //
3498   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3499   return x/TMath::Power(pass1,kxn);
3500 }
3501 Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3502 {
3503 // Beauty pT
3504 // Corresponding NLO total cross section is 0.384 mb
3505   const Double_t kpt0 = 8.4948;
3506   const Double_t kxn  = 3.2546;
3507   Double_t x=*px;
3508   //
3509   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3510   return x/TMath::Power(pass1,kxn);
3511 }
3512 Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3513 {
3514 // Beauty pT
3515 // Corresponding NLO total cross section is 0.648 mb
3516   const Double_t kpt0 = 7.6631;
3517   const Double_t kxn  = 3.1621;
3518   Double_t x=*px;
3519   //
3520   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3521   return x/TMath::Power(pass1,kxn);
3522 }
3523 Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3524 {
3525 // Beauty pT
3526 // Corresponding NLO total cross section is 0.294 mb
3527   const Double_t kpt0 = 8.7245;
3528   const Double_t kxn  = 3.2213;
3529   Double_t x=*px;
3530   //
3531   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3532   return x/TMath::Power(pass1,kxn);
3533 }
3534 Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3535 {
3536 // Beauty pT
3537 // Corresponding NLO total cross section is 0.475 mb
3538   const Double_t kpt0 = 8.5296;
3539   const Double_t kxn  = 3.2187;
3540   Double_t x=*px;
3541   //
3542   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3543   return x/TMath::Power(pass1,kxn);
3544 }
3545 Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3546 {
3547 // Beauty pT
3548 // Corresponding NLO total cross section is 0.324 mb
3549   const Double_t kpt0 = 7.9440;
3550   const Double_t kxn  = 3.1614;
3551   Double_t x=*px;
3552   //
3553   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3554   return x/TMath::Power(pass1,kxn);
3555 }
3556 Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3557 {
3558 // Beauty pT
3559 // Corresponding NLO total cross section is 0.536 mb
3560   const Double_t kpt0 = 8.2408;
3561   const Double_t kxn  = 3.3029;
3562   Double_t x=*px;
3563   //
3564   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3565   return x/TMath::Power(pass1,kxn);
3566 }
3567 Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3568 {
3569 // Beauty pT
3570 // Corresponding NLO total cross section is 0.420 mb
3571   const Double_t kpt0 = 7.8041;
3572   const Double_t kxn  = 3.2094;
3573   Double_t x=*px;
3574   //
3575   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3576   return x/TMath::Power(pass1,kxn);
3577 }
3578
3579 //                     y-distribution
3580 Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
3581 {
3582 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
3583 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3584 // shadowing + kt broadening 
3585
3586     Double_t x=px[0];
3587     Double_t c[2]={-1.27590e-02,-2.42731e-04};
3588     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
3589     Double_t ybeauty;
3590     
3591     if (TMath::Abs(x)>6) {
3592       ybeauty=0.;
3593     }
3594     else {
3595       ybeauty=TMath::Power(y,3);
3596     }
3597     
3598     return ybeauty;
3599 }
3600 Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3601 {
3602 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3603 // YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3604 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3605 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
3606 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
3607 // calculations for the following inputs: 
3608 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004 
3609 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV 
3610 // for j=0,1 & 2 respectively; 
3611 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
3612 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 
3613 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set 
3614 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3615 // June 2008, Smbat.Grigoryan@cern.ch
3616
3617 // Beauty y
3618 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3619 // for pp collisions at 14 TeV with one b-bbar pair per event.
3620 // Corresponding NLO total cross section is 0.494 mb
3621
3622
3623     Double_t x=px[0];
3624     Double_t c[2]={1.2350e-02,9.2667e-05};
3625     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3626     Double_t ybeauty;
3627     
3628     if (TMath::Abs(x)>7.6) {
3629       ybeauty=0.;
3630     }
3631     else {
3632       ybeauty=TMath::Power(y,3);
3633     }
3634     
3635     return ybeauty;
3636 }
3637 Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3638 {
3639 // Beauty y
3640 // Corresponding NLO total cross section is 0.445 mb
3641     Double_t x=px[0];
3642     Double_t c[2]={1.2292e-02,9.1847e-05};
3643     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3644     Double_t ybeauty;
3645     
3646     if (TMath::Abs(x)>7.6) {
3647       ybeauty=0.;
3648     }
3649     else {
3650       ybeauty=TMath::Power(y,3);
3651     }
3652     
3653     return ybeauty;
3654 }
3655 Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3656 {
3657 // Beauty y
3658 // Corresponding NLO total cross section is 0.445 mb
3659     Double_t x=px[0];
3660     Double_t c[2]={1.2436e-02,9.3709e-05};
3661     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3662     Double_t ybeauty;
3663     
3664     if (TMath::Abs(x)>7.6) {
3665       ybeauty=0.;
3666     }
3667     else {
3668       ybeauty=TMath::Power(y,3);
3669     }
3670     
3671     return ybeauty;
3672 }
3673 Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3674 {
3675 // Beauty y
3676 // Corresponding NLO total cross section is 0.518 mb
3677     Double_t x=px[0];
3678     Double_t c[2]={1.1714e-02,1.0068e-04};
3679     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3680     Double_t ybeauty;
3681     
3682     if (TMath::Abs(x)>7.6) {
3683       ybeauty=0.;
3684     }
3685     else {
3686       ybeauty=TMath::Power(y,3);
3687     }
3688     
3689     return ybeauty;
3690 }
3691 Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3692 {
3693 // Beauty y
3694 // Corresponding NLO total cross section is 0.384 mb
3695     Double_t x=px[0];
3696     Double_t c[2]={1.2944e-02,8.5500e-05};
3697     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3698     Double_t ybeauty;
3699     
3700     if (TMath::Abs(x)>7.6) {
3701       ybeauty=0.;
3702     }
3703     else {
3704       ybeauty=TMath::Power(y,3);
3705     }
3706     
3707     return ybeauty;
3708 }
3709 Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3710 {
3711 // Beauty y
3712 // Corresponding NLO total cross section is 0.648 mb
3713     Double_t x=px[0];
3714     Double_t c[2]={1.2455e-02,9.2713e-05};
3715     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3716     Double_t ybeauty;
3717     
3718     if (TMath::Abs(x)>7.6) {
3719       ybeauty=0.;
3720     }
3721     else {
3722       ybeauty=TMath::Power(y,3);
3723     }
3724     
3725     return ybeauty;
3726 }
3727 Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3728 {
3729 // Beauty y
3730 // Corresponding NLO total cross section is 0.294 mb
3731     Double_t x=px[0];
3732     Double_t c[2]={1.0897e-02,1.1878e-04};
3733     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3734     Double_t ybeauty;
3735     
3736     if (TMath::Abs(x)>7.6) {
3737       ybeauty=0.;
3738     }
3739     else {
3740       ybeauty=TMath::Power(y,3);
3741     }
3742     
3743     return ybeauty;
3744 }
3745 Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3746 {
3747 // Beauty y
3748 // Corresponding NLO total cross section is 0.475 mb
3749     Double_t x=px[0];
3750     Double_t c[2]={1.0912e-02,1.1858e-04};
3751     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3752     Double_t ybeauty;
3753     
3754     if (TMath::Abs(x)>7.6) {
3755       ybeauty=0.;
3756     }
3757     else {
3758       ybeauty=TMath::Power(y,3);
3759     }
3760     
3761     return ybeauty;
3762 }
3763 Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3764 {
3765 // Beauty y
3766 // Corresponding NLO total cross section is 0.324 mb
3767     Double_t x=px[0];
3768     Double_t c[2]={1.2378e-02,9.2490e-05};
3769     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3770     Double_t ybeauty;
3771     
3772     if (TMath::Abs(x)>7.6) {
3773       ybeauty=0.;
3774     }
3775     else {
3776       ybeauty=TMath::Power(y,3);
3777     }
3778     
3779     return ybeauty;
3780 }
3781 Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3782 {
3783 // Beauty y
3784 // Corresponding NLO total cross section is 0.536 mb
3785     Double_t x=px[0];
3786     Double_t c[2]={1.2886e-02,8.2912e-05};
3787     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3788     Double_t ybeauty;
3789     
3790     if (TMath::Abs(x)>7.6) {
3791       ybeauty=0.;
3792     }
3793     else {
3794       ybeauty=TMath::Power(y,3);
3795     }
3796     
3797     return ybeauty;
3798 }
3799 Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3800 {
3801 // Beauty y
3802 // Corresponding NLO total cross section is 0.420 mb
3803     Double_t x=px[0];
3804     Double_t c[2]={1.3106e-02,8.0115e-05};
3805     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3806     Double_t ybeauty;
3807     
3808     if (TMath::Abs(x)>7.6) {
3809       ybeauty=0.;
3810     }
3811     else {
3812       ybeauty=TMath::Power(y,3);
3813     }
3814     
3815     return ybeauty;
3816 }
3817
3818 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
3819 {  
3820 // Beauty Composition
3821     Float_t random;
3822     Int_t ip;
3823     random = ran->Rndm(); 
3824     
3825 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
3826 //  >>>>>&nbs