Recover mods from Rev. 1.8
[u/mrichter/AliRoot.git] / EVGEN / AliGenPHOSlib.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /*
18 $Log$
19 Revision 1.9  2002/04/23 12:54:29  morsch
20 New options kPi0Flat y kEtaFlat (Gustavo Conesa)
21
22 Revision 1.7  2001/03/09 13:01:41  morsch
23 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
24 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
25
26 Revision 1.6  2000/11/30 07:12:50  alibrary
27 Introducing new Rndm and QA classes
28
29 Revision 1.5  2000/06/29 21:08:27  morsch
30 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
31 This allows to pass a pointer to a library directly to AliGenParam and avoids the
32 use of function pointers in Config.C.
33
34 Revision 1.4  2000/06/14 15:21:05  morsch
35 Include clean-up (IH)
36
37 Revision 1.3  2000/06/09 20:32:54  morsch
38 All coding rule violations except RS3 corrected
39
40 Revision 1.2  1999/11/04 11:30:48  fca
41 Improve comments
42
43 Revision 1.1  1999/11/03 17:43:20  fca
44 New version from G.Martinez & A.Morsch
45
46 */
47
48 //======================================================================
49 //  AliGenPHOSlib class contains parameterizations of the
50 //  pion, kaon, eta, omega, etaprime, phi and baryon (proton, 
51 //  antiproton, neutron and anti-neutron) particles for the 
52 //  study of the neutral background in PHOS detector. 
53 //  These parameterizations are used by the 
54 //  AliGenParam  class:
55 //  AliGenParam(npar, param,  AliGenPHOSlib::GetPt(param),
56 //                            AliGenPHOSlib::GetY(param),
57 //                            AliGenPHOSlib::GetIp(param) )
58 //  param represents the particle to be simulated : 
59 //  Pion, Kaon, Eta, Omega, Etaprime, Phi or Baryon    
60 //  Pt distributions are calculated from the transverse mass scaling 
61 //  with Pions, using the PtScal function taken from AliGenMUONlib 
62 //  version aliroot 3.01
63 //
64 //     Gines MARTINEZ. Laurent APHECETCHE and Yves SCHUTZ
65 //      GPS @ SUBATECH,  Nantes , France  (October 1999)
66 //     http://www-subatech.in2p3.fr/~photons/subatech
67 //     martinez@subatech.in2p3.fr
68 //======================================================================
69
70 #include "TMath.h"
71 #include "TRandom.h"
72
73 #include "AliGenPHOSlib.h"
74
75 ClassImp(AliGenPHOSlib)
76
77 //======================================================================
78 //    P  I  O  N  S
79 //    (From GetPt, GetY and GetIp as param = Pion)
80 //    Transverse momentum distribution" PtPion 
81 //    Rapidity distribution YPion
82 //    Particle distribution IdPion  111, 211 and -211 (pi0, pi+ and pi-)
83 //
84  Double_t AliGenPHOSlib::PtPion(Double_t *px, Double_t *)
85 {
86 //     Pion transverse momentum distribtuion taken 
87 //     from AliGenMUONlib class, version 3.01 of aliroot
88 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
89 //     POWER LAW FOR PT > 500 MEV
90 //     MT SCALING BELOW (T=160 MEV)
91 //
92   const Double_t kp0 = 1.3;
93   const Double_t kxn = 8.28;
94   const Double_t kxlim=0.5;
95   const Double_t kt=0.160;
96   const Double_t kxmpi=0.139;
97   const Double_t kb=1.;
98   Double_t y, y1, kxmpi2, ynorm, a;
99   Double_t x=*px;
100   //
101   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
102   kxmpi2=kxmpi*kxmpi;
103   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+kxmpi2)/kt));
104   a=ynorm/y1;
105   if (x > kxlim)
106     y=a*TMath::Power(kp0/(kp0+x),kxn);
107   else
108     y=kb*TMath::Exp(-sqrt(x*x+kxmpi2)/kt);
109   return y*x;
110 }
111  Double_t AliGenPHOSlib::YPion( Double_t *py, Double_t *)
112 {
113 //
114 // pion y-distribution
115 //
116
117   const Double_t ka    = 7000.;   
118   const Double_t kdy   = 4.;
119
120   Double_t y=TMath::Abs(*py);
121   //
122   Double_t ex = y*y/(2*kdy*kdy);
123   return ka*TMath::Exp(-ex);
124 }
125
126  Int_t AliGenPHOSlib::IpPion(TRandom *ran)
127 {
128 //                 particle composition pi+, pi0, pi-
129 //
130
131      Float_t random = ran->Rndm();
132
133      if ( (3.*random)  < 1. ) 
134      {
135            return 211 ;
136      } 
137      else
138      {  
139        if ( (3.*random) >= 2.)
140        {
141           return -211 ;
142        }
143        else 
144        {
145         return 111  ;
146       }
147     }
148 }
149
150 //End Pions
151 //======================================================================
152 //    Pi 0 Flat Distribution
153 //    Transverse momentum distribution PtPi0Flat
154 //    Rapidity distribution YPi0Flat
155 //    Particle distribution IdPi0Flat  111 (pi0)
156 //
157
158 Double_t AliGenPHOSlib::PtPi0Flat(Double_t *px, Double_t *)
159 {
160 //     Pion transverse momentum flat distribution 
161
162 return 1;
163
164 }
165
166 Double_t AliGenPHOSlib::YPi0Flat( Double_t *py, Double_t *)
167 {
168
169 // pion y-distribution
170 //
171   return 1.;
172 }
173
174  Int_t AliGenPHOSlib::IpPi0Flat(TRandom *)
175 {
176
177 //                 particle composition pi0
178 //
179         return 111 ;
180 }
181 // End Pi0Flat
182 //============================================================= 
183 //
184  Double_t AliGenPHOSlib::PtScal(Double_t pt, Int_t np)
185 {
186 // Mt-scaling
187 // Fonction for the calculation of the Pt distribution for a 
188 // given particle np, from the pion Pt distribution using the 
189 // mt scaling. This function was taken from AliGenMUONlib 
190 // aliroot version 3.01, and was extended for baryons
191 // np = 1=>Pions 2=>Kaons 3=>Etas 4=>Omegas 5=>ETA' 6=>PHI
192 //      7=>BARYONS-BARYONBARS
193
194   //    SCALING EN MASSE PAR RAPPORT A PTPI
195   //    MASS                1=>PI,  2=>K, 3=>ETA, 4=>OMEGA,  5=>ETA',6=>PHI
196   const Double_t khm[10] = {0.1396, 0.494,  0.547,    0.782,   0.957,   1.02, 
197   //    MASS               7=>BARYON-BARYONBAR  
198                                          0.938, 0. , 0., 0.};
199   //     VALUE MESON/PI AT 5 GEV
200   const Double_t kfmax[10]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
201   np--;
202   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
203   Double_t kfmax2=f5/kfmax[np];
204   // PIONS
205   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
206   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
207                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
208   return fmtscal*ptpion;
209
210 }
211 // End Scaling
212 //============================================================================
213 //    K  A  O  N  S
214  Double_t AliGenPHOSlib::PtKaon( Double_t *px, Double_t *)
215 {
216 //                kaon
217 //                pt-distribution
218 //____________________________________________________________
219
220   return PtScal(*px,2);  //  2==> Kaon in the PtScal function
221 }
222
223  Double_t AliGenPHOSlib::YKaon( Double_t *py, Double_t *)
224 {
225 // y-distribution
226 //____________________________________________________________
227
228   const Double_t ka    = 1000.;
229   const Double_t kdy   = 4.;
230
231
232   Double_t y=TMath::Abs(*py);
233   //
234   Double_t ex = y*y/(2*kdy*kdy);
235   return ka*TMath::Exp(-ex);
236 }
237
238  Int_t AliGenPHOSlib::IpKaon(TRandom *ran)
239 {
240 //                 particle composition
241 //
242
243     Float_t random = ran->Rndm();
244     Float_t random2 = ran->Rndm();
245     if (random2 < 0.5) 
246     {
247       if (random < 0.5) {       
248         return  321;   //   K+
249       } else {
250         return -321;   // K-
251       }
252     }
253     else
254     {  
255       if (random < 0.5) {       
256         return  130;   // K^0 short
257       } else {  
258         return  310;   // K^0 long
259       }
260     }
261 }
262 // End Kaons
263 //============================================================================
264 //============================================================================
265 //   E  T  A  S
266  Double_t AliGenPHOSlib::PtEta( Double_t *px, Double_t *)
267 {
268 //                etas
269 //                pt-distribution
270 //____________________________________________________________
271
272   return PtScal(*px,3);  //  3==> Eta in the PtScal function
273 }
274
275  Double_t AliGenPHOSlib::YEta( Double_t *py, Double_t *)
276 {
277 // y-distribution
278 //____________________________________________________________
279
280   const Double_t ka    = 1000.;
281   const Double_t kdy   = 4.;
282
283
284   Double_t y=TMath::Abs(*py);
285   //
286   Double_t ex = y*y/(2*kdy*kdy);
287   return ka*TMath::Exp(-ex);
288 }
289
290  Int_t AliGenPHOSlib::IpEta(TRandom *)
291 {
292 //                 particle composition
293 //
294
295         return  221;   //   eta
296 }
297 // End Etas
298
299 //======================================================================
300 //    Eta Flat Distribution
301 //    Transverse momentum distribution PtEtaFlat
302 //    Rapidity distribution YEtaFlat
303 //    Particle distribution IdEtaFlat  111 (pi0)
304 //
305
306 Double_t AliGenPHOSlib::PtEtaFlat(Double_t *px, Double_t *)
307 {
308 //     Eta transverse momentum flat distribution 
309
310   return 1;
311
312 }
313
314 Double_t AliGenPHOSlib::YEtaFlat( Double_t *py, Double_t *)
315 {
316 //
317 // pion y-distribution
318 //
319   return 1.;
320 }
321
322  Int_t AliGenPHOSlib::IpEtaFlat(TRandom *)
323 {
324 //
325 //                 particle composition eta
326 //
327         return 221 ;
328 }
329 // End Pi0Flat
330 //============================================================================
331 //============================================================================
332 //    O  M  E  G  A  S
333  Double_t AliGenPHOSlib::PtOmega( Double_t *px, Double_t *)
334 {
335 // omegas
336 //                pt-distribution
337 //____________________________________________________________
338
339   return PtScal(*px,4);  //  4==> Omega in the PtScal function
340 }
341
342  Double_t AliGenPHOSlib::YOmega( Double_t *py, Double_t *)
343 {
344 // y-distribution
345 //____________________________________________________________
346
347   const Double_t ka    = 1000.;
348   const Double_t kdy   = 4.;
349
350
351   Double_t y=TMath::Abs(*py);
352   //
353   Double_t ex = y*y/(2*kdy*kdy);
354   return ka*TMath::Exp(-ex);
355 }
356
357  Int_t AliGenPHOSlib::IpOmega(TRandom *)
358 {
359 //                 particle composition
360 //
361
362         return  223;   // Omega
363 }
364 // End Omega
365 //============================================================================
366 //============================================================================
367 //    E  T  A  P  R  I  M  E
368  Double_t AliGenPHOSlib::PtEtaprime( Double_t *px, Double_t *)
369 {
370 // etaprime
371 //                pt-distribution
372 //____________________________________________________________
373
374   return PtScal(*px,5);  //  5==> Etaprime in the PtScal function
375 }
376
377  Double_t AliGenPHOSlib::YEtaprime( Double_t *py, Double_t *)
378 {
379 // y-distribution
380 //____________________________________________________________
381
382   const Double_t ka    = 1000.;
383   const Double_t kdy   = 4.;
384
385
386   Double_t y=TMath::Abs(*py);
387   //
388   Double_t ex = y*y/(2*kdy*kdy);
389   return ka*TMath::Exp(-ex);
390 }
391
392  Int_t AliGenPHOSlib::IpEtaprime(TRandom *)
393 {
394 //                 particle composition
395 //
396
397         return  331;   //   Etaprime
398 }
399 // End EtaPrime
400 //===================================================================
401 //============================================================================
402 //    P  H  I   S
403  Double_t AliGenPHOSlib::PtPhi( Double_t *px, Double_t *)
404 {
405 // phi
406 //                pt-distribution
407 //____________________________________________________________
408
409   return PtScal(*px,6);  //  6==> Phi in the PtScal function
410 }
411
412  Double_t AliGenPHOSlib::YPhi( Double_t *py, Double_t *)
413 {
414 // y-distribution
415 //____________________________________________________________
416
417   const Double_t ka    = 1000.;
418   const Double_t kdy   = 4.;
419
420
421   Double_t y=TMath::Abs(*py);
422   //
423   Double_t ex = y*y/(2*kdy*kdy);
424   return ka*TMath::Exp(-ex);
425 }
426
427  Int_t AliGenPHOSlib::IpPhi(TRandom *)
428 {
429 //                 particle composition
430 //
431     
432         return  333;   //   Phi      
433 }
434 // End Phis
435 //===================================================================
436 //============================================================================
437 //    B  A  R  Y  O  N  S  == protons, protonsbar, neutrons, and neutronsbars
438  Double_t AliGenPHOSlib::PtBaryon( Double_t *px, Double_t *)
439 {
440 // baryons
441 //                pt-distribution
442 //____________________________________________________________
443
444   return PtScal(*px,7);  //  7==> Baryon in the PtScal function
445 }
446
447  Double_t AliGenPHOSlib::YBaryon( Double_t *py, Double_t *)
448 {
449 // y-distribution
450 //____________________________________________________________
451
452   const Double_t ka    = 1000.;
453   const Double_t kdy   = 4.;
454
455
456   Double_t y=TMath::Abs(*py);
457   //
458   Double_t ex = y*y/(2*kdy*kdy);
459   return ka*TMath::Exp(-ex);
460 }
461
462  Int_t AliGenPHOSlib::IpBaryon(TRandom *ran)
463 {
464 //                 particle composition
465 //
466
467     Float_t random = ran->Rndm();
468     Float_t random2 = ran->Rndm();
469     if (random2 < 0.5) 
470     {
471       if (random < 0.5) {       
472         return  2212;   //   p
473       } else {
474         return -2212;   // pbar
475       }
476     }
477     else
478     {  
479       if (random < 0.5) {       
480         return  2112;   // n
481       } else {  
482         return -2112;   // n bar
483       }
484     }
485 }
486 // End Baryons
487 //===================================================================
488
489
490 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
491  GenFunc AliGenPHOSlib::GetPt(Int_t param, const char* tname) const
492 {
493 // Return pinter to pT parameterisation
494     GenFunc func;
495     
496     switch (param)
497     {
498     case kPion:     
499         func=PtPion;
500         break;
501     case kPi0Flat:     
502         func=PtPi0Flat;
503         break;
504     case kKaon:
505         func=PtKaon;
506         break;
507     case kEta:
508         func=PtEta;
509         break;
510     case kEtaFlat:
511         func=PtEtaFlat;
512         break;
513     case kOmega:
514         func=PtOmega;
515         break;
516     case kEtaPrime:
517         func=PtEtaprime;
518         break;
519     case kBaryon:
520         func=PtBaryon;
521         break;
522     default:
523         func=0;
524         printf("<AliGenPHOSlib::GetPt> unknown parametrisationn");
525     }
526     return func;
527 }
528
529  GenFunc AliGenPHOSlib::GetY(Int_t param, const char* tname) const
530 {
531 // Return pointer to Y parameterisation
532     GenFunc func;
533     switch (param)
534     {
535     case kPion:
536         func=YPion;
537         break;
538     case kPi0Flat:
539         func=YPi0Flat;
540         break;
541     case kKaon:
542         func=YKaon;
543         break;
544     case kEta:
545         func=YEta;
546         break;
547    case kEtaFlat:
548         func=YEtaFlat;
549         break;
550     case kOmega:
551         func=YOmega;
552         break;
553     case kEtaPrime:
554         func=YEtaprime;
555         break;
556     case kPhi:
557         func=YPhi;
558         break;
559     case kBaryon:
560         func=YBaryon;
561         break;
562     default:
563         func=0;
564         printf("<AliGenPHOSlib::GetY> unknown parametrisationn");
565     }
566     return func;
567 }
568 typedef Int_t (*GenFuncIp) (TRandom *);
569  GenFuncIp AliGenPHOSlib::GetIp(Int_t param,  const char* tname) const
570 {
571 // Return pointer to particle composition
572     GenFuncIp func;
573     switch (param)
574     {
575     case kPion:       
576         func=IpPion;
577         break;
578     case kPi0Flat:       
579         func=IpPi0Flat;
580         break;
581     case kKaon:
582         func=IpKaon;
583         break;
584     case kEta:
585         func=IpEta;
586         break;
587     case kEtaFlat:
588         func=IpEtaFlat;
589         break;
590         
591     case kOmega:
592         func=IpOmega;
593         break;
594     case kEtaPrime:
595         func=IpEtaprime;
596         break;
597     case kPhi:
598         func=IpPhi;
599         break;
600     case kBaryon:
601         func=IpBaryon;
602         break;
603     default:
604         func=0;
605         printf("<AliGenPHOSlib::GetIp> unknown parametrisationn");
606     }
607     return func;
608 }
609