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