Fix for event mixing, when it was selecting events out of range of multiplicity cut
[u/mrichter/AliRoot.git] / EVGEN / AliGenPHOSlib.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 //======================================================================
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 
31 //  with Pions, using the PtScal function taken from AliGenMUONlib 
32 //  version aliroot 3.01
33 //
34 //     Gines MARTINEZ. Laurent APHECETCHE and Yves SCHUTZ
35 //      GPS @ SUBATECH,  Nantes , France  (October 1999)
36 //     http://www-subatech.in2p3.fr/~photons/subatech
37 //     martinez@subatech.in2p3.fr
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
41 //  Add flat Omega(782) distribution in Nov. 2010 by Renzhuo WAN
42 //======================================================================
43
44 #include "TMath.h"
45 #include "TRandom.h"
46
47 #include "AliGenPHOSlib.h"
48
49 ClassImp(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 //
58  Double_t AliGenPHOSlib::PtPion(const Double_t *px, const Double_t *)
59 {
60 //     Pion transverse momentum distribtuion taken 
61 //     from AliGenMUONlib class, version 3.01 of aliroot
62 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
63 //     POWER LAW FOR PT > 500 MEV
64 //     MT SCALING BELOW (T=160 MEV)
65 //
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;
73   Double_t x=*px;
74   //
75   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
76   kxmpi2=kxmpi*kxmpi;
77   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+kxmpi2)/kt));
78   a=ynorm/y1;
79   if (x > kxlim)
80     y=a*TMath::Power(kp0/(kp0+x),kxn);
81   else
82     y=kb*TMath::Exp(-sqrt(x*x+kxmpi2)/kt);
83   return y*x;
84 }
85  Double_t AliGenPHOSlib::YPion( const Double_t *py, const Double_t *)
86 {
87 //
88 // pion y-distribution
89 //
90
91   const Double_t ka    = 7000.;   
92   const Double_t kdy   = 4.;
93
94   Double_t y=TMath::Abs(*py);
95   //
96   Double_t ex = y*y/(2*kdy*kdy);
97   return ka*TMath::Exp(-ex);
98 }
99
100 Int_t AliGenPHOSlib::IpPion(TRandom */*ran*/)
101 {
102 //                 particle composition pi+, pi0, pi-
103 //
104
105         return 111  ;
106 }
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 }
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
132 Double_t AliGenPHOSlib::PtPi0Flat(const Double_t */*px*/, const Double_t *)
133 {
134 //     Pion transverse momentum flat distribution 
135
136 return 1;
137
138 }
139
140 Double_t AliGenPHOSlib::YPi0Flat( const Double_t */*py*/, const Double_t *)
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
156 //============================================================= 
157 //
158  Double_t AliGenPHOSlib::PtScal(Double_t pt, Int_t np)
159 {
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
167
168   //    SCALING EN MASSE PAR RAPPORT A PTPI
169   //    MASS                0=>PI,  1=>K, 2=>ETA, 3=>OMEGA,  4=>ETA',5=>PHI
170   const Double_t khm[10] = {0.1396, 0.494,  0.547,    0.782,   0.957,   1.02, 
171   //    MASS               6=>BARYON-BARYONBAR  
172                                          0.938, 0. , 0., 0.};
173   //     VALUE MESON/PI AT 5 GEV
174   const Double_t kfmax[10]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
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];
177   // PIONS
178   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
179   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
180                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
181   return fmtscal*ptpion;
182
183 }
184 // End Scaling
185 //============================================================================
186 //    K  A  O  N  S
187  Double_t AliGenPHOSlib::PtKaon( const Double_t *px, const Double_t *)
188 {
189 //                kaon
190 //                pt-distribution
191 //____________________________________________________________
192
193   return PtScal(*px,1);  //  1==> Kaon in the PtScal function
194 }
195
196  Double_t AliGenPHOSlib::YKaon( const Double_t *py, const Double_t *)
197 {
198 // y-distribution
199 //____________________________________________________________
200
201   const Double_t ka    = 1000.;
202   const Double_t kdy   = 4.;
203
204
205   Double_t y=TMath::Abs(*py);
206   //
207   Double_t ex = y*y/(2*kdy*kdy);
208   return ka*TMath::Exp(-ex);
209 }
210
211  Int_t AliGenPHOSlib::IpKaon(TRandom *ran)
212 {
213 //                 particle composition
214 //
215
216     Float_t random = ran->Rndm();
217     Float_t random2 = ran->Rndm();
218     if (random2 < 0.5) 
219     {
220       if (random < 0.5) {       
221         return  321;   //   K+
222       } else {
223         return -321;   // K-
224       }
225     }
226     else
227     {  
228       if (random < 0.5) {       
229         return  130;   // K^0 short
230       } else {  
231         return  310;   // K^0 long
232       }
233     }
234 }
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 }
251 Int_t AliGenPHOSlib::IpKaon0L(TRandom *)
252 {
253   //                 particle composition
254   //
255   
256         return  130;   // K^0 long
257 }
258 // End Kaons
259 //============================================================================
260 //============================================================================
261 //   E  T  A  S
262  Double_t AliGenPHOSlib::PtEta( const Double_t *px, const Double_t *)
263 {
264 //                etas
265 //                pt-distribution
266 //____________________________________________________________
267
268   return PtScal(*px,2);  //  2==> Eta in the PtScal function
269 }
270
271  Double_t AliGenPHOSlib::YEta( const Double_t *py, const Double_t *)
272 {
273 // y-distribution
274 //____________________________________________________________
275
276   const Double_t ka    = 1000.;
277   const Double_t kdy   = 4.;
278
279
280   Double_t y=TMath::Abs(*py);
281   //
282   Double_t ex = y*y/(2*kdy*kdy);
283   return ka*TMath::Exp(-ex);
284 }
285
286  Int_t AliGenPHOSlib::IpEta(TRandom *)
287 {
288 //                 particle composition
289 //
290
291         return  221;   //   eta
292 }
293 // End Etas
294
295 //======================================================================
296 //    Eta Flat Distribution
297 //    Transverse momentum distribution PtEtaFlat
298 //    Rapidity distribution YEtaFlat
299 //    Particle distribution IdEtaFlat  111 (pi0)
300 //
301
302 Double_t AliGenPHOSlib::PtEtaFlat(const Double_t */*px*/, const Double_t *)
303 {
304 //     Eta transverse momentum flat distribution 
305
306   return 1;
307
308 }
309
310 Double_t AliGenPHOSlib::YEtaFlat( const Double_t */*py*/, const Double_t *)
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 }
325 // End EtaFlat
326 //============================================================================
327 //============================================================================
328 //    O  M  E  G  A  S
329  Double_t AliGenPHOSlib::PtOmega( const Double_t *px, const Double_t *)
330 {
331 // omegas
332 //                pt-distribution
333 //____________________________________________________________
334
335   return PtScal(*px,3);  //  3==> Omega in the PtScal function
336 }
337
338  Double_t AliGenPHOSlib::YOmega( const Double_t *py, const Double_t *)
339 {
340 // y-distribution
341 //____________________________________________________________
342
343   const Double_t ka    = 1000.;
344   const Double_t kdy   = 4.;
345
346
347   Double_t y=TMath::Abs(*py);
348   //
349   Double_t ex = y*y/(2*kdy*kdy);
350   return ka*TMath::Exp(-ex);
351 }
352
353  Int_t AliGenPHOSlib::IpOmega(TRandom *)
354 {
355 //                 particle composition
356 //
357
358         return  223;   // Omega
359 }
360 // End Omega
361 //============================================================================
362 //======================================================================
363 //    Omega(782) Flat Distribution
364 //    Transverse momentum distribution PtOmegaFlat
365 //    Rapidity distribution YOmegaFlat
366 //    Particle distribution IdOmegaFlat  223(0mega)
367 //
368
369 Double_t AliGenPHOSlib::PtOmegaFlat(const Double_t */*px*/, const Double_t *)
370 {
371 //     omega transverse momentum flat distribution 
372
373 return 1;
374
375 }
376
377 Double_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
395 //============================================================================
396 //    E  T  A  P  R  I  M  E
397  Double_t AliGenPHOSlib::PtEtaprime( const Double_t *px, const Double_t *)
398 {
399 // etaprime
400 //                pt-distribution
401 //____________________________________________________________
402
403   return PtScal(*px,4);  //  4==> Etaprime in the PtScal function
404 }
405
406  Double_t AliGenPHOSlib::YEtaprime( const Double_t *py, const Double_t *)
407 {
408 // y-distribution
409 //____________________________________________________________
410
411   const Double_t ka    = 1000.;
412   const Double_t kdy   = 4.;
413
414
415   Double_t y=TMath::Abs(*py);
416   //
417   Double_t ex = y*y/(2*kdy*kdy);
418   return ka*TMath::Exp(-ex);
419 }
420
421  Int_t AliGenPHOSlib::IpEtaprime(TRandom *)
422 {
423 //                 particle composition
424 //
425
426         return  331;   //   Etaprime
427 }
428 // End EtaPrime
429 //===================================================================
430 //============================================================================
431 //    P  H  I   S
432  Double_t AliGenPHOSlib::PtPhi( const Double_t *px, const Double_t *)
433 {
434 // phi
435 //                pt-distribution
436 //____________________________________________________________
437
438   return PtScal(*px,5);  //  5==> Phi in the PtScal function
439 }
440
441  Double_t AliGenPHOSlib::YPhi( const Double_t *py, const Double_t *)
442 {
443 // y-distribution
444 //____________________________________________________________
445
446   const Double_t ka    = 1000.;
447   const Double_t kdy   = 4.;
448
449
450   Double_t y=TMath::Abs(*py);
451   //
452   Double_t ex = y*y/(2*kdy*kdy);
453   return ka*TMath::Exp(-ex);
454 }
455
456  Int_t AliGenPHOSlib::IpPhi(TRandom *)
457 {
458 //                 particle composition
459 //
460     
461         return  333;   //   Phi      
462 }
463 // End Phis
464 //===================================================================
465 //============================================================================
466 //    B  A  R  Y  O  N  S  == protons, protonsbar, neutrons, and neutronsbars
467  Double_t AliGenPHOSlib::PtBaryon( const Double_t *px, const Double_t *)
468 {
469 // baryons
470 //                pt-distribution
471 //____________________________________________________________
472
473   return PtScal(*px,6);  //  6==> Baryon in the PtScal function
474 }
475
476  Double_t AliGenPHOSlib::YBaryon( const Double_t *py, const Double_t *)
477 {
478 // y-distribution
479 //____________________________________________________________
480
481   const Double_t ka    = 1000.;
482   const Double_t kdy   = 4.;
483
484
485   Double_t y=TMath::Abs(*py);
486   //
487   Double_t ex = y*y/(2*kdy*kdy);
488   return ka*TMath::Exp(-ex);
489 }
490
491  Int_t AliGenPHOSlib::IpBaryon(TRandom *ran)
492 {
493 //                 particle composition
494 //
495
496     Float_t random = ran->Rndm();
497     Float_t random2 = ran->Rndm();
498     if (random2 < 0.5) 
499     {
500       if (random < 0.5) {       
501         return  2212;   //   p
502       } else {
503         return -2212;   // pbar
504       }
505     }
506     else
507     {  
508       if (random < 0.5) {       
509         return  2112;   // n
510       } else {  
511         return -2112;   // n bar
512       }
513     }
514 }
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 }
545 // End Baryons
546 //===================================================================
547
548 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
549 GenFunc AliGenPHOSlib::GetPt(Int_t param, const char* /*tname*/) const
550 {
551 // Return pinter to pT parameterisation
552     GenFunc func;
553     
554     switch (param)
555       {
556       case kPion:     
557         func=PtPion;
558         break;
559       case kPi0Flat:     
560         func=PtPi0Flat;
561         break;
562       case kKaon:
563         func=PtKaon;
564         break;
565       case kEta:
566         func=PtEta;
567         break;
568       case kEtaFlat:
569         func=PtEtaFlat;
570         break;
571       case kOmega:
572         func=PtOmega;
573         break;
574       case kOmegaFlat:
575         func=PtOmegaFlat;
576       break;
577       case kEtaPrime:
578         func=PtEtaprime;
579         break;
580       case kBaryon:
581         func=PtBaryon;
582         break;
583       default:
584         func=0;
585         printf("<AliGenPHOSlib::GetPt> unknown parametrisationn");
586       }
587     return func;
588 }
589
590 GenFunc AliGenPHOSlib::GetY(Int_t param, const char* /*tname*/) const
591 {
592   // Return pointer to Y parameterisation
593   GenFunc func;
594   switch (param)
595     {
596     case kPion:
597       func=YPion;
598       break;
599     case kPi0Flat:
600       func=YPi0Flat;
601       break;
602     case kKaon:
603       func=YKaon;
604       break;
605     case kEta:
606       func=YEta;
607       break;
608     case kEtaFlat:
609       func=YEtaFlat;
610       break;
611     case kOmega:
612       func=YOmega;
613       break;
614     case kOmegaFlat:
615       func=YOmegaFlat;
616       break;
617     case kEtaPrime:
618       func=YEtaprime;
619       break;
620     case kPhi:
621       func=YPhi;
622       break;
623     case kBaryon:
624       func=YBaryon;
625       break;
626     default:
627       func=0;
628       printf("<AliGenPHOSlib::GetY> unknown parametrisationn");
629     }
630   return func;
631 }
632 typedef Int_t (*GenFuncIp) (TRandom *);
633 GenFuncIp AliGenPHOSlib::GetIp(Int_t param,  const char* /*tname*/) const
634 {
635   // Return pointer to particle composition
636   GenFuncIp func;
637   switch (param)
638     {
639     case kPion:       
640       func=IpPion;
641       break;
642     case kChargedPion:       
643       func=IpChargedPion;
644       break;
645     case kPi0Flat:       
646       func=IpPi0Flat;
647       break;
648     case kKaon:
649       func=IpKaon;
650       break;
651     case kChargedKaon:
652       func=IpChargedKaon;
653       break;
654     case kKaon0L:
655       func=IpKaon0L;
656       break;
657     case kEta:
658       func=IpEta;
659       break;
660     case kEtaFlat:
661       func=IpEtaFlat;
662       break;
663       
664     case kOmega:
665       func=IpOmega;
666       break;
667     case kOmegaFlat:
668       func=IpOmegaFlat;
669       break;
670     case kEtaPrime:
671       func=IpEtaprime;
672       break;
673     case kPhi:
674       func=IpPhi;
675       break;
676     case kBaryon:
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       
692     default:
693       func=0;
694       printf("<AliGenPHOSlib::GetIp> unknown parametrisationn");
695     }
696   return func;
697 }
698