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
  30 #include "AliGenMUONlib.h"
  31
  32 ClassImp(AliGenMUONlib)
  33 //
  34 //  Pions
  35 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t* /*dummy*/)
  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 //
  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.;
  48   Double_t y, y1, xmpi2, ynorm, a;
  49   Double_t x=*px;
  50   //
  51   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  52   xmpi2=kxmpi*kxmpi;
  53   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
  54   a=ynorm/y1;
  55   if (x > kxlim)
  56     y=a*TMath::Power(kp0/(kp0+x),kxn);
  57   else
  58     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
  59   return y*x;
  60 }
  61 //
  62 // y-distribution
  63 //
  64 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t */*dummy*/)
  65 {
  66 // Pion y
  67   Double_t y=TMath::Abs(*py);
  68 /*
  69   const Double_t ka    = 7000.;
  70   const Double_t kdy   = 4.;
  71   Double_t ex = y*y/(2*kdy*kdy);
  72   return ka*TMath::Exp(-ex);
  73 */
  74   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
  75
  76 }
  77 //                 particle composition
  78 //
  79 Int_t AliGenMUONlib::IpPion(TRandom *ran)
  80 {
  81 // Pion composition
  82     if (ran->Rndm() < 0.5) {
  83         return  211;
  84     } else {
  85         return -211;
  86     }
  87 }
  88
  89 //____________________________________________________________
  90 //
  91 // Mt-scaling
  92
  93 Double_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
  97   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
  98   //     VALUE MESON/PI AT 5 GEV
  99   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 100   np--;
 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];
 103   // PIONS
 104   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 105   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 106                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 107   return fmtscal*ptpion;
 108 }
 109 //
 110 // kaon
 111 //
 112 //                pt-distribution
 113 //____________________________________________________________
 114 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t */*dummy*/)
 115 {
 116 // Kaon pT
 117   return PtScal(*px,2);
 118 }
 119
 120 // y-distribution
 121 //____________________________________________________________
 122 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t */*dummy*/)
 123 {
 124 // Kaon y
 125   Double_t y=TMath::Abs(*py);
 126 /*
 127   const Double_t ka    = 1000.;
 128   const Double_t kdy   = 4.;
 129   //
 130   Double_t ex = y*y/(2*kdy*kdy);
 131   return ka*TMath::Exp(-ex);
 132 */
 133
 134   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 135 }
 136
 137 //                 particle composition
 138 //
 139 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
 140 {
 141 // Kaon composition
 142     if (ran->Rndm() < 0.5) {
 143         return  321;
 144     } else {
 145         return -321;
 146     }
 147 }
 148
 149 //                    J/Psi
 150 //
 151 //
 152 //                pt-distribution
 153 //____________________________________________________________
 154 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t */*dummy*/)
 155 {
 156 // J/Psi pT
 157   const Double_t kpt0 = 4.;
 158   const Double_t kxn  = 3.6;
 159   Double_t x=*px;
 160   //
 161   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 162   return x/TMath::Power(pass1,kxn);
 163 }
 164
 165 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t */*dummy*/)
 166 {
 167 // J/Psi pT spectrum
 168 //
 169 // R. Vogt 2002
 170 // PbPb 5.5 TeV
 171 // MRST HO
 172 // mc = 1.4 GeV, pt-kick 1 GeV
 173 //
 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
 178     };
 179
 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     }
 189     return y;
 190 }
 191 Double_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
 215 //
 216 //               y-distribution
 217 //____________________________________________________________
 218 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t */*dummy*/)
 219 {
 220 // J/psi y
 221   const Double_t ky0 = 4.;
 222   const Double_t kb=1.;
 223   Double_t yj;
 224   Double_t y=TMath::Abs(*py);
 225   //
 226   if (y < ky0)
 227     yj=kb;
 228   else
 229     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 230   return yj;
 231 }
 232
 233
 234 Double_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 //
 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;
 261 }
 262
 263 Double_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
 293 //                 particle composition
 294 //
 295 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 296 {
 297 // J/Psi composition
 298     return 443;
 299 }
 300
 301 //                      Upsilon
 302 //
 303 //
 304 //                  pt-distribution
 305 //____________________________________________________________
 306 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t */*dummy*/ )
 307 {
 308 // Upsilon pT
 309   const Double_t kpt0 = 5.3;
 310   const Double_t kxn  = 2.5;
 311   Double_t x=*px;
 312   //
 313   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 314   return x/TMath::Power(pass1,kxn);
 315 }
 316
 317 Double_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 //
 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     };
 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     }
 343     return y;
 344 }
 345
 346 Double_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
 374 //
 375 //                    y-distribution
 376 //
 377 //____________________________________________________________
 378 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t */*dummy*/)
 379 {
 380 // Upsilon y
 381   const Double_t ky0 = 3.;
 382   const Double_t kb=1.;
 383   Double_t yu;
 384   Double_t y=TMath::Abs(*py);
 385   //
 386   if (y < ky0)
 387     yu=kb;
 388   else
 389     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 390   return yu;
 391 }
 392
 393
 394 Double_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
 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;
 416 }
 417
 418 Double_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
 441 //                 particle composition
 442 //
 443 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 444 {
 445 // y composition
 446     return 553;
 447 }
 448
 449 //
 450 //                        Phi
 451 //
 452 //
 453 //    pt-distribution (by scaling of pion distribution)
 454 //____________________________________________________________
 455 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t */*dummy*/)
 456 {
 457 // Phi pT
 458   return PtScal(*px,7);
 459 }
 460 //    y-distribution
 461 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t */*dummy*/)
 462 {
 463 // Phi y
 464     Double_t *dum=0;
 465     return YJpsi(px,dum);
 466 }
 467 //                 particle composition
 468 //
 469 Int_t AliGenMUONlib::IpPhi(TRandom *)
 470 {
 471 // Phi composition
 472     return 333;
 473 }
 474
 475 //
 476 //                        omega
 477 //
 478 //
 479 //    pt-distribution (by scaling of pion distribution)
 480 //____________________________________________________________
 481 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t */*dummy*/)
 482 {
 483 // Omega pT
 484   return PtScal(*px,5);
 485 }
 486 //    y-distribution
 487 Double_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 //
 495 Int_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 //____________________________________________________________
 508 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t */*dummy*/)
 509 {
 510 // Eta pT
 511   return PtScal(*px,3);
 512 }
 513 //    y-distribution
 514 Double_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 //
 522 Int_t AliGenMUONlib::IpEta(TRandom *)
 523 {
 524 // Eta composition
 525     return 221;
 526 }
 527
 528 //
 529 //                        Charm
 530 //
 531 //
 532 //                    pt-distribution
 533 //____________________________________________________________
 534 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t */*dummy*/)
 535 {
 536 // Charm pT
 537   const Double_t kpt0 = 4.08;
 538   const Double_t kxn  = 9.40;
 539
 540   Double_t x=*px;
 541   //
 542   Double_t pass1 = 1.+(x/kpt0);
 543   return x/TMath::Power(pass1,kxn);
 544 }
 545 //                  y-distribution
 546 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t */*dummy*/)
 547 {
 548 // Charm y
 549     Double_t *dum=0;
 550     return YJpsi(px,dum);
 551 }
 552
 553 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 554 {
 555 // Charm composition
 556     Float_t random;
 557     Int_t ip;
 558 //    411,421,431,4122
 559     random = ran->Rndm();
 560     if (random < 0.5) {
 561         ip=411;
 562     } else if (random < 0.75) {
 563         ip=421;
 564     } else if (random < 0.90) {
 565         ip=431;
 566     } else {
 567         ip=4122;
 568     }
 569     if (ran->Rndm() < 0.5) {ip=-ip;}
 570
 571     return ip;
 572 }
 573
 574
 575 //
 576 //                        Beauty
 577 //
 578 //
 579 //                    pt-distribution
 580 //____________________________________________________________
 581 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t */*dummy*/)
 582 {
 583 // Beauty pT
 584   const Double_t kpt0 = 4.;
 585   const Double_t kxn  = 3.6;
 586   Double_t x=*px;
 587   //
 588   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 589   return x/TMath::Power(pass1,kxn);
 590 }
 591 //                     y-distribution
 592 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t */*dummy*/)
 593 {
 594 // Beauty y
 595     Double_t *dum=0;
 596     return YJpsi(px,dum);
 597 }
 598
 599 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 600 {
 601 // Beauty Composition
 602     Float_t random;
 603     Int_t ip;
 604     random = ran->Rndm();
 605     if (random < 0.5) {
 606         ip=511;
 607     } else if (random < 0.75) {
 608         ip=521;
 609     } else if (random < 0.90) {
 610         ip=531;
 611     } else {
 612         ip=5122;
 613     }
 614     if (ran->Rndm() < 0.5) {ip=-ip;}
 615
 616     return ip;
 617 }
 618
 619 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 620 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 621 {
 622 // Return pointer to pT parameterisation
 623     TString sname = TString(tname);
 624     GenFunc func;
 625     switch (param)
 626     {
 627     case kPhi:
 628         func=PtPhi;
 629         break;
 630     case kOmega:
 631         func=PtOmega;
 632         break;
 633     case kEta:
 634         func=PtEta;
 635         break;
 636     case kJpsi:
 637         if (sname == "Vogt" || sname == "Vogt PbPb") {
 638             func=PtJpsiPbPb;
 639         } else if (sname == "Vogt pp") {
 640             func=PtJpsiPP;
 641         } else {
 642             func=PtJpsi;
 643         }
 644         break;
 645     case kUpsilon:
 646         if (sname == "Vogt" || sname == "Vogt PbPb") {
 647             func=PtUpsilonPbPb;
 648         } else if (sname == "Vogt pp") {
 649             func=PtUpsilonPP;
 650         } else {
 651             func=PtUpsilon;
 652         }
 653         break;
 654     case kCharm:
 655         func=PtCharm;
 656         break;
 657     case kBeauty:
 658         func=PtBeauty;
 659         break;
 660     case kPion:
 661         func=PtPion;
 662         break;
 663     case kKaon:
 664         func=PtKaon;
 665         break;
 666     default:
 667         func=0;
 668         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 669     }
 670     return func;
 671 }
 672
 673 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 674 {
 675     TString sname = TString(tname);
 676
 677 // Return pointer to y- parameterisation
 678     GenFunc func;
 679     switch (param)
 680     {
 681     case kPhi:
 682         func=YPhi;
 683         break;
 684     case kEta:
 685         func=YEta;
 686         break;
 687     case kOmega:
 688         func=YOmega;
 689         break;
 690     case kJpsi:
 691         if (sname == "Vogt" || sname == "Vogt PbPb") {
 692             func=YJpsiPbPb;
 693         } else if (sname == "Vogt pp"){
 694             func=YJpsiPP;
 695         } else {
 696             func=YJpsi;
 697         }
 698
 699         break;
 700     case kUpsilon:
 701         if (sname == "Vogt" || sname == "Vogt PbPb") {
 702             func=YUpsilonPbPb;
 703         } else if (sname == "Vogt pp") {
 704             func = YUpsilonPP;
 705         } else {
 706             func=YUpsilon;
 707         }
 708         break;
 709     case kCharm:
 710         func=YCharm;
 711         break;
 712     case kBeauty:
 713         func=YBeauty;
 714         break;
 715     case kPion:
 716         func=YPion;
 717         break;
 718     case kKaon:
 719         func=YKaon;
 720         break;
 721     default:
 722         func=0;
 723         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 724     }
 725     return func;
 726 }
 727 typedef Int_t (*GenFuncIp) (TRandom *);
 728 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
 729 {
 730 // Return pointer to particle type parameterisation
 731     GenFuncIp func;
 732     switch (param)
 733     {
 734     case kPhi:
 735         func=IpPhi;
 736         break;
 737     case kEta:
 738         func=IpEta;
 739         break;
 740     case kOmega:
 741         func=IpOmega;
 742         break;
 743     case kJpsi:
 744         func=IpJpsi;
 745         break;
 746     case kUpsilon:
 747         func=IpUpsilon;
 748         break;
 749     case kCharm:
 750         func=IpCharm;
 751         break;
 752     case kBeauty:
 753         func=IpBeauty;
 754         break;
 755     case kPion:
 756         func=IpPion;
 757         break;
 758     case kKaon:
 759         func=IpKaon;
 760         break;
 761     default:
 762         func=0;
 763         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
 764     }
 765     return func;
 766 }
 767
 768
 769
 770 Float_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
 820