Limited pT range for parameterized Upsilon and J/Psi pT distributions.
[u/mrichter/AliRoot.git] / 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 /*
17 $Log$
18 Revision 1.19  2003/02/24 16:46:11  morsch
19 New parameterisation for Psi and Upsilon (PbPb)
20
21 Revision 1.18  2002/11/07 09:13:42  morsch
22 Use "Vogt" to label new distributions.
23
24 Revision 1.17  2002/11/07 09:06:10  morsch
25 J/Psi and Upsilon pt and y distributions from R. Vogt 2002 added.
26
27 Revision 1.16  2002/10/14 14:55:35  hristov
28 Merging the VirtualMC branch to the main development branch (HEAD)
29
30 Revision 1.14.6.1  2002/06/10 14:57:41  hristov
31 Merged with v3-08-02
32
33 Revision 1.15  2002/04/17 10:11:51  morsch
34 Coding Rule violations corrected.
35
36 Revision 1.14  2002/02/22 17:26:43  morsch
37 Eta and omega added.
38
39 Revision 1.13  2001/03/27 11:01:04  morsch
40 Charm pt-distribution corrected. More realistic y-distribution for pi and K.
41
42 Revision 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
46 Revision 1.11  2000/11/30 07:12:50  alibrary
47 Introducing new Rndm and QA classes
48
49 Revision 1.10  2000/06/29 21:08:27  morsch
50 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
51 This allows to pass a pointer to a library directly to AliGenParam and avoids the
52 use of function pointers in Config.C.
53
54 Revision 1.9  2000/06/14 15:20:56  morsch
55 Include clean-up (IH)
56
57 Revision 1.8  2000/06/09 20:32:11  morsch
58 All coding rule violations except RS3 corrected
59
60 Revision 1.7  2000/05/02 08:12:13  morsch
61 Coding rule violations corrected.
62
63 Revision 1.6  1999/09/29 09:24:14  fca
64 Introduction of the Copyright and cvs Log
65
66 */
67
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
77 #include "TMath.h"
78 #include "TRandom.h"
79
80 #include "AliGenMUONlib.h"
81
82 ClassImp(AliGenMUONlib)
83 //
84 //  Pions
85 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
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 //
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.;
98   Double_t y, y1, xmpi2, ynorm, a;
99   Double_t x=*px;
100   //
101   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
102   xmpi2=kxmpi*kxmpi;
103   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
104   a=ynorm/y1;
105   if (x > kxlim)
106     y=a*TMath::Power(kp0/(kp0+x),kxn);
107   else
108     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
109   return y*x;
110 }
111 //
112 // y-distribution
113 //
114 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
115 {
116 // Pion y
117   Double_t y=TMath::Abs(*py);
118 /*
119   const Double_t ka    = 7000.;
120   const Double_t kdy   = 4.;
121   Double_t ex = y*y/(2*kdy*kdy);
122   return ka*TMath::Exp(-ex);
123 */
124   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
125   
126 }
127 //                 particle composition
128 //
129 Int_t AliGenMUONlib::IpPion(TRandom *ran)
130 {
131 // Pion composition 
132     if (ran->Rndm() < 0.5) {
133         return  211;
134     } else {
135         return -211;
136     }
137 }
138
139 //____________________________________________________________
140 //
141 // Mt-scaling
142
143 Double_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
147   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
148   //     VALUE MESON/PI AT 5 GEV
149   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
150   np--;
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];
153   // PIONS
154   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
155   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
156                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
157   return fmtscal*ptpion;
158 }
159 //
160 // kaon
161 //
162 //                pt-distribution
163 //____________________________________________________________
164 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
165 {
166 // Kaon pT
167   return PtScal(*px,2);
168 }
169
170 // y-distribution
171 //____________________________________________________________
172 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
173 {
174 // Kaon y
175   Double_t y=TMath::Abs(*py);
176 /*
177   const Double_t ka    = 1000.;
178   const Double_t kdy   = 4.;
179   //
180   Double_t ex = y*y/(2*kdy*kdy);
181   return ka*TMath::Exp(-ex);
182 */
183
184   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
185 }
186
187 //                 particle composition
188 //
189 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
190 {
191 // Kaon composition
192     if (ran->Rndm() < 0.5) {
193         return  321;
194     } else {
195         return -321;
196     }
197 }
198
199 //                    J/Psi 
200 //
201 //
202 //                pt-distribution
203 //____________________________________________________________
204 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
205 {
206 // J/Psi pT
207   const Double_t kpt0 = 4.;
208   const Double_t kxn  = 3.6;
209   Double_t x=*px;
210   //
211   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
212   return x/TMath::Power(pass1,kxn);
213 }
214
215 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
216 {
217 // J/Psi pT spectrum
218 //
219 // R. Vogt 2002
220 // PbPb 5.5 TeV
221 // MRST HO
222 // mc = 1.4 GeV, pt-kick 1 GeV
223 //
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
228     };
229     
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     }
239     return y;
240 }
241 //
242 //               y-distribution
243 //____________________________________________________________
244 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
245 {
246 // J/psi y
247   const Double_t ky0 = 4.;
248   const Double_t kb=1.;
249   Double_t yj;
250   Double_t y=TMath::Abs(*py);
251   //
252   if (y < ky0)
253     yj=kb;
254   else
255     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
256   return yj;
257 }
258
259
260 Double_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 //
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;
287 }
288
289 //                 particle composition
290 //
291 Int_t AliGenMUONlib::IpJpsi(TRandom *)
292 {
293 // J/Psi composition
294     return 443;
295 }
296
297 //                      Upsilon
298 //
299 //
300 //                  pt-distribution
301 //____________________________________________________________
302 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
303 {
304 // Upsilon pT
305   const Double_t kpt0 = 5.3;
306   const Double_t kxn  = 2.5;
307   Double_t x=*px;
308   //
309   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
310   return x/TMath::Power(pass1,kxn);
311 }
312
313 Double_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 //
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     };
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     }
339     return y;
340 }
341
342 //
343 //                    y-distribution
344 //
345 //____________________________________________________________
346 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
347 {
348 // Upsilon y
349   const Double_t ky0 = 3.;
350   const Double_t kb=1.;
351   Double_t yu;
352   Double_t y=TMath::Abs(*py);
353   //
354   if (y < ky0)
355     yu=kb;
356   else
357     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
358   return yu;
359 }
360
361
362 Double_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
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;
384 }
385
386 //                 particle composition
387 //
388 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
389 {
390 // y composition
391     return 553;
392 }
393
394 //
395 //                        Phi
396 //
397 //
398 //    pt-distribution (by scaling of pion distribution)
399 //____________________________________________________________
400 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
401 {
402 // Phi pT
403   return PtScal(*px,7);
404 }
405 //    y-distribution
406 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
407 {
408 // Phi y
409     Double_t *dum=0;
410     return YJpsi(px,dum);
411 }
412 //                 particle composition
413 //
414 Int_t AliGenMUONlib::IpPhi(TRandom *)
415 {
416 // Phi composition
417     return 333;
418 }
419
420 //
421 //                        omega
422 //
423 //
424 //    pt-distribution (by scaling of pion distribution)
425 //____________________________________________________________
426 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
427 {
428 // Omega pT
429   return PtScal(*px,5);
430 }
431 //    y-distribution
432 Double_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 //
440 Int_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 //____________________________________________________________
453 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
454 {
455 // Eta pT
456   return PtScal(*px,3);
457 }
458 //    y-distribution
459 Double_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 //
467 Int_t AliGenMUONlib::IpEta(TRandom *)
468 {
469 // Eta composition
470     return 221;
471 }
472
473 //
474 //                        Charm
475 //
476 //
477 //                    pt-distribution
478 //____________________________________________________________
479 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
480 {
481 // Charm pT
482   const Double_t kpt0 = 4.08;
483   const Double_t kxn  = 9.40;
484
485   Double_t x=*px;
486   //
487   Double_t pass1 = 1.+(x/kpt0);
488   return x/TMath::Power(pass1,kxn);
489 }
490 //                  y-distribution
491 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
492 {
493 // Charm y
494     Double_t *dum=0;
495     return YJpsi(px,dum);
496 }
497
498 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
499 {  
500 // Charm composition
501     Float_t random;
502     Int_t ip;
503 //    411,421,431,4122
504     random = ran->Rndm();
505     if (random < 0.5) {
506         ip=411;
507     } else if (random < 0.75) {
508         ip=421;
509     } else if (random < 0.90) {
510         ip=431;
511     } else {
512         ip=4122;
513     }
514     if (ran->Rndm() < 0.5) {ip=-ip;}
515     
516     return ip;
517 }
518
519
520 //
521 //                        Beauty
522 //
523 //
524 //                    pt-distribution
525 //____________________________________________________________
526 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
527 {
528 // Beauty pT
529   const Double_t kpt0 = 4.;
530   const Double_t kxn  = 3.6;
531   Double_t x=*px;
532   //
533   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
534   return x/TMath::Power(pass1,kxn);
535 }
536 //                     y-distribution
537 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
538 {
539 // Beauty y
540     Double_t *dum=0;
541     return YJpsi(px,dum);
542 }
543
544 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
545 {  
546 // Beauty Composition
547     Float_t random;
548     Int_t ip;
549     random = ran->Rndm();
550     if (random < 0.5) {
551         ip=511;
552     } else if (random < 0.75) {
553         ip=521;
554     } else if (random < 0.90) {
555         ip=531;
556     } else {
557         ip=5122;
558     }
559     if (ran->Rndm() < 0.5) {ip=-ip;}
560     
561     return ip;
562 }
563
564 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
565 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
566 {
567 // Return pointer to pT parameterisation
568     TString sname = TString(tname);
569     GenFunc func;
570     switch (param) 
571     {
572     case kPhi:
573         func=PtPhi;
574         break;
575     case kOmega:
576         func=PtOmega;
577         break;
578     case kEta:
579         func=PtEta;
580         break;
581     case kJpsi:
582         if (sname == "Vogt") {
583             func=PtJpsiPbPb;
584         } else {
585             func=PtJpsi;
586         }
587         break;
588     case kUpsilon:
589         if (sname == "Vogt") {
590             func=PtUpsilonPbPb;
591         } else {
592             func=PtUpsilon;
593         }
594         break;
595     case kCharm:
596         func=PtCharm;
597         break;
598     case kBeauty:
599         func=PtBeauty;
600         break;
601     case kPion:
602         func=PtPion;
603         break;
604     case kKaon:
605         func=PtKaon;
606         break;
607     default:
608         func=0;
609         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
610     }
611     return func;
612 }
613
614 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
615 {
616     TString sname = TString(tname);
617     
618 // Return pointer to y- parameterisation
619     GenFunc func;
620     switch (param) 
621     {
622     case kPhi:
623         func=YPhi;
624         break;
625     case kEta:
626         func=YEta;
627         break;
628     case kOmega:
629         func=YOmega;
630         break;
631     case kJpsi:
632         if (sname == "Vogt") {
633             func=YJpsiPbPb;
634         } else {
635             func=YJpsi;
636         }
637         break;
638     case kUpsilon:
639         if (sname == "Vogt") {
640             func=YUpsilonPbPb;
641         } else {
642             func=YUpsilon;
643         }
644         break;
645     case kCharm:
646         func=YCharm;
647         break;
648     case kBeauty:
649         func=YBeauty;
650         break;
651     case kPion:
652         func=YPion;
653         break;
654     case kKaon:
655         func=YKaon;
656         break;
657     default:
658         func=0;
659         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
660     }
661     return func;
662 }
663 typedef Int_t (*GenFuncIp) (TRandom *);
664 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
665 {
666 // Return pointer to particle type parameterisation
667     GenFuncIp func;
668     switch (param) 
669     {
670     case kPhi:
671         func=IpPhi;
672         break;
673     case kEta:
674         func=IpEta;
675         break;
676     case kOmega:
677         func=IpOmega;
678         break;
679     case kJpsi:
680         func=IpJpsi;
681         break;
682     case kUpsilon:
683         func=IpUpsilon;
684         break;
685     case kCharm:
686         func=IpCharm;
687         break;
688     case kBeauty:
689         func=IpBeauty;
690         break;
691     case kPion:
692         func=IpPion;
693         break;
694     case kKaon:
695         func=IpKaon;
696         break;
697     default:
698         func=0;
699         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
700     }
701     return func;
702 }
703
704
705
706 Float_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
756