Changes needed in the case that the serial event number in a file does not start...
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONlib.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17$Log$
113049ac 18Revision 1.17 2002/11/07 09:06:10 morsch
19J/Psi and Upsilon pt and y distributions from R. Vogt 2002 added.
20
05932df6 21Revision 1.16 2002/10/14 14:55:35 hristov
22Merging the VirtualMC branch to the main development branch (HEAD)
23
b9d0a01d 24Revision 1.14.6.1 2002/06/10 14:57:41 hristov
25Merged with v3-08-02
26
27Revision 1.15 2002/04/17 10:11:51 morsch
28Coding Rule violations corrected.
29
53904666 30Revision 1.14 2002/02/22 17:26:43 morsch
31Eta and omega added.
32
89512a3b 33Revision 1.13 2001/03/27 11:01:04 morsch
34Charm pt-distribution corrected. More realistic y-distribution for pi and K.
35
2280e6af 36Revision 1.12 2001/03/09 13:01:41 morsch
37- enum constants for paramterisation type (particle family) moved to AliGen*lib.h
38- use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
39
34f60c01 40Revision 1.11 2000/11/30 07:12:50 alibrary
41Introducing new Rndm and QA classes
42
65fb704d 43Revision 1.10 2000/06/29 21:08:27 morsch
44All paramatrisation libraries derive from the pure virtual base class AliGenLib.
45This allows to pass a pointer to a library directly to AliGenParam and avoids the
46use of function pointers in Config.C.
47
b22ee262 48Revision 1.9 2000/06/14 15:20:56 morsch
49Include clean-up (IH)
50
5c3fd7ea 51Revision 1.8 2000/06/09 20:32:11 morsch
52All coding rule violations except RS3 corrected
53
f87cfe57 54Revision 1.7 2000/05/02 08:12:13 morsch
55Coding rule violations corrected.
56
d90f80fd 57Revision 1.6 1999/09/29 09:24:14 fca
58Introduction of the Copyright and cvs Log
59
4c039060 60*/
61
53904666 62// Library class for particle pt and y distributions used for
63// muon spectrometer simulations.
64// To be used with AliGenParam.
65// The following particle typed can be simulated:
66// pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
67//
68// andreas.morsch@cern.ch
69//
70
65fb704d 71#include "TMath.h"
72#include "TRandom.h"
73
fe4da5cc 74#include "AliGenMUONlib.h"
5c3fd7ea 75
fe4da5cc 76ClassImp(AliGenMUONlib)
77//
78// Pions
d90f80fd 79Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
fe4da5cc 80{
81//
82// PT-PARAMETERIZATION CDF, PRL 61(88) 1819
83// POWER LAW FOR PT > 500 MEV
84// MT SCALING BELOW (T=160 MEV)
85//
d90f80fd 86 const Double_t kp0 = 1.3;
87 const Double_t kxn = 8.28;
88 const Double_t kxlim=0.5;
89 const Double_t kt=0.160;
90 const Double_t kxmpi=0.139;
91 const Double_t kb=1.;
fe4da5cc 92 Double_t y, y1, xmpi2, ynorm, a;
93 Double_t x=*px;
94 //
d90f80fd 95 y1=TMath::Power(kp0/(kp0+kxlim),kxn);
96 xmpi2=kxmpi*kxmpi;
97 ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
fe4da5cc 98 a=ynorm/y1;
d90f80fd 99 if (x > kxlim)
100 y=a*TMath::Power(kp0/(kp0+x),kxn);
fe4da5cc 101 else
d90f80fd 102 y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
fe4da5cc 103 return y*x;
104}
753690b0 105//
106// y-distribution
107//
d90f80fd 108Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
753690b0 109{
d90f80fd 110// Pion y
2280e6af 111 Double_t y=TMath::Abs(*py);
112/*
d90f80fd 113 const Double_t ka = 7000.;
114 const Double_t kdy = 4.;
d90f80fd 115 Double_t ex = y*y/(2*kdy*kdy);
116 return ka*TMath::Exp(-ex);
2280e6af 117*/
118 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
119
753690b0 120}
121// particle composition
122//
65fb704d 123Int_t AliGenMUONlib::IpPion(TRandom *ran)
753690b0 124{
d90f80fd 125// Pion composition
65fb704d 126 if (ran->Rndm() < 0.5) {
753690b0 127 return 211;
128 } else {
129 return -211;
130 }
131}
fe4da5cc 132
133//____________________________________________________________
134//
135// Mt-scaling
136
137Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
138{
139 // SCALING EN MASSE PAR RAPPORT A PTPI
140 // MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
d90f80fd 141 const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
fe4da5cc 142 // VALUE MESON/PI AT 5 GEV
d90f80fd 143 const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
fe4da5cc 144 np--;
d90f80fd 145 Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
146 Double_t fmax2=f5/kfmax[np];
fe4da5cc 147 // PIONS
148 Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
149 Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
d90f80fd 150 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
fe4da5cc 151 return fmtscal*ptpion;
152}
153//
753690b0 154// kaon
155//
156// pt-distribution
157//____________________________________________________________
d90f80fd 158Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
753690b0 159{
d90f80fd 160// Kaon pT
753690b0 161 return PtScal(*px,2);
162}
163
164// y-distribution
fe4da5cc 165//____________________________________________________________
d90f80fd 166Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
fe4da5cc 167{
d90f80fd 168// Kaon y
2280e6af 169 Double_t y=TMath::Abs(*py);
170/*
d90f80fd 171 const Double_t ka = 1000.;
172 const Double_t kdy = 4.;
fe4da5cc 173 //
d90f80fd 174 Double_t ex = y*y/(2*kdy*kdy);
175 return ka*TMath::Exp(-ex);
2280e6af 176*/
177
178 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
753690b0 179}
180
181// particle composition
182//
65fb704d 183Int_t AliGenMUONlib::IpKaon(TRandom *ran)
753690b0 184{
d90f80fd 185// Kaon composition
65fb704d 186 if (ran->Rndm() < 0.5) {
753690b0 187 return 321;
188 } else {
189 return -321;
190 }
fe4da5cc 191}
753690b0 192
fe4da5cc 193// J/Psi
194//
195//
196// pt-distribution
197//____________________________________________________________
d90f80fd 198Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
fe4da5cc 199{
d90f80fd 200// J/Psi pT
201 const Double_t kpt0 = 4.;
202 const Double_t kxn = 3.6;
fe4da5cc 203 Double_t x=*px;
204 //
d90f80fd 205 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
206 return x/TMath::Power(pass1,kxn);
fe4da5cc 207}
05932df6 208
209Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
210{
211// J/Psi pT
212//
213// R. Vogt 2002
214// PbPb 5.5 TeV
215// MRST HO
216// mc = 1.4 GeV, pt-kick 1 GeV
217//
218
219 Float_t ptJpsi[100] = {
220 0.0000e-01, 4.5870e+01, 6.5200e+01, 7.1740e+01, 6.5090e+01,
221 5.5070e+01, 4.9420e+01, 3.9780e+01, 3.2390e+01, 2.8120e+01,
222 2.3870e+01, 1.9540e+01, 1.6510e+01, 1.4180e+01, 1.2050e+01,
223 1.0390e+01, 8.7970e+00, 7.8680e+00, 6.7710e+00, 5.9360e+00,
224 5.3460e+00, 4.5670e+00, 4.6500e+00, 3.9360e+00, 3.5070e+00,
225 3.2070e+00, 2.8310e+00, 2.6340e+00, 2.4900e+00, 2.2410e+00,
226 2.1090e+00, 1.9070e+00, 1.7360e+00, 1.6120e+00, 1.5450e+00,
227 1.4350e+00, 1.3890e+00, 1.2610e+00, 1.0880e+00, 1.0930e+00,
228 1.0680e+00, 9.2500e-01, 8.6790e-01, 8.1790e-01, 7.9770e-01,
229 7.4660e-01, 7.3110e-01, 6.5120e-01, 6.8140e-01, 5.7960e-01,
230 5.8210e-01, 5.4640e-01, 5.1700e-01, 5.0760e-01, 4.8280e-01,
231 4.5360e-01, 4.4910e-01, 4.2410e-01, 4.2100e-01, 3.9530e-01,
232 3.7220e-01, 3.4840e-01, 3.4550e-01, 3.3000e-01, 3.1670e-01,
233 3.1470e-01, 2.8920e-01, 2.7650e-01, 2.6860e-01, 2.5390e-01,
234 2.4190e-01, 2.5200e-01, 2.2960e-01, 2.2540e-01, 2.0950e-01,
235 2.0250e-01, 1.8720e-01, 1.8200e-01, 1.7860e-01, 1.8290e-01,
236 1.6970e-01, 1.7130e-01, 1.6310e-01, 1.5500e-01, 1.5100e-01,
237 1.5770e-01, 1.4240e-01, 1.4560e-01, 1.3330e-01, 1.4190e-01,
238 1.2010e-01, 1.2430e-01, 1.2430e-01, 1.1340e-01, 1.1840e-01,
239 1.1380e-01, 1.0330e-01, 1.0130e-01, 1.0390e-01, 9.5810e-02
240 };
241 Float_t x = px[0] * px[0];
242 if (x < 1.5 || x > 100) {
243 return 0.0;
244 } else {
245 Float_t y = Interpolate(x, ptJpsi, 0.5, 1., 100, 4);
246 return px[0] * y;
247 }
248}
fe4da5cc 249//
250// y-distribution
251//____________________________________________________________
d90f80fd 252Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
fe4da5cc 253{
d90f80fd 254// J/psi y
255 const Double_t ky0 = 4.;
256 const Double_t kb=1.;
fe4da5cc 257 Double_t yj;
258 Double_t y=TMath::Abs(*py);
259 //
d90f80fd 260 if (y < ky0)
261 yj=kb;
fe4da5cc 262 else
d90f80fd 263 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
fe4da5cc 264 return yj;
265}
05932df6 266
267
268Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
269{
270
271//
272// J/Psi y
273//
274//
275// R. Vogt 2002
276// PbPb 5.5 TeV
277// MRST HO
278// mc = 1.4 GeV, pt-kick 1 GeV
279//
280
281 Float_t yJpsi[62] =
282 {
283 0.3981E-03, 0.1169E-01, 0.6143E-01, 0.3554E+00, 0.1249E+01
284 , 0.1677E+01, 0.3634E+01, 0.5414E+01, 0.8242E+01, 0.1102E+02
285 , 0.1353E+02, 0.1964E+02, 0.2357E+02, 0.2662E+02, 0.3023E+02
286 , 0.3250E+02, 0.3137E+02, 0.3243E+02, 0.3120E+02, 0.3249E+02
287 , 0.3166E+02, 0.3104E+02, 0.3203E+02, 0.3149E+02, 0.3117E+02
288 , 0.3210E+02, 0.3170E+02, 0.3279E+02, 0.3079E+02, 0.3208E+02
289 , 0.3218E+02, 0.3218E+02, 0.3208E+02, 0.3079E+02, 0.3279E+02
290 , 0.3170E+02, 0.3210E+02, 0.3118E+02, 0.3149E+02, 0.3203E+02
291 , 0.3104E+02, 0.3167E+02, 0.3250E+02, 0.3120E+02, 0.3243E+02
292 , 0.3137E+02, 0.3250E+02, 0.3022E+02, 0.2662E+02, 0.2357E+02
293 , 0.1964E+02, 0.1353E+02, 0.1102E+02, 0.8242E+01, 0.5414E+01
294 , 0.3634E+01, 0.1677E+01, 0.1249E+01, 0.3554E+00, 0.6142E-01
295 , 0.1169E-01, 0.3981E-03};
296
297 return Interpolate(px[0], yJpsi, -7.625, 0.25, 62, 2);
298}
299
fe4da5cc 300// particle composition
301//
65fb704d 302Int_t AliGenMUONlib::IpJpsi(TRandom *)
fe4da5cc 303{
d90f80fd 304// J/Psi composition
fe4da5cc 305 return 443;
306}
307
308// Upsilon
309//
310//
311// pt-distribution
312//____________________________________________________________
d90f80fd 313Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
fe4da5cc 314{
d90f80fd 315// Upsilon pT
316 const Double_t kpt0 = 5.3;
317 const Double_t kxn = 2.5;
fe4da5cc 318 Double_t x=*px;
319 //
d90f80fd 320 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
321 return x/TMath::Power(pass1,kxn);
fe4da5cc 322}
05932df6 323
324Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
325{
326
327//
328// Upsilon pT
329//
330//
331// R. Vogt 2002
332// PbPb 5.5 TeV
333// MRST HO
334// mc = 1.4 GeV, pt-kick 1 GeV
335//
336
337 Float_t ptUps[100] =
338 {
339 0.0000e-01, -1.5290e-02, 2.6020e-01, 3.4220e-01, 3.3710e-01,
340 3.1880e-01, 3.3420e-01, 2.7740e-01, 2.3730e-01, 2.0640e-01,
341 1.7690e-01, 1.6190e-01, 1.4500e-01, 1.3310e-01, 1.1440e-01,
342 1.0800e-01, 1.0210e-01, 8.4690e-02, 8.0050e-02, 7.0710e-02,
343 6.4160e-02, 6.5200e-02, 6.6890e-02, 6.0600e-02, 5.4030e-02,
344 5.1140e-02, 4.6120e-02, 4.4800e-02, 4.2490e-02, 4.1440e-02,
345 4.0310e-02, 3.7110e-02, 3.5890e-02, 3.5420e-02, 3.0370e-02,
346 2.9970e-02, 3.0770e-02, 2.6380e-02, 2.7740e-02, 2.6690e-02,
347 2.4210e-02, 2.5200e-02, 2.3760e-02, 2.1370e-02, 2.2290e-02,
348 2.2700e-02, 2.0110e-02, 1.9320e-02, 1.8830e-02, 1.9910e-02,
349 1.9740e-02, 1.8460e-02, 1.8240e-02, 1.6740e-02, 1.6140e-02,
350 1.7340e-02, 1.5950e-02, 1.5430e-02, 1.4780e-02, 1.2750e-02,
351 1.4370e-02, 1.2810e-02, 1.2900e-02, 1.1070e-02, 1.1830e-02,
352 1.1150e-02, 1.1260e-02, 1.1610e-02, 1.0700e-02, 1.1600e-02,
353 1.0390e-02, 1.0280e-02, 1.0180e-02, 1.0030e-02, 9.6050e-03,
354 8.8050e-03, 8.9680e-03, 9.0120e-03, 8.4110e-03, 8.6660e-03,
355 8.3060e-03, 8.5850e-03, 8.2600e-03, 8.3800e-03, 8.4200e-03,
356 7.5690e-03, 7.2100e-03, 7.1230e-03, 7.3350e-03, 7.1980e-03,
357 6.7500e-03, 6.6190e-03, 6.3370e-03, 6.6270e-03, 6.8290e-03,
358 6.0880e-03, 6.6310e-03, 6.0490e-03, 5.8900e-03, 5.6100e-03
359 };
360 Float_t x = px[0] * px[0];
361 if (x < 1.5 || x > 100) {
362 return 0.0;
363 } else {
364 Float_t y = Interpolate(x, ptUps, 0.5, 1., 100, 4);
365 return px[0] * y;
366 }
367}
368
fe4da5cc 369//
370// y-distribution
371//
372//____________________________________________________________
d90f80fd 373Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
fe4da5cc 374{
d90f80fd 375// Upsilon y
376 const Double_t ky0 = 3.;
377 const Double_t kb=1.;
fe4da5cc 378 Double_t yu;
379 Double_t y=TMath::Abs(*py);
380 //
d90f80fd 381 if (y < ky0)
382 yu=kb;
fe4da5cc 383 else
d90f80fd 384 yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
fe4da5cc 385 return yu;
386}
05932df6 387
388
389Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
390{
391
392//
393// Upsilon y
394//
395//
396// R. Vogt 2002
397// PbPb 5.5 TeV
398// MRST HO
399// mc = 1.4 GeV, pt-kick 1 GeV
400//
401
402 Float_t yUps[52] =
403 {
404 0.000000, 0.000065, 0.000767, 0.004332, 0.014790,
405 0.029370, 0.052060, 0.077930, 0.122500, 0.135800,
406 0.184000, 0.207900, 0.228300, 0.258500, 0.269500,
407 0.288500, 0.316200, 0.304100, 0.315800, 0.323300,
408 0.322400, 0.322600, 0.345500, 0.338100, 0.331900,
409 0.343700, 0.343700, 0.331900, 0.338100, 0.345500,
410 0.322600, 0.322400, 0.323300, 0.315800, 0.304100,
411 0.316200, 0.288500, 0.269500, 0.258500, 0.228300,
412 0.207900, 0.184000, 0.135800, 0.122500, 0.077930,
413 0.052060, 0.029380, 0.014780, 0.004332, 0.000767,
414 0.6479E-04, 0.1013E-06
415 };
416
417
418
419 return Interpolate(px[0], yUps, -6.5, 0.25, 52, 2);
420}
421
fe4da5cc 422// particle composition
423//
65fb704d 424Int_t AliGenMUONlib::IpUpsilon(TRandom *)
fe4da5cc 425{
d90f80fd 426// y composition
fe4da5cc 427 return 553;
428}
429
430//
431// Phi
432//
433//
434// pt-distribution (by scaling of pion distribution)
435//____________________________________________________________
d90f80fd 436Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
fe4da5cc 437{
d90f80fd 438// Phi pT
fe4da5cc 439 return PtScal(*px,7);
440}
441// y-distribution
d90f80fd 442Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
fe4da5cc 443{
d90f80fd 444// Phi y
445 Double_t *dum=0;
446 return YJpsi(px,dum);
fe4da5cc 447}
448// particle composition
449//
65fb704d 450Int_t AliGenMUONlib::IpPhi(TRandom *)
fe4da5cc 451{
d90f80fd 452// Phi composition
89512a3b 453 return 333;
454}
455
456//
457// omega
458//
459//
460// pt-distribution (by scaling of pion distribution)
461//____________________________________________________________
462Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
463{
464// Omega pT
465 return PtScal(*px,5);
466}
467// y-distribution
468Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
469{
470// Omega y
471 Double_t *dum=0;
472 return YJpsi(px,dum);
473}
474// particle composition
475//
476Int_t AliGenMUONlib::IpOmega(TRandom *)
477{
478// Omega composition
479 return 223;
480}
481
482
483//
484// Eta
485//
486//
487// pt-distribution (by scaling of pion distribution)
488//____________________________________________________________
489Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
490{
491// Eta pT
492 return PtScal(*px,3);
493}
494// y-distribution
495Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
496{
497// Eta y
498 Double_t *dum=0;
499 return YJpsi(px,dum);
500}
501// particle composition
502//
503Int_t AliGenMUONlib::IpEta(TRandom *)
504{
505// Eta composition
506 return 221;
fe4da5cc 507}
508
509//
510// Charm
511//
512//
513// pt-distribution
514//____________________________________________________________
d90f80fd 515Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
fe4da5cc 516{
d90f80fd 517// Charm pT
518 const Double_t kpt0 = 4.08;
519 const Double_t kxn = 9.40;
2280e6af 520
fe4da5cc 521 Double_t x=*px;
522 //
2280e6af 523 Double_t pass1 = 1.+(x/kpt0);
d90f80fd 524 return x/TMath::Power(pass1,kxn);
fe4da5cc 525}
526// y-distribution
d90f80fd 527Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
fe4da5cc 528{
d90f80fd 529// Charm y
530 Double_t *dum=0;
531 return YJpsi(px,dum);
fe4da5cc 532}
533
65fb704d 534Int_t AliGenMUONlib::IpCharm(TRandom *ran)
fe4da5cc 535{
d90f80fd 536// Charm composition
65fb704d 537 Float_t random;
fe4da5cc 538 Int_t ip;
539// 411,421,431,4122
65fb704d 540 random = ran->Rndm();
541 if (random < 0.5) {
fe4da5cc 542 ip=411;
65fb704d 543 } else if (random < 0.75) {
fe4da5cc 544 ip=421;
65fb704d 545 } else if (random < 0.90) {
fe4da5cc 546 ip=431;
547 } else {
548 ip=4122;
549 }
65fb704d 550 if (ran->Rndm() < 0.5) {ip=-ip;}
fe4da5cc 551
552 return ip;
553}
554
555
556//
557// Beauty
558//
559//
560// pt-distribution
561//____________________________________________________________
d90f80fd 562Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
fe4da5cc 563{
d90f80fd 564// Beauty pT
565 const Double_t kpt0 = 4.;
566 const Double_t kxn = 3.6;
fe4da5cc 567 Double_t x=*px;
568 //
d90f80fd 569 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
570 return x/TMath::Power(pass1,kxn);
fe4da5cc 571}
572// y-distribution
d90f80fd 573Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
fe4da5cc 574{
d90f80fd 575// Beauty y
576 Double_t *dum=0;
577 return YJpsi(px,dum);
fe4da5cc 578}
579
65fb704d 580Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
fe4da5cc 581{
d90f80fd 582// Beauty Composition
65fb704d 583 Float_t random;
fe4da5cc 584 Int_t ip;
65fb704d 585 random = ran->Rndm();
586 if (random < 0.5) {
fe4da5cc 587 ip=511;
65fb704d 588 } else if (random < 0.75) {
fe4da5cc 589 ip=521;
65fb704d 590 } else if (random < 0.90) {
fe4da5cc 591 ip=531;
592 } else {
593 ip=5122;
594 }
65fb704d 595 if (ran->Rndm() < 0.5) {ip=-ip;}
fe4da5cc 596
597 return ip;
598}
599
600typedef Double_t (*GenFunc) (Double_t*, Double_t*);
53904666 601GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname) const
fe4da5cc 602{
d90f80fd 603// Return pointer to pT parameterisation
05932df6 604 TString sname = TString(tname);
fe4da5cc 605 GenFunc func;
753690b0 606 switch (param)
fe4da5cc 607 {
34f60c01 608 case kPhi:
fe4da5cc 609 func=PtPhi;
610 break;
89512a3b 611 case kOmega:
612 func=PtOmega;
613 break;
614 case kEta:
615 func=PtEta;
616 break;
34f60c01 617 case kJpsi:
113049ac 618 if (sname == "Vogt") {
05932df6 619 func=PtJpsiPbPb;
620 } else {
621 func=PtJpsi;
622 }
fe4da5cc 623 break;
34f60c01 624 case kUpsilon:
113049ac 625 if (sname == "Vogt") {
05932df6 626 func=PtUpsilonPbPb;
627 } else {
628 func=PtUpsilon;
629 }
fe4da5cc 630 break;
34f60c01 631 case kCharm:
fe4da5cc 632 func=PtCharm;
633 break;
34f60c01 634 case kBeauty:
fe4da5cc 635 func=PtBeauty;
636 break;
34f60c01 637 case kPion:
753690b0 638 func=PtPion;
639 break;
34f60c01 640 case kKaon:
753690b0 641 func=PtKaon;
642 break;
119b35c7 643 default:
644 func=0;
645 printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
fe4da5cc 646 }
647 return func;
648}
649
53904666 650GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
fe4da5cc 651{
05932df6 652 TString sname = TString(tname);
653
d90f80fd 654// Return pointer to y- parameterisation
fe4da5cc 655 GenFunc func;
753690b0 656 switch (param)
fe4da5cc 657 {
34f60c01 658 case kPhi:
fe4da5cc 659 func=YPhi;
660 break;
89512a3b 661 case kEta:
662 func=YEta;
663 break;
664 case kOmega:
665 func=YOmega;
666 break;
34f60c01 667 case kJpsi:
113049ac 668 if (sname == "Vogt") {
05932df6 669 func=YJpsiPbPb;
670 } else {
671 func=YJpsi;
672 }
fe4da5cc 673 break;
34f60c01 674 case kUpsilon:
113049ac 675 if (sname == "Vogt") {
05932df6 676 func=YUpsilonPbPb;
677 } else {
678 func=YUpsilon;
679 }
fe4da5cc 680 break;
34f60c01 681 case kCharm:
fe4da5cc 682 func=YCharm;
683 break;
34f60c01 684 case kBeauty:
fe4da5cc 685 func=YBeauty;
686 break;
34f60c01 687 case kPion:
753690b0 688 func=YPion;
689 break;
34f60c01 690 case kKaon:
753690b0 691 func=YKaon;
692 break;
119b35c7 693 default:
694 func=0;
695 printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
fe4da5cc 696 }
697 return func;
698}
65fb704d 699typedef Int_t (*GenFuncIp) (TRandom *);
53904666 700GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* tname) const
fe4da5cc 701{
d90f80fd 702// Return pointer to particle type parameterisation
fe4da5cc 703 GenFuncIp func;
753690b0 704 switch (param)
fe4da5cc 705 {
34f60c01 706 case kPhi:
fe4da5cc 707 func=IpPhi;
708 break;
89512a3b 709 case kEta:
710 func=IpEta;
711 break;
712 case kOmega:
713 func=IpOmega;
714 break;
34f60c01 715 case kJpsi:
fe4da5cc 716 func=IpJpsi;
717 break;
34f60c01 718 case kUpsilon:
fe4da5cc 719 func=IpUpsilon;
720 break;
34f60c01 721 case kCharm:
fe4da5cc 722 func=IpCharm;
723 break;
34f60c01 724 case kBeauty:
fe4da5cc 725 func=IpBeauty;
726 break;
34f60c01 727 case kPion:
753690b0 728 func=IpPion;
729 break;
34f60c01 730 case kKaon:
753690b0 731 func=IpKaon;
732 break;
119b35c7 733 default:
734 func=0;
735 printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
fe4da5cc 736 }
737 return func;
738}
739
740
753690b0 741
05932df6 742Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
743 Float_t dx,
744 Int_t n, Int_t no)
745{
746//
747// Neville's alorithm for interpolation
748//
749// x: x-value
750// y: Input array
751// x0: minimum x
752// dx: step size
753// n: number of data points
754// no: order of polynom
755//
756 Float_t* c = new Float_t[n];
757 Float_t* d = new Float_t[n];
758 Int_t m, i;
759 for (i = 0; i < n; i++) {
760 c[i] = y[i];
761 d[i] = y[i];
762 }
763
764 Int_t ns = int((x - x0)/dx);
765
766 Float_t y1 = y[ns];
767 ns--;
768 for (m = 0; m < no; m++) {
769 for (i = 0; i < n-m; i++) {
770 Float_t ho = x0 + Float_t(i) * dx - x;
771 Float_t hp = x0 + Float_t(i+m+1) * dx - x;
772 Float_t w = c[i+1] - d[i];
773 Float_t den = ho-hp;
774 den = w/den;
775 d[i] = hp * den;
776 c[i] = ho * den;
777 }
778 Float_t dy;
779
780 if (2*ns < (n-m-1)) {
781 dy = c[ns+1];
782 } else {
783 dy = d[ns--];
784 }
785 y1 += dy;}
786 delete[] c;
787 delete[] d;
788
789 return y1;
790}
791
753690b0 792