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