small update
[u/mrichter/AliRoot.git] / EVGEN / AliGenPHOSlib.cxx
CommitLineData
886b6f73 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$ */
886b6f73 17
18//======================================================================
19// AliGenPHOSlib class contains parameterizations of the
20// pion, kaon, eta, omega, etaprime, phi and baryon (proton,
21// antiproton, neutron and anti-neutron) particles for the
22// study of the neutral background in PHOS detector.
23// These parameterizations are used by the
24// AliGenParam class:
25// AliGenParam(npar, param, AliGenPHOSlib::GetPt(param),
26// AliGenPHOSlib::GetY(param),
27// AliGenPHOSlib::GetIp(param) )
28// param represents the particle to be simulated :
29// Pion, Kaon, Eta, Omega, Etaprime, Phi or Baryon
30// Pt distributions are calculated from the transverse mass scaling
9ae59e17 31// with Pions, using the PtScal function taken from AliGenMUONlib
886b6f73 32// version aliroot 3.01
33//
9ae59e17 34// Gines MARTINEZ. Laurent APHECETCHE and Yves SCHUTZ
35// GPS @ SUBATECH, Nantes , France (October 1999)
886b6f73 36// http://www-subatech.in2p3.fr/~photons/subatech
37// martinez@subatech.in2p3.fr
35ead6c1 38// Additional particle species simulation options has been added:
39// Charged Pion, Charged Kaons, KLong Proton, Anti-Proton, Neutron,
40// Anti-Neutron --> Changes made by Gustavo Conesa in November 2004
8fc6910b 41// Add flat Omega(782) distribution in Nov. 2010 by Renzhuo WAN
886b6f73 42//======================================================================
43
65fb704d 44#include "TMath.h"
45#include "TRandom.h"
46
886b6f73 47#include "AliGenPHOSlib.h"
886b6f73 48
49ClassImp(AliGenPHOSlib)
50
51//======================================================================
52// P I O N S
53// (From GetPt, GetY and GetIp as param = Pion)
54// Transverse momentum distribution" PtPion
55// Rapidity distribution YPion
56// Particle distribution IdPion 111, 211 and -211 (pi0, pi+ and pi-)
57//
75e0cc59 58 Double_t AliGenPHOSlib::PtPion(const Double_t *px, const Double_t *)
886b6f73 59{
60// Pion transverse momentum distribtuion taken
9ae59e17 61// from AliGenMUONlib class, version 3.01 of aliroot
886b6f73 62// PT-PARAMETERIZATION CDF, PRL 61(88) 1819
63// POWER LAW FOR PT > 500 MEV
64// MT SCALING BELOW (T=160 MEV)
65//
f87cfe57 66 const Double_t kp0 = 1.3;
67 const Double_t kxn = 8.28;
68 const Double_t kxlim=0.5;
69 const Double_t kt=0.160;
70 const Double_t kxmpi=0.139;
71 const Double_t kb=1.;
72 Double_t y, y1, kxmpi2, ynorm, a;
886b6f73 73 Double_t x=*px;
74 //
f87cfe57 75 y1=TMath::Power(kp0/(kp0+kxlim),kxn);
76 kxmpi2=kxmpi*kxmpi;
77 ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+kxmpi2)/kt));
886b6f73 78 a=ynorm/y1;
f87cfe57 79 if (x > kxlim)
80 y=a*TMath::Power(kp0/(kp0+x),kxn);
886b6f73 81 else
f87cfe57 82 y=kb*TMath::Exp(-sqrt(x*x+kxmpi2)/kt);
886b6f73 83 return y*x;
84}
75e0cc59 85 Double_t AliGenPHOSlib::YPion( const Double_t *py, const Double_t *)
886b6f73 86{
f87cfe57 87//
88// pion y-distribution
89//
90
91 const Double_t ka = 7000.;
92 const Double_t kdy = 4.;
886b6f73 93
94 Double_t y=TMath::Abs(*py);
95 //
f87cfe57 96 Double_t ex = y*y/(2*kdy*kdy);
97 return ka*TMath::Exp(-ex);
886b6f73 98}
f87cfe57 99
078ff23f 100Int_t AliGenPHOSlib::IpPion(TRandom */*ran*/)
886b6f73 101{
f87cfe57 102// particle composition pi+, pi0, pi-
103//
104
886b6f73 105 return 111 ;
886b6f73 106}
35ead6c1 107 Int_t AliGenPHOSlib::IpChargedPion(TRandom *ran)
108{
109// particle composition pi+, pi0, pi-
110//
111
112 Float_t random = ran->Rndm();
113
114 if ( (2.*random) < 1. )
115 {
116 return 211 ;
117 }
118 else
119 {
120 return -211 ;
121 }
122}
74c62c73 123
124//End Pions
125//======================================================================
126// Pi 0 Flat Distribution
127// Transverse momentum distribution PtPi0Flat
128// Rapidity distribution YPi0Flat
129// Particle distribution IdPi0Flat 111 (pi0)
130//
131
6784ce25 132Double_t AliGenPHOSlib::PtPi0(const Double_t * px, const Double_t *)
133{
134// Pion transverse momentum
135 const Double_t kp0 =1.35;
136 const Double_t kxn= 6.18;
137 return TMath::Power(kp0 /(kp0 + px[0]), kxn);
138}
139
75e0cc59 140Double_t AliGenPHOSlib::PtPi0Flat(const Double_t */*px*/, const Double_t *)
74c62c73 141{
142// Pion transverse momentum flat distribution
143
144return 1;
145
146}
147
75e0cc59 148Double_t AliGenPHOSlib::YPi0Flat( const Double_t */*py*/, const Double_t *)
74c62c73 149{
150
151// pion y-distribution
152//
153 return 1.;
154}
155
156 Int_t AliGenPHOSlib::IpPi0Flat(TRandom *)
157{
158
159// particle composition pi0
160//
161 return 111 ;
162}
163// End Pi0Flat
886b6f73 164//=============================================================
165//
f87cfe57 166 Double_t AliGenPHOSlib::PtScal(Double_t pt, Int_t np)
167{
886b6f73 168// Mt-scaling
169// Fonction for the calculation of the Pt distribution for a
170// given particle np, from the pion Pt distribution using the
171// mt scaling. This function was taken from AliGenMUONlib
172// aliroot version 3.01, and was extended for baryons
173// np = 1=>Pions 2=>Kaons 3=>Etas 4=>Omegas 5=>ETA' 6=>PHI
174// 7=>BARYONS-BARYONBARS
f87cfe57 175
886b6f73 176 // SCALING EN MASSE PAR RAPPORT A PTPI
0db2f441 177 // MASS 0=>PI, 1=>K, 2=>ETA, 3=>OMEGA, 4=>ETA',5=>PHI
f87cfe57 178 const Double_t khm[10] = {0.1396, 0.494, 0.547, 0.782, 0.957, 1.02,
0db2f441 179 // MASS 6=>BARYON-BARYONBAR
886b6f73 180 0.938, 0. , 0., 0.};
181 // VALUE MESON/PI AT 5 GEV
f87cfe57 182 const Double_t kfmax[10]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
f87cfe57 183 Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
184 Double_t kfmax2=f5/kfmax[np];
886b6f73 185 // PIONS
186 Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
187 Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
f87cfe57 188 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
886b6f73 189 return fmtscal*ptpion;
74c62c73 190
886b6f73 191}
192// End Scaling
193//============================================================================
194// K A O N S
75e0cc59 195 Double_t AliGenPHOSlib::PtKaon( const Double_t *px, const Double_t *)
886b6f73 196{
f87cfe57 197// kaon
198// pt-distribution
199//____________________________________________________________
200
0db2f441 201 return PtScal(*px,1); // 1==> Kaon in the PtScal function
886b6f73 202}
203
75e0cc59 204 Double_t AliGenPHOSlib::YKaon( const Double_t *py, const Double_t *)
886b6f73 205{
f87cfe57 206// y-distribution
207//____________________________________________________________
208
209 const Double_t ka = 1000.;
210 const Double_t kdy = 4.;
886b6f73 211
212
213 Double_t y=TMath::Abs(*py);
214 //
f87cfe57 215 Double_t ex = y*y/(2*kdy*kdy);
216 return ka*TMath::Exp(-ex);
886b6f73 217}
218
65fb704d 219 Int_t AliGenPHOSlib::IpKaon(TRandom *ran)
886b6f73 220{
f87cfe57 221// particle composition
222//
223
65fb704d 224 Float_t random = ran->Rndm();
225 Float_t random2 = ran->Rndm();
226 if (random2 < 0.5)
886b6f73 227 {
65fb704d 228 if (random < 0.5) {
886b6f73 229 return 321; // K+
230 } else {
231 return -321; // K-
232 }
233 }
234 else
235 {
65fb704d 236 if (random < 0.5) {
9ae59e17 237 return 130; // K^0 short
886b6f73 238 } else {
9ae59e17 239 return 310; // K^0 long
886b6f73 240 }
241 }
242}
35ead6c1 243
244 Int_t AliGenPHOSlib::IpChargedKaon(TRandom *ran)
245{
246// particle composition
247//
248
249 Float_t random = ran->Rndm();
250
251 if (random < 0.5) {
252 return 321; // K+
253 } else {
254 return -321; // K-
255 }
256
257
258}
259Int_t AliGenPHOSlib::IpKaon0L(TRandom *)
260{
261 // particle composition
262 //
263
264 return 130; // K^0 long
265}
886b6f73 266// End Kaons
267//============================================================================
268//============================================================================
269// E T A S
75e0cc59 270 Double_t AliGenPHOSlib::PtEta( const Double_t *px, const Double_t *)
886b6f73 271{
f87cfe57 272// etas
273// pt-distribution
274//____________________________________________________________
275
0db2f441 276 return PtScal(*px,2); // 2==> Eta in the PtScal function
886b6f73 277}
278
75e0cc59 279 Double_t AliGenPHOSlib::YEta( const Double_t *py, const Double_t *)
886b6f73 280{
f87cfe57 281// y-distribution
282//____________________________________________________________
283
284 const Double_t ka = 1000.;
285 const Double_t kdy = 4.;
886b6f73 286
287
288 Double_t y=TMath::Abs(*py);
289 //
f87cfe57 290 Double_t ex = y*y/(2*kdy*kdy);
291 return ka*TMath::Exp(-ex);
886b6f73 292}
293
65fb704d 294 Int_t AliGenPHOSlib::IpEta(TRandom *)
886b6f73 295{
f87cfe57 296// particle composition
297//
298
886b6f73 299 return 221; // eta
300}
301// End Etas
74c62c73 302
303//======================================================================
304// Eta Flat Distribution
305// Transverse momentum distribution PtEtaFlat
306// Rapidity distribution YEtaFlat
307// Particle distribution IdEtaFlat 111 (pi0)
308//
309
75e0cc59 310Double_t AliGenPHOSlib::PtEtaFlat(const Double_t */*px*/, const Double_t *)
74c62c73 311{
312// Eta transverse momentum flat distribution
313
314 return 1;
315
316}
317
75e0cc59 318Double_t AliGenPHOSlib::YEtaFlat( const Double_t */*py*/, const Double_t *)
74c62c73 319{
320//
321// pion y-distribution
322//
323 return 1.;
324}
325
326 Int_t AliGenPHOSlib::IpEtaFlat(TRandom *)
327{
328//
329// particle composition eta
330//
331 return 221 ;
332}
35ead6c1 333// End EtaFlat
886b6f73 334//============================================================================
335//============================================================================
336// O M E G A S
75e0cc59 337 Double_t AliGenPHOSlib::PtOmega( const Double_t *px, const Double_t *)
f87cfe57 338{
886b6f73 339// omegas
340// pt-distribution
341//____________________________________________________________
f87cfe57 342
0db2f441 343 return PtScal(*px,3); // 3==> Omega in the PtScal function
886b6f73 344}
345
75e0cc59 346 Double_t AliGenPHOSlib::YOmega( const Double_t *py, const Double_t *)
886b6f73 347{
f87cfe57 348// y-distribution
349//____________________________________________________________
350
351 const Double_t ka = 1000.;
352 const Double_t kdy = 4.;
886b6f73 353
354
355 Double_t y=TMath::Abs(*py);
356 //
f87cfe57 357 Double_t ex = y*y/(2*kdy*kdy);
358 return ka*TMath::Exp(-ex);
886b6f73 359}
360
65fb704d 361 Int_t AliGenPHOSlib::IpOmega(TRandom *)
886b6f73 362{
f87cfe57 363// particle composition
364//
365
886b6f73 366 return 223; // Omega
367}
368// End Omega
369//============================================================================
8fc6910b 370//======================================================================
371// Omega(782) Flat Distribution
372// Transverse momentum distribution PtOmegaFlat
373// Rapidity distribution YOmegaFlat
374// Particle distribution IdOmegaFlat 223(0mega)
375//
376
377Double_t AliGenPHOSlib::PtOmegaFlat(const Double_t */*px*/, const Double_t *)
378{
379// omega transverse momentum flat distribution
380
381return 1;
382
383}
384
385Double_t AliGenPHOSlib::YOmegaFlat( const Double_t */*py*/, const Double_t *)
386{
387
388// omega y-distribution
389//
390 return 1.;
391}
392
393 Int_t AliGenPHOSlib::IpOmegaFlat(TRandom *)
394{
395
396// particle composition omega
397//
398 return 223 ;
399}
400// End OmegaFlat
401
402
886b6f73 403//============================================================================
404// E T A P R I M E
75e0cc59 405 Double_t AliGenPHOSlib::PtEtaprime( const Double_t *px, const Double_t *)
f87cfe57 406{
886b6f73 407// etaprime
408// pt-distribution
409//____________________________________________________________
f87cfe57 410
0db2f441 411 return PtScal(*px,4); // 4==> Etaprime in the PtScal function
886b6f73 412}
413
75e0cc59 414 Double_t AliGenPHOSlib::YEtaprime( const Double_t *py, const Double_t *)
886b6f73 415{
f87cfe57 416// y-distribution
417//____________________________________________________________
418
419 const Double_t ka = 1000.;
420 const Double_t kdy = 4.;
886b6f73 421
422
423 Double_t y=TMath::Abs(*py);
424 //
f87cfe57 425 Double_t ex = y*y/(2*kdy*kdy);
426 return ka*TMath::Exp(-ex);
886b6f73 427}
428
65fb704d 429 Int_t AliGenPHOSlib::IpEtaprime(TRandom *)
886b6f73 430{
f87cfe57 431// particle composition
432//
433
886b6f73 434 return 331; // Etaprime
435}
436// End EtaPrime
437//===================================================================
438//============================================================================
439// P H I S
75e0cc59 440 Double_t AliGenPHOSlib::PtPhi( const Double_t *px, const Double_t *)
f87cfe57 441{
886b6f73 442// phi
443// pt-distribution
444//____________________________________________________________
f87cfe57 445
0db2f441 446 return PtScal(*px,5); // 5==> Phi in the PtScal function
886b6f73 447}
448
75e0cc59 449 Double_t AliGenPHOSlib::YPhi( const Double_t *py, const Double_t *)
886b6f73 450{
f87cfe57 451// y-distribution
452//____________________________________________________________
453
454 const Double_t ka = 1000.;
455 const Double_t kdy = 4.;
886b6f73 456
457
458 Double_t y=TMath::Abs(*py);
459 //
f87cfe57 460 Double_t ex = y*y/(2*kdy*kdy);
461 return ka*TMath::Exp(-ex);
886b6f73 462}
463
65fb704d 464 Int_t AliGenPHOSlib::IpPhi(TRandom *)
f87cfe57 465{
886b6f73 466// particle composition
467//
f87cfe57 468
886b6f73 469 return 333; // Phi
470}
471// End Phis
472//===================================================================
473//============================================================================
474// B A R Y O N S == protons, protonsbar, neutrons, and neutronsbars
75e0cc59 475 Double_t AliGenPHOSlib::PtBaryon( const Double_t *px, const Double_t *)
f87cfe57 476{
886b6f73 477// baryons
478// pt-distribution
479//____________________________________________________________
f87cfe57 480
0db2f441 481 return PtScal(*px,6); // 6==> Baryon in the PtScal function
886b6f73 482}
483
75e0cc59 484 Double_t AliGenPHOSlib::YBaryon( const Double_t *py, const Double_t *)
886b6f73 485{
f87cfe57 486// y-distribution
487//____________________________________________________________
488
489 const Double_t ka = 1000.;
490 const Double_t kdy = 4.;
886b6f73 491
492
493 Double_t y=TMath::Abs(*py);
494 //
f87cfe57 495 Double_t ex = y*y/(2*kdy*kdy);
496 return ka*TMath::Exp(-ex);
886b6f73 497}
498
65fb704d 499 Int_t AliGenPHOSlib::IpBaryon(TRandom *ran)
886b6f73 500{
f87cfe57 501// particle composition
502//
503
65fb704d 504 Float_t random = ran->Rndm();
505 Float_t random2 = ran->Rndm();
506 if (random2 < 0.5)
886b6f73 507 {
65fb704d 508 if (random < 0.5) {
886b6f73 509 return 2212; // p
510 } else {
511 return -2212; // pbar
512 }
513 }
514 else
515 {
65fb704d 516 if (random < 0.5) {
886b6f73 517 return 2112; // n
518 } else {
519 return -2112; // n bar
520 }
521 }
522}
35ead6c1 523
524 Int_t AliGenPHOSlib::IpProton(TRandom *)
525{
526// particle composition
527//
528 return 2212; // p
529
530}
531 Int_t AliGenPHOSlib::IpAProton(TRandom *)
532{
533// particle composition
534//
535 return -2212; // p bar
536
537}
538
539 Int_t AliGenPHOSlib::IpNeutron(TRandom *)
540{
541// particle composition
542//
543 return 2112; // n
544
545}
546 Int_t AliGenPHOSlib::IpANeutron(TRandom *)
547{
548// particle composition
549//
550 return -2112; // n
551
552}
886b6f73 553// End Baryons
554//===================================================================
555
75e0cc59 556typedef Double_t (*GenFunc) (const Double_t*, const Double_t*);
198bb1c7 557GenFunc AliGenPHOSlib::GetPt(Int_t param, const char* /*tname*/) const
886b6f73 558{
f87cfe57 559// Return pinter to pT parameterisation
886b6f73 560 GenFunc func;
561
562 switch (param)
35ead6c1 563 {
564 case kPion:
886b6f73 565 func=PtPion;
566 break;
6784ce25 567 case kPi0:
568 func=PtPi0;
569 break;
35ead6c1 570 case kPi0Flat:
74c62c73 571 func=PtPi0Flat;
572 break;
35ead6c1 573 case kKaon:
886b6f73 574 func=PtKaon;
575 break;
35ead6c1 576 case kEta:
886b6f73 577 func=PtEta;
578 break;
35ead6c1 579 case kEtaFlat:
74c62c73 580 func=PtEtaFlat;
581 break;
35ead6c1 582 case kOmega:
886b6f73 583 func=PtOmega;
584 break;
8fc6910b 585 case kOmegaFlat:
586 func=PtOmegaFlat;
587 break;
35ead6c1 588 case kEtaPrime:
886b6f73 589 func=PtEtaprime;
590 break;
35ead6c1 591 case kBaryon:
886b6f73 592 func=PtBaryon;
593 break;
35ead6c1 594 default:
886b6f73 595 func=0;
596 printf("<AliGenPHOSlib::GetPt> unknown parametrisationn");
35ead6c1 597 }
886b6f73 598 return func;
599}
600
198bb1c7 601GenFunc AliGenPHOSlib::GetY(Int_t param, const char* /*tname*/) const
886b6f73 602{
35ead6c1 603 // Return pointer to Y parameterisation
604 GenFunc func;
605 switch (param)
886b6f73 606 {
34f60c01 607 case kPion:
35ead6c1 608 func=YPion;
609 break;
6784ce25 610 case kPi0:
74c62c73 611 case kPi0Flat:
35ead6c1 612 func=YPi0Flat;
613 break;
34f60c01 614 case kKaon:
35ead6c1 615 func=YKaon;
616 break;
34f60c01 617 case kEta:
35ead6c1 618 func=YEta;
619 break;
620 case kEtaFlat:
621 func=YEtaFlat;
622 break;
34f60c01 623 case kOmega:
35ead6c1 624 func=YOmega;
625 break;
8fc6910b 626 case kOmegaFlat:
627 func=YOmegaFlat;
628 break;
34f60c01 629 case kEtaPrime:
35ead6c1 630 func=YEtaprime;
631 break;
34f60c01 632 case kPhi:
35ead6c1 633 func=YPhi;
634 break;
34f60c01 635 case kBaryon:
35ead6c1 636 func=YBaryon;
637 break;
886b6f73 638 default:
35ead6c1 639 func=0;
640 printf("<AliGenPHOSlib::GetY> unknown parametrisationn");
886b6f73 641 }
35ead6c1 642 return func;
886b6f73 643}
65fb704d 644typedef Int_t (*GenFuncIp) (TRandom *);
198bb1c7 645GenFuncIp AliGenPHOSlib::GetIp(Int_t param, const char* /*tname*/) const
886b6f73 646{
35ead6c1 647 // Return pointer to particle composition
648 GenFuncIp func;
649 switch (param)
886b6f73 650 {
74c62c73 651 case kPion:
35ead6c1 652 func=IpPion;
653 break;
654 case kChargedPion:
655 func=IpChargedPion;
656 break;
6784ce25 657 case kPi0:
74c62c73 658 case kPi0Flat:
35ead6c1 659 func=IpPi0Flat;
660 break;
34f60c01 661 case kKaon:
35ead6c1 662 func=IpKaon;
663 break;
664 case kChargedKaon:
665 func=IpChargedKaon;
666 break;
667 case kKaon0L:
668 func=IpKaon0L;
669 break;
34f60c01 670 case kEta:
35ead6c1 671 func=IpEta;
672 break;
74c62c73 673 case kEtaFlat:
35ead6c1 674 func=IpEtaFlat;
675 break;
676
34f60c01 677 case kOmega:
35ead6c1 678 func=IpOmega;
679 break;
8fc6910b 680 case kOmegaFlat:
681 func=IpOmegaFlat;
682 break;
34f60c01 683 case kEtaPrime:
35ead6c1 684 func=IpEtaprime;
685 break;
34f60c01 686 case kPhi:
35ead6c1 687 func=IpPhi;
688 break;
34f60c01 689 case kBaryon:
35ead6c1 690 func=IpBaryon;
691 break;
692 case kProton:
693 func=IpProton;
694 break;
695 case kAProton:
696 func=IpAProton;
697 break;
698 case kNeutron:
699 func=IpNeutron;
700 break;
701 case kANeutron:
702 func=IpANeutron;
703 break;
704
886b6f73 705 default:
35ead6c1 706 func=0;
707 printf("<AliGenPHOSlib::GetIp> unknown parametrisationn");
886b6f73 708 }
35ead6c1 709 return func;
886b6f73 710}
711