Transition to NewIO
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONlib.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
88cb7938 16/* $Id$ */
4c039060 17
53904666 18// Library class for particle pt and y distributions used for
19// muon spectrometer simulations.
20// To be used with AliGenParam.
21// The following particle typed can be simulated:
22// pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
23//
24// andreas.morsch@cern.ch
25//
26
65fb704d 27#include "TMath.h"
28#include "TRandom.h"
29
fe4da5cc 30#include "AliGenMUONlib.h"
5c3fd7ea 31
fe4da5cc 32ClassImp(AliGenMUONlib)
33//
34// Pions
d90f80fd 35Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
fe4da5cc 36{
37//
38// PT-PARAMETERIZATION CDF, PRL 61(88) 1819
39// POWER LAW FOR PT > 500 MEV
40// MT SCALING BELOW (T=160 MEV)
41//
d90f80fd 42 const Double_t kp0 = 1.3;
43 const Double_t kxn = 8.28;
44 const Double_t kxlim=0.5;
45 const Double_t kt=0.160;
46 const Double_t kxmpi=0.139;
47 const Double_t kb=1.;
fe4da5cc 48 Double_t y, y1, xmpi2, ynorm, a;
49 Double_t x=*px;
50 //
d90f80fd 51 y1=TMath::Power(kp0/(kp0+kxlim),kxn);
52 xmpi2=kxmpi*kxmpi;
53 ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
fe4da5cc 54 a=ynorm/y1;
d90f80fd 55 if (x > kxlim)
56 y=a*TMath::Power(kp0/(kp0+x),kxn);
fe4da5cc 57 else
d90f80fd 58 y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
fe4da5cc 59 return y*x;
60}
753690b0 61//
62// y-distribution
63//
d90f80fd 64Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
753690b0 65{
d90f80fd 66// Pion y
2280e6af 67 Double_t y=TMath::Abs(*py);
68/*
d90f80fd 69 const Double_t ka = 7000.;
70 const Double_t kdy = 4.;
d90f80fd 71 Double_t ex = y*y/(2*kdy*kdy);
72 return ka*TMath::Exp(-ex);
2280e6af 73*/
74 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
75
753690b0 76}
77// particle composition
78//
65fb704d 79Int_t AliGenMUONlib::IpPion(TRandom *ran)
753690b0 80{
d90f80fd 81// Pion composition
65fb704d 82 if (ran->Rndm() < 0.5) {
753690b0 83 return 211;
84 } else {
85 return -211;
86 }
87}
fe4da5cc 88
89//____________________________________________________________
90//
91// Mt-scaling
92
93Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
94{
95 // SCALING EN MASSE PAR RAPPORT A PTPI
96 // MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
d90f80fd 97 const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
fe4da5cc 98 // VALUE MESON/PI AT 5 GEV
d90f80fd 99 const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
fe4da5cc 100 np--;
d90f80fd 101 Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
102 Double_t fmax2=f5/kfmax[np];
fe4da5cc 103 // PIONS
104 Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
105 Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
d90f80fd 106 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
fe4da5cc 107 return fmtscal*ptpion;
108}
109//
753690b0 110// kaon
111//
112// pt-distribution
113//____________________________________________________________
d90f80fd 114Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
753690b0 115{
d90f80fd 116// Kaon pT
753690b0 117 return PtScal(*px,2);
118}
119
120// y-distribution
fe4da5cc 121//____________________________________________________________
d90f80fd 122Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
fe4da5cc 123{
d90f80fd 124// Kaon y
2280e6af 125 Double_t y=TMath::Abs(*py);
126/*
d90f80fd 127 const Double_t ka = 1000.;
128 const Double_t kdy = 4.;
fe4da5cc 129 //
d90f80fd 130 Double_t ex = y*y/(2*kdy*kdy);
131 return ka*TMath::Exp(-ex);
2280e6af 132*/
133
134 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
753690b0 135}
136
137// particle composition
138//
65fb704d 139Int_t AliGenMUONlib::IpKaon(TRandom *ran)
753690b0 140{
d90f80fd 141// Kaon composition
65fb704d 142 if (ran->Rndm() < 0.5) {
753690b0 143 return 321;
144 } else {
145 return -321;
146 }
fe4da5cc 147}
753690b0 148
fe4da5cc 149// J/Psi
150//
151//
152// pt-distribution
153//____________________________________________________________
d90f80fd 154Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
fe4da5cc 155{
d90f80fd 156// J/Psi pT
157 const Double_t kpt0 = 4.;
158 const Double_t kxn = 3.6;
fe4da5cc 159 Double_t x=*px;
160 //
d90f80fd 161 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
162 return x/TMath::Power(pass1,kxn);
fe4da5cc 163}
05932df6 164
165Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
166{
1af7144e 167// J/Psi pT spectrum
05932df6 168//
169// R. Vogt 2002
170// PbPb 5.5 TeV
171// MRST HO
172// mc = 1.4 GeV, pt-kick 1 GeV
173//
1af7144e 174 Float_t x = px[0];
175 Float_t c[8] = {
176 -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
177 -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
05932df6 178 };
1af7144e 179
3d905dd7 180 Double_t y;
181 if (x < 10.) {
182 Int_t j;
183 y = c[j = 7];
184 while (j > 0) y = y * x +c[--j];
185 y = x * TMath::Exp(y);
186 } else {
187 y = 0.;
188 }
1af7144e 189 return y;
05932df6 190}
bb6e81ac 191Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t *dummy)
192{
193// J/Psi pT spectrum
194//
195// R. Vogt 2002
196// pp 14 TeV
197// MRST HO
198// mc = 1.4 GeV, pt-kick 1 GeV
199//
200 Float_t x = px[0];
201 Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
202
203 Double_t y;
204 if (x < 10.) {
205 Int_t j;
206 y = c[j = 3];
207 while (j > 0) y = y * x +c[--j];
208 y = x * TMath::Exp(y);
209 } else {
210 y = 0.;
211 }
212 return y;
213}
214
fe4da5cc 215//
216// y-distribution
217//____________________________________________________________
d90f80fd 218Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
fe4da5cc 219{
d90f80fd 220// J/psi y
221 const Double_t ky0 = 4.;
222 const Double_t kb=1.;
fe4da5cc 223 Double_t yj;
224 Double_t y=TMath::Abs(*py);
225 //
d90f80fd 226 if (y < ky0)
227 yj=kb;
fe4da5cc 228 else
d90f80fd 229 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
fe4da5cc 230 return yj;
231}
05932df6 232
233
234Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
235{
236
237//
238// J/Psi y
239//
240//
241// R. Vogt 2002
242// PbPb 5.5 TeV
243// MRST HO
244// mc = 1.4 GeV, pt-kick 1 GeV
245//
1af7144e 246 Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
247 Double_t x = TMath::Abs(px[0]);
248 Double_t y;
249
250 if (x < 4.) {
251 y = 31.754;
252 } else if (x < 6) {
253 Int_t j;
254 y = c[j = 4];
255 while (j > 0) y = y * x + c[--j];
256 } else {
257 y =0.;
258 }
259
260 return y;
05932df6 261}
262
bb6e81ac 263Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t *dummy)
264{
265
266//
267// J/Psi y
268//
269//
270// R. Vogt 2002
271// pp 14 TeV
272// MRST HO
273// mc = 1.4 GeV, pt-kick 1 GeV
274//
275
276 Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
277 Double_t x = TMath::Abs(px[0]);
278 Double_t y;
279
280 if (x < 2.5) {
281 y = 96.455 - 0.8483 * x * x;
282 } else if (x < 7.9) {
283 Int_t j;
284 y = c[j = 4];
285 while (j > 0) y = y * x + c[--j];
286 } else {
287 y =0.;
288 }
289
290 return y;
291}
292
fe4da5cc 293// particle composition
294//
65fb704d 295Int_t AliGenMUONlib::IpJpsi(TRandom *)
fe4da5cc 296{
d90f80fd 297// J/Psi composition
88cb7938 298 return 443;
fe4da5cc 299}
300
301// Upsilon
302//
303//
304// pt-distribution
305//____________________________________________________________
d90f80fd 306Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
fe4da5cc 307{
d90f80fd 308// Upsilon pT
309 const Double_t kpt0 = 5.3;
310 const Double_t kxn = 2.5;
fe4da5cc 311 Double_t x=*px;
312 //
d90f80fd 313 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
314 return x/TMath::Power(pass1,kxn);
fe4da5cc 315}
05932df6 316
317Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
318{
319
320//
321// Upsilon pT
322//
323//
324// R. Vogt 2002
325// PbPb 5.5 TeV
326// MRST HO
327// mc = 1.4 GeV, pt-kick 1 GeV
328//
1af7144e 329 Float_t x = px[0];
330 Double_t c[8] = {
331 -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
332 -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
333 };
3d905dd7 334 Double_t y;
335 if (x < 10.) {
336 Int_t j;
337 y = c[j = 7];
338 while (j > 0) y = y * x +c[--j];
339 y = x * TMath::Exp(y);
340 } else {
341 y = 0.;
342 }
1af7144e 343 return y;
05932df6 344}
345
bb6e81ac 346Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t *dummy)
347{
348
349//
350// Upsilon pT
351//
352//
353// R. Vogt 2002
354// pp 14 TeV
355// MRST HO
356// mc = 1.4 GeV, pt-kick 1 GeV
357//
358 Float_t x = px[0];
359 Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
360 -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
361
362 Double_t y;
363 if (x < 10.) {
364 Int_t j;
365 y = c[j = 7];
366 while (j > 0) y = y * x +c[--j];
367 y = x * TMath::Exp(y);
368 } else {
369 y = 0.;
370 }
371 return y;
372}
373
fe4da5cc 374//
375// y-distribution
376//
377//____________________________________________________________
d90f80fd 378Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
fe4da5cc 379{
d90f80fd 380// Upsilon y
381 const Double_t ky0 = 3.;
382 const Double_t kb=1.;
fe4da5cc 383 Double_t yu;
384 Double_t y=TMath::Abs(*py);
385 //
d90f80fd 386 if (y < ky0)
387 yu=kb;
fe4da5cc 388 else
d90f80fd 389 yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
fe4da5cc 390 return yu;
391}
05932df6 392
393
394Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
395{
396
397//
398// Upsilon y
399//
400//
401// R. Vogt 2002
402// PbPb 5.5 TeV
403// MRST HO
404// mc = 1.4 GeV, pt-kick 1 GeV
405//
406
1af7144e 407 Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
408 -2.99753e-09, 1.28895e-05};
409
410 Double_t x = px[0];
411 if (TMath::Abs(x) > 5.55) return 0.;
412 Int_t j;
413 Double_t y = c[j = 6];
414 while (j > 0) y = y * x +c[--j];
415 return y;
05932df6 416}
417
bb6e81ac 418Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t *dummy)
419{
420
421//
422// Upsilon y
423//
424//
425// R. Vogt 2002
426// p p 14. TeV
427// MRST HO
428// mc = 1.4 GeV, pt-kick 1 GeV
429//
430 Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
431 -6.20539e-10, 1.29943e-05};
432
433 Double_t x = px[0];
434 if (TMath::Abs(x) > 6.2) return 0.;
435 Int_t j;
436 Double_t y = c[j = 6];
437 while (j > 0) y = y * x +c[--j];
438 return y;
439}
440
fe4da5cc 441// particle composition
442//
65fb704d 443Int_t AliGenMUONlib::IpUpsilon(TRandom *)
fe4da5cc 444{
d90f80fd 445// y composition
88cb7938 446 return 553;
fe4da5cc 447}
448
449//
450// Phi
451//
452//
453// pt-distribution (by scaling of pion distribution)
454//____________________________________________________________
d90f80fd 455Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
fe4da5cc 456{
d90f80fd 457// Phi pT
fe4da5cc 458 return PtScal(*px,7);
459}
460// y-distribution
d90f80fd 461Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
fe4da5cc 462{
d90f80fd 463// Phi y
464 Double_t *dum=0;
465 return YJpsi(px,dum);
fe4da5cc 466}
467// particle composition
468//
65fb704d 469Int_t AliGenMUONlib::IpPhi(TRandom *)
fe4da5cc 470{
d90f80fd 471// Phi composition
89512a3b 472 return 333;
473}
474
475//
476// omega
477//
478//
479// pt-distribution (by scaling of pion distribution)
480//____________________________________________________________
481Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
482{
483// Omega pT
484 return PtScal(*px,5);
485}
486// y-distribution
487Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
488{
489// Omega y
490 Double_t *dum=0;
491 return YJpsi(px,dum);
492}
493// particle composition
494//
495Int_t AliGenMUONlib::IpOmega(TRandom *)
496{
497// Omega composition
498 return 223;
499}
500
501
502//
503// Eta
504//
505//
506// pt-distribution (by scaling of pion distribution)
507//____________________________________________________________
508Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
509{
510// Eta pT
511 return PtScal(*px,3);
512}
513// y-distribution
514Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
515{
516// Eta y
517 Double_t *dum=0;
518 return YJpsi(px,dum);
519}
520// particle composition
521//
522Int_t AliGenMUONlib::IpEta(TRandom *)
523{
524// Eta composition
525 return 221;
fe4da5cc 526}
527
528//
529// Charm
530//
531//
532// pt-distribution
533//____________________________________________________________
d90f80fd 534Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
fe4da5cc 535{
d90f80fd 536// Charm pT
537 const Double_t kpt0 = 4.08;
538 const Double_t kxn = 9.40;
2280e6af 539
fe4da5cc 540 Double_t x=*px;
541 //
2280e6af 542 Double_t pass1 = 1.+(x/kpt0);
d90f80fd 543 return x/TMath::Power(pass1,kxn);
fe4da5cc 544}
545// y-distribution
d90f80fd 546Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
fe4da5cc 547{
d90f80fd 548// Charm y
549 Double_t *dum=0;
550 return YJpsi(px,dum);
fe4da5cc 551}
552
65fb704d 553Int_t AliGenMUONlib::IpCharm(TRandom *ran)
fe4da5cc 554{
d90f80fd 555// Charm composition
65fb704d 556 Float_t random;
fe4da5cc 557 Int_t ip;
558// 411,421,431,4122
65fb704d 559 random = ran->Rndm();
560 if (random < 0.5) {
fe4da5cc 561 ip=411;
65fb704d 562 } else if (random < 0.75) {
fe4da5cc 563 ip=421;
65fb704d 564 } else if (random < 0.90) {
fe4da5cc 565 ip=431;
566 } else {
567 ip=4122;
568 }
65fb704d 569 if (ran->Rndm() < 0.5) {ip=-ip;}
fe4da5cc 570
571 return ip;
572}
573
574
575//
576// Beauty
577//
578//
579// pt-distribution
580//____________________________________________________________
d90f80fd 581Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
fe4da5cc 582{
d90f80fd 583// Beauty pT
584 const Double_t kpt0 = 4.;
585 const Double_t kxn = 3.6;
fe4da5cc 586 Double_t x=*px;
587 //
d90f80fd 588 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
589 return x/TMath::Power(pass1,kxn);
fe4da5cc 590}
591// y-distribution
d90f80fd 592Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
fe4da5cc 593{
d90f80fd 594// Beauty y
595 Double_t *dum=0;
596 return YJpsi(px,dum);
fe4da5cc 597}
598
65fb704d 599Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
fe4da5cc 600{
d90f80fd 601// Beauty Composition
65fb704d 602 Float_t random;
fe4da5cc 603 Int_t ip;
65fb704d 604 random = ran->Rndm();
605 if (random < 0.5) {
fe4da5cc 606 ip=511;
65fb704d 607 } else if (random < 0.75) {
fe4da5cc 608 ip=521;
65fb704d 609 } else if (random < 0.90) {
fe4da5cc 610 ip=531;
611 } else {
612 ip=5122;
613 }
65fb704d 614 if (ran->Rndm() < 0.5) {ip=-ip;}
fe4da5cc 615
616 return ip;
617}
618
619typedef Double_t (*GenFunc) (Double_t*, Double_t*);
53904666 620GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname) const
fe4da5cc 621{
d90f80fd 622// Return pointer to pT parameterisation
05932df6 623 TString sname = TString(tname);
fe4da5cc 624 GenFunc func;
753690b0 625 switch (param)
fe4da5cc 626 {
34f60c01 627 case kPhi:
fe4da5cc 628 func=PtPhi;
629 break;
89512a3b 630 case kOmega:
631 func=PtOmega;
632 break;
633 case kEta:
634 func=PtEta;
635 break;
34f60c01 636 case kJpsi:
bb6e81ac 637 if (sname == "Vogt" || sname == "Vogt PbPb") {
05932df6 638 func=PtJpsiPbPb;
bb6e81ac 639 } else if (sname == "Vogt pp") {
640 func=PtJpsiPP;
05932df6 641 } else {
642 func=PtJpsi;
643 }
fe4da5cc 644 break;
34f60c01 645 case kUpsilon:
bb6e81ac 646 if (sname == "Vogt" || sname == "Vogt PbPb") {
05932df6 647 func=PtUpsilonPbPb;
bb6e81ac 648 } else if (sname == "Vogt pp") {
649 func=PtUpsilonPP;
05932df6 650 } else {
651 func=PtUpsilon;
652 }
fe4da5cc 653 break;
34f60c01 654 case kCharm:
fe4da5cc 655 func=PtCharm;
656 break;
34f60c01 657 case kBeauty:
fe4da5cc 658 func=PtBeauty;
659 break;
34f60c01 660 case kPion:
753690b0 661 func=PtPion;
662 break;
34f60c01 663 case kKaon:
753690b0 664 func=PtKaon;
665 break;
119b35c7 666 default:
667 func=0;
668 printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
fe4da5cc 669 }
670 return func;
671}
672
53904666 673GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
fe4da5cc 674{
05932df6 675 TString sname = TString(tname);
676
d90f80fd 677// Return pointer to y- parameterisation
fe4da5cc 678 GenFunc func;
753690b0 679 switch (param)
fe4da5cc 680 {
34f60c01 681 case kPhi:
fe4da5cc 682 func=YPhi;
683 break;
89512a3b 684 case kEta:
685 func=YEta;
686 break;
687 case kOmega:
688 func=YOmega;
689 break;
34f60c01 690 case kJpsi:
bb6e81ac 691 if (sname == "Vogt" || sname == "Vogt PbPb") {
05932df6 692 func=YJpsiPbPb;
bb6e81ac 693 } else if (sname == "Vogt pp"){
694 func=YJpsiPP;
05932df6 695 } else {
696 func=YJpsi;
697 }
bb6e81ac 698
fe4da5cc 699 break;
34f60c01 700 case kUpsilon:
bb6e81ac 701 if (sname == "Vogt" || sname == "Vogt PbPb") {
05932df6 702 func=YUpsilonPbPb;
bb6e81ac 703 } else if (sname == "Vogt pp") {
704 func = YUpsilonPP;
05932df6 705 } else {
706 func=YUpsilon;
707 }
fe4da5cc 708 break;
34f60c01 709 case kCharm:
fe4da5cc 710 func=YCharm;
711 break;
34f60c01 712 case kBeauty:
fe4da5cc 713 func=YBeauty;
714 break;
34f60c01 715 case kPion:
753690b0 716 func=YPion;
717 break;
34f60c01 718 case kKaon:
753690b0 719 func=YKaon;
720 break;
119b35c7 721 default:
722 func=0;
723 printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
fe4da5cc 724 }
725 return func;
726}
65fb704d 727typedef Int_t (*GenFuncIp) (TRandom *);
53904666 728GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* tname) const
fe4da5cc 729{
d90f80fd 730// Return pointer to particle type parameterisation
fe4da5cc 731 GenFuncIp func;
753690b0 732 switch (param)
fe4da5cc 733 {
34f60c01 734 case kPhi:
fe4da5cc 735 func=IpPhi;
736 break;
89512a3b 737 case kEta:
738 func=IpEta;
739 break;
740 case kOmega:
741 func=IpOmega;
742 break;
34f60c01 743 case kJpsi:
fe4da5cc 744 func=IpJpsi;
745 break;
34f60c01 746 case kUpsilon:
fe4da5cc 747 func=IpUpsilon;
748 break;
34f60c01 749 case kCharm:
fe4da5cc 750 func=IpCharm;
751 break;
34f60c01 752 case kBeauty:
fe4da5cc 753 func=IpBeauty;
754 break;
34f60c01 755 case kPion:
753690b0 756 func=IpPion;
757 break;
34f60c01 758 case kKaon:
753690b0 759 func=IpKaon;
760 break;
119b35c7 761 default:
762 func=0;
763 printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
fe4da5cc 764 }
765 return func;
766}
767
768
753690b0 769
05932df6 770Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
771 Float_t dx,
772 Int_t n, Int_t no)
773{
774//
775// Neville's alorithm for interpolation
776//
777// x: x-value
778// y: Input array
779// x0: minimum x
780// dx: step size
781// n: number of data points
782// no: order of polynom
783//
784 Float_t* c = new Float_t[n];
785 Float_t* d = new Float_t[n];
786 Int_t m, i;
787 for (i = 0; i < n; i++) {
788 c[i] = y[i];
789 d[i] = y[i];
790 }
791
792 Int_t ns = int((x - x0)/dx);
793
794 Float_t y1 = y[ns];
795 ns--;
796 for (m = 0; m < no; m++) {
797 for (i = 0; i < n-m; i++) {
798 Float_t ho = x0 + Float_t(i) * dx - x;
799 Float_t hp = x0 + Float_t(i+m+1) * dx - x;
800 Float_t w = c[i+1] - d[i];
801 Float_t den = ho-hp;
802 den = w/den;
803 d[i] = hp * den;
804 c[i] = ho * den;
805 }
806 Float_t dy;
807
808 if (2*ns < (n-m-1)) {
809 dy = c[ns+1];
810 } else {
811 dy = d[ns--];
812 }
813 y1 += dy;}
814 delete[] c;
815 delete[] d;
816
817 return y1;
818}
819
753690b0 820