Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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::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
140 Double_t AliGenPHOSlib::PtPi0Flat(const Double_t */*px*/, const Double_t *)
141 {
142 //     Pion transverse momentum flat distribution 
143
144 return 1;
145
146 }
147
148 Double_t AliGenPHOSlib::YPi0Flat( const Double_t */*py*/, const Double_t *)
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
164 //============================================================= 
165 //
166  Double_t AliGenPHOSlib::PtScal(Double_t pt, Int_t np)
167 {
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
175
176   //    SCALING EN MASSE PAR RAPPORT A PTPI
177   //    MASS                0=>PI,  1=>K, 2=>ETA, 3=>OMEGA,  4=>ETA',5=>PHI
178   const Double_t khm[10] = {0.1396, 0.494,  0.547,    0.782,   0.957,   1.02, 
179   //    MASS               6=>BARYON-BARYONBAR  
180                                          0.938, 0. , 0., 0.};
181   //     VALUE MESON/PI AT 5 GEV
182   const Double_t kfmax[10]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
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];
185   // PIONS
186   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
187   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
188                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
189   return fmtscal*ptpion;
190
191 }
192 // End Scaling
193 //============================================================================
194 //    K  A  O  N  S
195  Double_t AliGenPHOSlib::PtKaon( const Double_t *px, const Double_t *)
196 {
197 //                kaon
198 //                pt-distribution
199 //____________________________________________________________
200
201   return PtScal(*px,1);  //  1==> Kaon in the PtScal function
202 }
203
204  Double_t AliGenPHOSlib::YKaon( const Double_t *py, const Double_t *)
205 {
206 // y-distribution
207 //____________________________________________________________
208
209   const Double_t ka    = 1000.;
210   const Double_t kdy   = 4.;
211
212
213   Double_t y=TMath::Abs(*py);
214   //
215   Double_t ex = y*y/(2*kdy*kdy);
216   return ka*TMath::Exp(-ex);
217 }
218
219  Int_t AliGenPHOSlib::IpKaon(TRandom *ran)
220 {
221 //                 particle composition
222 //
223
224     Float_t random = ran->Rndm();
225     Float_t random2 = ran->Rndm();
226     if (random2 < 0.5) 
227     {
228       if (random < 0.5) {       
229         return  321;   //   K+
230       } else {
231         return -321;   // K-
232       }
233     }
234     else
235     {  
236       if (random < 0.5) {       
237         return  130;   // K^0 short
238       } else {  
239         return  310;   // K^0 long
240       }
241     }
242 }
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 }
259 Int_t AliGenPHOSlib::IpKaon0L(TRandom *)
260 {
261   //                 particle composition
262   //
263   
264         return  130;   // K^0 long
265 }
266 // End Kaons
267 //============================================================================
268 //============================================================================
269 //   E  T  A  S
270  Double_t AliGenPHOSlib::PtEta( const Double_t *px, const Double_t *)
271 {
272 //                etas
273 //                pt-distribution
274 //____________________________________________________________
275
276   return PtScal(*px,2);  //  2==> Eta in the PtScal function
277 }
278
279  Double_t AliGenPHOSlib::YEta( const Double_t *py, const Double_t *)
280 {
281 // y-distribution
282 //____________________________________________________________
283
284   const Double_t ka    = 1000.;
285   const Double_t kdy   = 4.;
286
287
288   Double_t y=TMath::Abs(*py);
289   //
290   Double_t ex = y*y/(2*kdy*kdy);
291   return ka*TMath::Exp(-ex);
292 }
293
294  Int_t AliGenPHOSlib::IpEta(TRandom *)
295 {
296 //                 particle composition
297 //
298
299         return  221;   //   eta
300 }
301 // End Etas
302
303 //======================================================================
304 //    Eta Flat Distribution
305 //    Transverse momentum distribution PtEtaFlat
306 //    Rapidity distribution YEtaFlat
307 //    Particle distribution IdEtaFlat  111 (pi0)
308 //
309
310 Double_t AliGenPHOSlib::PtEtaFlat(const Double_t */*px*/, const Double_t *)
311 {
312 //     Eta transverse momentum flat distribution 
313
314   return 1;
315
316 }
317
318 Double_t AliGenPHOSlib::YEtaFlat( const Double_t */*py*/, const Double_t *)
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 }
333 // End EtaFlat
334 //============================================================================
335 //============================================================================
336 //    O  M  E  G  A  S
337  Double_t AliGenPHOSlib::PtOmega( const Double_t *px, const Double_t *)
338 {
339 // omegas
340 //                pt-distribution
341 //____________________________________________________________
342
343   return PtScal(*px,3);  //  3==> Omega in the PtScal function
344 }
345
346  Double_t AliGenPHOSlib::YOmega( const Double_t *py, const Double_t *)
347 {
348 // y-distribution
349 //____________________________________________________________
350
351   const Double_t ka    = 1000.;
352   const Double_t kdy   = 4.;
353
354
355   Double_t y=TMath::Abs(*py);
356   //
357   Double_t ex = y*y/(2*kdy*kdy);
358   return ka*TMath::Exp(-ex);
359 }
360
361  Int_t AliGenPHOSlib::IpOmega(TRandom *)
362 {
363 //                 particle composition
364 //
365
366         return  223;   // Omega
367 }
368 // End Omega
369 //============================================================================
370 //======================================================================
371 //    Omega(782) Flat Distribution
372 //    Transverse momentum distribution PtOmegaFlat
373 //    Rapidity distribution YOmegaFlat
374 //    Particle distribution IdOmegaFlat  223(0mega)
375 //
376
377 Double_t AliGenPHOSlib::PtOmegaFlat(const Double_t */*px*/, const Double_t *)
378 {
379 //     omega transverse momentum flat distribution 
380
381 return 1;
382
383 }
384
385 Double_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
403 //============================================================================
404 //    E  T  A  P  R  I  M  E
405  Double_t AliGenPHOSlib::PtEtaprime( const Double_t *px, const Double_t *)
406 {
407 // etaprime
408 //                pt-distribution
409 //____________________________________________________________
410
411   return PtScal(*px,4);  //  4==> Etaprime in the PtScal function
412 }
413
414  Double_t AliGenPHOSlib::YEtaprime( const Double_t *py, const Double_t *)
415 {
416 // y-distribution
417 //____________________________________________________________
418
419   const Double_t ka    = 1000.;
420   const Double_t kdy   = 4.;
421
422
423   Double_t y=TMath::Abs(*py);
424   //
425   Double_t ex = y*y/(2*kdy*kdy);
426   return ka*TMath::Exp(-ex);
427 }
428
429  Int_t AliGenPHOSlib::IpEtaprime(TRandom *)
430 {
431 //                 particle composition
432 //
433
434         return  331;   //   Etaprime
435 }
436 // End EtaPrime
437 //===================================================================
438 //============================================================================
439 //    P  H  I   S
440  Double_t AliGenPHOSlib::PtPhi( const Double_t *px, const Double_t *)
441 {
442 // phi
443 //                pt-distribution
444 //____________________________________________________________
445
446   return PtScal(*px,5);  //  5==> Phi in the PtScal function
447 }
448
449  Double_t AliGenPHOSlib::YPhi( const Double_t *py, const Double_t *)
450 {
451 // y-distribution
452 //____________________________________________________________
453
454   const Double_t ka    = 1000.;
455   const Double_t kdy   = 4.;
456
457
458   Double_t y=TMath::Abs(*py);
459   //
460   Double_t ex = y*y/(2*kdy*kdy);
461   return ka*TMath::Exp(-ex);
462 }
463
464  Int_t AliGenPHOSlib::IpPhi(TRandom *)
465 {
466 //                 particle composition
467 //
468     
469         return  333;   //   Phi      
470 }
471 // End Phis
472 //===================================================================
473 //============================================================================
474 //    B  A  R  Y  O  N  S  == protons, protonsbar, neutrons, and neutronsbars
475  Double_t AliGenPHOSlib::PtBaryon( const Double_t *px, const Double_t *)
476 {
477 // baryons
478 //                pt-distribution
479 //____________________________________________________________
480
481   return PtScal(*px,6);  //  6==> Baryon in the PtScal function
482 }
483
484  Double_t AliGenPHOSlib::YBaryon( const Double_t *py, const Double_t *)
485 {
486 // y-distribution
487 //____________________________________________________________
488
489   const Double_t ka    = 1000.;
490   const Double_t kdy   = 4.;
491
492
493   Double_t y=TMath::Abs(*py);
494   //
495   Double_t ex = y*y/(2*kdy*kdy);
496   return ka*TMath::Exp(-ex);
497 }
498
499  Int_t AliGenPHOSlib::IpBaryon(TRandom *ran)
500 {
501 //                 particle composition
502 //
503
504     Float_t random = ran->Rndm();
505     Float_t random2 = ran->Rndm();
506     if (random2 < 0.5) 
507     {
508       if (random < 0.5) {       
509         return  2212;   //   p
510       } else {
511         return -2212;   // pbar
512       }
513     }
514     else
515     {  
516       if (random < 0.5) {       
517         return  2112;   // n
518       } else {  
519         return -2112;   // n bar
520       }
521     }
522 }
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 }
553 // End Baryons
554 //===================================================================
555
556 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
557 GenFunc AliGenPHOSlib::GetPt(Int_t param, const char* /*tname*/) const
558 {
559 // Return pinter to pT parameterisation
560     GenFunc func;
561     
562     switch (param)
563       {
564       case kPion:     
565         func=PtPion;
566         break;
567       case kPi0:     
568         func=PtPi0;
569         break;
570       case kPi0Flat:     
571         func=PtPi0Flat;
572         break;
573       case kKaon:
574         func=PtKaon;
575         break;
576       case kEta:
577         func=PtEta;
578         break;
579       case kEtaFlat:
580         func=PtEtaFlat;
581         break;
582       case kOmega:
583         func=PtOmega;
584         break;
585       case kOmegaFlat:
586         func=PtOmegaFlat;
587       break;
588       case kEtaPrime:
589         func=PtEtaprime;
590         break;
591       case kBaryon:
592         func=PtBaryon;
593         break;
594       default:
595         func=0;
596         printf("<AliGenPHOSlib::GetPt> unknown parametrisationn");
597       }
598     return func;
599 }
600
601 GenFunc AliGenPHOSlib::GetY(Int_t param, const char* /*tname*/) const
602 {
603   // Return pointer to Y parameterisation
604   GenFunc func;
605   switch (param)
606     {
607     case kPion:
608       func=YPion;
609       break;
610     case kPi0:     
611     case kPi0Flat:
612       func=YPi0Flat;
613       break;
614     case kKaon:
615       func=YKaon;
616       break;
617     case kEta:
618       func=YEta;
619       break;
620     case kEtaFlat:
621       func=YEtaFlat;
622       break;
623     case kOmega:
624       func=YOmega;
625       break;
626     case kOmegaFlat:
627       func=YOmegaFlat;
628       break;
629     case kEtaPrime:
630       func=YEtaprime;
631       break;
632     case kPhi:
633       func=YPhi;
634       break;
635     case kBaryon:
636       func=YBaryon;
637       break;
638     default:
639       func=0;
640       printf("<AliGenPHOSlib::GetY> unknown parametrisationn");
641     }
642   return func;
643 }
644 typedef Int_t (*GenFuncIp) (TRandom *);
645 GenFuncIp AliGenPHOSlib::GetIp(Int_t param,  const char* /*tname*/) const
646 {
647   // Return pointer to particle composition
648   GenFuncIp func;
649   switch (param)
650     {
651     case kPion:       
652       func=IpPion;
653       break;
654     case kChargedPion:       
655       func=IpChargedPion;
656       break;
657     case kPi0:     
658     case kPi0Flat:       
659       func=IpPi0Flat;
660       break;
661     case kKaon:
662       func=IpKaon;
663       break;
664     case kChargedKaon:
665       func=IpChargedKaon;
666       break;
667     case kKaon0L:
668       func=IpKaon0L;
669       break;
670     case kEta:
671       func=IpEta;
672       break;
673     case kEtaFlat:
674       func=IpEtaFlat;
675       break;
676       
677     case kOmega:
678       func=IpOmega;
679       break;
680     case kOmegaFlat:
681       func=IpOmegaFlat;
682       break;
683     case kEtaPrime:
684       func=IpEtaprime;
685       break;
686     case kPhi:
687       func=IpPhi;
688       break;
689     case kBaryon:
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       
705     default:
706       func=0;
707       printf("<AliGenPHOSlib::GetIp> unknown parametrisationn");
708     }
709   return func;
710 }
711