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