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