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