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