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