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