Fix for event mixing, when it was selecting events out of range of multiplicity cut
[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
75e0cc59 132Double_t AliGenPHOSlib::PtPi0Flat(const Double_t */*px*/, const Double_t *)
74c62c73 133{
134// Pion transverse momentum flat distribution
135
136return 1;
137
138}
139
75e0cc59 140Double_t AliGenPHOSlib::YPi0Flat( const Double_t */*py*/, const Double_t *)
74c62c73 141{
142
143// pion y-distribution
144//
145 return 1.;
146}
147
148 Int_t AliGenPHOSlib::IpPi0Flat(TRandom *)
149{
150
151// particle composition pi0
152//
153 return 111 ;
154}
155// End Pi0Flat
886b6f73 156//=============================================================
157//
f87cfe57 158 Double_t AliGenPHOSlib::PtScal(Double_t pt, Int_t np)
159{
886b6f73 160// Mt-scaling
161// Fonction for the calculation of the Pt distribution for a
162// given particle np, from the pion Pt distribution using the
163// mt scaling. This function was taken from AliGenMUONlib
164// aliroot version 3.01, and was extended for baryons
165// np = 1=>Pions 2=>Kaons 3=>Etas 4=>Omegas 5=>ETA' 6=>PHI
166// 7=>BARYONS-BARYONBARS
f87cfe57 167
886b6f73 168 // SCALING EN MASSE PAR RAPPORT A PTPI
0db2f441 169 // MASS 0=>PI, 1=>K, 2=>ETA, 3=>OMEGA, 4=>ETA',5=>PHI
f87cfe57 170 const Double_t khm[10] = {0.1396, 0.494, 0.547, 0.782, 0.957, 1.02,
0db2f441 171 // MASS 6=>BARYON-BARYONBAR
886b6f73 172 0.938, 0. , 0., 0.};
173 // VALUE MESON/PI AT 5 GEV
f87cfe57 174 const Double_t kfmax[10]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
f87cfe57 175 Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
176 Double_t kfmax2=f5/kfmax[np];
886b6f73 177 // PIONS
178 Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
179 Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
f87cfe57 180 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
886b6f73 181 return fmtscal*ptpion;
74c62c73 182
886b6f73 183}
184// End Scaling
185//============================================================================
186// K A O N S
75e0cc59 187 Double_t AliGenPHOSlib::PtKaon( const Double_t *px, const Double_t *)
886b6f73 188{
f87cfe57 189// kaon
190// pt-distribution
191//____________________________________________________________
192
0db2f441 193 return PtScal(*px,1); // 1==> Kaon in the PtScal function
886b6f73 194}
195
75e0cc59 196 Double_t AliGenPHOSlib::YKaon( const Double_t *py, const Double_t *)
886b6f73 197{
f87cfe57 198// y-distribution
199//____________________________________________________________
200
201 const Double_t ka = 1000.;
202 const Double_t kdy = 4.;
886b6f73 203
204
205 Double_t y=TMath::Abs(*py);
206 //
f87cfe57 207 Double_t ex = y*y/(2*kdy*kdy);
208 return ka*TMath::Exp(-ex);
886b6f73 209}
210
65fb704d 211 Int_t AliGenPHOSlib::IpKaon(TRandom *ran)
886b6f73 212{
f87cfe57 213// particle composition
214//
215
65fb704d 216 Float_t random = ran->Rndm();
217 Float_t random2 = ran->Rndm();
218 if (random2 < 0.5)
886b6f73 219 {
65fb704d 220 if (random < 0.5) {
886b6f73 221 return 321; // K+
222 } else {
223 return -321; // K-
224 }
225 }
226 else
227 {
65fb704d 228 if (random < 0.5) {
9ae59e17 229 return 130; // K^0 short
886b6f73 230 } else {
9ae59e17 231 return 310; // K^0 long
886b6f73 232 }
233 }
234}
35ead6c1 235
236 Int_t AliGenPHOSlib::IpChargedKaon(TRandom *ran)
237{
238// particle composition
239//
240
241 Float_t random = ran->Rndm();
242
243 if (random < 0.5) {
244 return 321; // K+
245 } else {
246 return -321; // K-
247 }
248
249
250}
251Int_t AliGenPHOSlib::IpKaon0L(TRandom *)
252{
253 // particle composition
254 //
255
256 return 130; // K^0 long
257}
886b6f73 258// End Kaons
259//============================================================================
260//============================================================================
261// E T A S
75e0cc59 262 Double_t AliGenPHOSlib::PtEta( const Double_t *px, const Double_t *)
886b6f73 263{
f87cfe57 264// etas
265// pt-distribution
266//____________________________________________________________
267
0db2f441 268 return PtScal(*px,2); // 2==> Eta in the PtScal function
886b6f73 269}
270
75e0cc59 271 Double_t AliGenPHOSlib::YEta( const Double_t *py, const Double_t *)
886b6f73 272{
f87cfe57 273// y-distribution
274//____________________________________________________________
275
276 const Double_t ka = 1000.;
277 const Double_t kdy = 4.;
886b6f73 278
279
280 Double_t y=TMath::Abs(*py);
281 //
f87cfe57 282 Double_t ex = y*y/(2*kdy*kdy);
283 return ka*TMath::Exp(-ex);
886b6f73 284}
285
65fb704d 286 Int_t AliGenPHOSlib::IpEta(TRandom *)
886b6f73 287{
f87cfe57 288// particle composition
289//
290
886b6f73 291 return 221; // eta
292}
293// End Etas
74c62c73 294
295//======================================================================
296// Eta Flat Distribution
297// Transverse momentum distribution PtEtaFlat
298// Rapidity distribution YEtaFlat
299// Particle distribution IdEtaFlat 111 (pi0)
300//
301
75e0cc59 302Double_t AliGenPHOSlib::PtEtaFlat(const Double_t */*px*/, const Double_t *)
74c62c73 303{
304// Eta transverse momentum flat distribution
305
306 return 1;
307
308}
309
75e0cc59 310Double_t AliGenPHOSlib::YEtaFlat( const Double_t */*py*/, const Double_t *)
74c62c73 311{
312//
313// pion y-distribution
314//
315 return 1.;
316}
317
318 Int_t AliGenPHOSlib::IpEtaFlat(TRandom *)
319{
320//
321// particle composition eta
322//
323 return 221 ;
324}
35ead6c1 325// End EtaFlat
886b6f73 326//============================================================================
327//============================================================================
328// O M E G A S
75e0cc59 329 Double_t AliGenPHOSlib::PtOmega( const Double_t *px, const Double_t *)
f87cfe57 330{
886b6f73 331// omegas
332// pt-distribution
333//____________________________________________________________
f87cfe57 334
0db2f441 335 return PtScal(*px,3); // 3==> Omega in the PtScal function
886b6f73 336}
337
75e0cc59 338 Double_t AliGenPHOSlib::YOmega( const Double_t *py, const Double_t *)
886b6f73 339{
f87cfe57 340// y-distribution
341//____________________________________________________________
342
343 const Double_t ka = 1000.;
344 const Double_t kdy = 4.;
886b6f73 345
346
347 Double_t y=TMath::Abs(*py);
348 //
f87cfe57 349 Double_t ex = y*y/(2*kdy*kdy);
350 return ka*TMath::Exp(-ex);
886b6f73 351}
352
65fb704d 353 Int_t AliGenPHOSlib::IpOmega(TRandom *)
886b6f73 354{
f87cfe57 355// particle composition
356//
357
886b6f73 358 return 223; // Omega
359}
360// End Omega
361//============================================================================
8fc6910b 362//======================================================================
363// Omega(782) Flat Distribution
364// Transverse momentum distribution PtOmegaFlat
365// Rapidity distribution YOmegaFlat
366// Particle distribution IdOmegaFlat 223(0mega)
367//
368
369Double_t AliGenPHOSlib::PtOmegaFlat(const Double_t */*px*/, const Double_t *)
370{
371// omega transverse momentum flat distribution
372
373return 1;
374
375}
376
377Double_t AliGenPHOSlib::YOmegaFlat( const Double_t */*py*/, const Double_t *)
378{
379
380// omega y-distribution
381//
382 return 1.;
383}
384
385 Int_t AliGenPHOSlib::IpOmegaFlat(TRandom *)
386{
387
388// particle composition omega
389//
390 return 223 ;
391}
392// End OmegaFlat
393
394
886b6f73 395//============================================================================
396// E T A P R I M E
75e0cc59 397 Double_t AliGenPHOSlib::PtEtaprime( const Double_t *px, const Double_t *)
f87cfe57 398{
886b6f73 399// etaprime
400// pt-distribution
401//____________________________________________________________
f87cfe57 402
0db2f441 403 return PtScal(*px,4); // 4==> Etaprime in the PtScal function
886b6f73 404}
405
75e0cc59 406 Double_t AliGenPHOSlib::YEtaprime( const Double_t *py, const Double_t *)
886b6f73 407{
f87cfe57 408// y-distribution
409//____________________________________________________________
410
411 const Double_t ka = 1000.;
412 const Double_t kdy = 4.;
886b6f73 413
414
415 Double_t y=TMath::Abs(*py);
416 //
f87cfe57 417 Double_t ex = y*y/(2*kdy*kdy);
418 return ka*TMath::Exp(-ex);
886b6f73 419}
420
65fb704d 421 Int_t AliGenPHOSlib::IpEtaprime(TRandom *)
886b6f73 422{
f87cfe57 423// particle composition
424//
425
886b6f73 426 return 331; // Etaprime
427}
428// End EtaPrime
429//===================================================================
430//============================================================================
431// P H I S
75e0cc59 432 Double_t AliGenPHOSlib::PtPhi( const Double_t *px, const Double_t *)
f87cfe57 433{
886b6f73 434// phi
435// pt-distribution
436//____________________________________________________________
f87cfe57 437
0db2f441 438 return PtScal(*px,5); // 5==> Phi in the PtScal function
886b6f73 439}
440
75e0cc59 441 Double_t AliGenPHOSlib::YPhi( const Double_t *py, const Double_t *)
886b6f73 442{
f87cfe57 443// y-distribution
444//____________________________________________________________
445
446 const Double_t ka = 1000.;
447 const Double_t kdy = 4.;
886b6f73 448
449
450 Double_t y=TMath::Abs(*py);
451 //
f87cfe57 452 Double_t ex = y*y/(2*kdy*kdy);
453 return ka*TMath::Exp(-ex);
886b6f73 454}
455
65fb704d 456 Int_t AliGenPHOSlib::IpPhi(TRandom *)
f87cfe57 457{
886b6f73 458// particle composition
459//
f87cfe57 460
886b6f73 461 return 333; // Phi
462}
463// End Phis
464//===================================================================
465//============================================================================
466// B A R Y O N S == protons, protonsbar, neutrons, and neutronsbars
75e0cc59 467 Double_t AliGenPHOSlib::PtBaryon( const Double_t *px, const Double_t *)
f87cfe57 468{
886b6f73 469// baryons
470// pt-distribution
471//____________________________________________________________
f87cfe57 472
0db2f441 473 return PtScal(*px,6); // 6==> Baryon in the PtScal function
886b6f73 474}
475
75e0cc59 476 Double_t AliGenPHOSlib::YBaryon( const Double_t *py, const Double_t *)
886b6f73 477{
f87cfe57 478// y-distribution
479//____________________________________________________________
480
481 const Double_t ka = 1000.;
482 const Double_t kdy = 4.;
886b6f73 483
484
485 Double_t y=TMath::Abs(*py);
486 //
f87cfe57 487 Double_t ex = y*y/(2*kdy*kdy);
488 return ka*TMath::Exp(-ex);
886b6f73 489}
490
65fb704d 491 Int_t AliGenPHOSlib::IpBaryon(TRandom *ran)
886b6f73 492{
f87cfe57 493// particle composition
494//
495
65fb704d 496 Float_t random = ran->Rndm();
497 Float_t random2 = ran->Rndm();
498 if (random2 < 0.5)
886b6f73 499 {
65fb704d 500 if (random < 0.5) {
886b6f73 501 return 2212; // p
502 } else {
503 return -2212; // pbar
504 }
505 }
506 else
507 {
65fb704d 508 if (random < 0.5) {
886b6f73 509 return 2112; // n
510 } else {
511 return -2112; // n bar
512 }
513 }
514}
35ead6c1 515
516 Int_t AliGenPHOSlib::IpProton(TRandom *)
517{
518// particle composition
519//
520 return 2212; // p
521
522}
523 Int_t AliGenPHOSlib::IpAProton(TRandom *)
524{
525// particle composition
526//
527 return -2212; // p bar
528
529}
530
531 Int_t AliGenPHOSlib::IpNeutron(TRandom *)
532{
533// particle composition
534//
535 return 2112; // n
536
537}
538 Int_t AliGenPHOSlib::IpANeutron(TRandom *)
539{
540// particle composition
541//
542 return -2112; // n
543
544}
886b6f73 545// End Baryons
546//===================================================================
547
75e0cc59 548typedef Double_t (*GenFunc) (const Double_t*, const Double_t*);
198bb1c7 549GenFunc AliGenPHOSlib::GetPt(Int_t param, const char* /*tname*/) const
886b6f73 550{
f87cfe57 551// Return pinter to pT parameterisation
886b6f73 552 GenFunc func;
553
554 switch (param)
35ead6c1 555 {
556 case kPion:
886b6f73 557 func=PtPion;
558 break;
35ead6c1 559 case kPi0Flat:
74c62c73 560 func=PtPi0Flat;
561 break;
35ead6c1 562 case kKaon:
886b6f73 563 func=PtKaon;
564 break;
35ead6c1 565 case kEta:
886b6f73 566 func=PtEta;
567 break;
35ead6c1 568 case kEtaFlat:
74c62c73 569 func=PtEtaFlat;
570 break;
35ead6c1 571 case kOmega:
886b6f73 572 func=PtOmega;
573 break;
8fc6910b 574 case kOmegaFlat:
575 func=PtOmegaFlat;
576 break;
35ead6c1 577 case kEtaPrime:
886b6f73 578 func=PtEtaprime;
579 break;
35ead6c1 580 case kBaryon:
886b6f73 581 func=PtBaryon;
582 break;
35ead6c1 583 default:
886b6f73 584 func=0;
585 printf("<AliGenPHOSlib::GetPt> unknown parametrisationn");
35ead6c1 586 }
886b6f73 587 return func;
588}
589
198bb1c7 590GenFunc AliGenPHOSlib::GetY(Int_t param, const char* /*tname*/) const
886b6f73 591{
35ead6c1 592 // Return pointer to Y parameterisation
593 GenFunc func;
594 switch (param)
886b6f73 595 {
34f60c01 596 case kPion:
35ead6c1 597 func=YPion;
598 break;
74c62c73 599 case kPi0Flat:
35ead6c1 600 func=YPi0Flat;
601 break;
34f60c01 602 case kKaon:
35ead6c1 603 func=YKaon;
604 break;
34f60c01 605 case kEta:
35ead6c1 606 func=YEta;
607 break;
608 case kEtaFlat:
609 func=YEtaFlat;
610 break;
34f60c01 611 case kOmega:
35ead6c1 612 func=YOmega;
613 break;
8fc6910b 614 case kOmegaFlat:
615 func=YOmegaFlat;
616 break;
34f60c01 617 case kEtaPrime:
35ead6c1 618 func=YEtaprime;
619 break;
34f60c01 620 case kPhi:
35ead6c1 621 func=YPhi;
622 break;
34f60c01 623 case kBaryon:
35ead6c1 624 func=YBaryon;
625 break;
886b6f73 626 default:
35ead6c1 627 func=0;
628 printf("<AliGenPHOSlib::GetY> unknown parametrisationn");
886b6f73 629 }
35ead6c1 630 return func;
886b6f73 631}
65fb704d 632typedef Int_t (*GenFuncIp) (TRandom *);
198bb1c7 633GenFuncIp AliGenPHOSlib::GetIp(Int_t param, const char* /*tname*/) const
886b6f73 634{
35ead6c1 635 // Return pointer to particle composition
636 GenFuncIp func;
637 switch (param)
886b6f73 638 {
74c62c73 639 case kPion:
35ead6c1 640 func=IpPion;
641 break;
642 case kChargedPion:
643 func=IpChargedPion;
644 break;
74c62c73 645 case kPi0Flat:
35ead6c1 646 func=IpPi0Flat;
647 break;
34f60c01 648 case kKaon:
35ead6c1 649 func=IpKaon;
650 break;
651 case kChargedKaon:
652 func=IpChargedKaon;
653 break;
654 case kKaon0L:
655 func=IpKaon0L;
656 break;
34f60c01 657 case kEta:
35ead6c1 658 func=IpEta;
659 break;
74c62c73 660 case kEtaFlat:
35ead6c1 661 func=IpEtaFlat;
662 break;
663
34f60c01 664 case kOmega:
35ead6c1 665 func=IpOmega;
666 break;
8fc6910b 667 case kOmegaFlat:
668 func=IpOmegaFlat;
669 break;
34f60c01 670 case kEtaPrime:
35ead6c1 671 func=IpEtaprime;
672 break;
34f60c01 673 case kPhi:
35ead6c1 674 func=IpPhi;
675 break;
34f60c01 676 case kBaryon:
35ead6c1 677 func=IpBaryon;
678 break;
679 case kProton:
680 func=IpProton;
681 break;
682 case kAProton:
683 func=IpAProton;
684 break;
685 case kNeutron:
686 func=IpNeutron;
687 break;
688 case kANeutron:
689 func=IpANeutron;
690 break;
691
886b6f73 692 default:
35ead6c1 693 func=0;
694 printf("<AliGenPHOSlib::GetIp> unknown parametrisationn");
886b6f73 695 }
35ead6c1 696 return func;
886b6f73 697}
698