Saver calculation of rapdity.
[u/mrichter/AliRoot.git] / EVGEN / AliGenMC.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  2002/03/12 17:02:20  morsch
19 Change in calculation of rapidity, include case in which numerically e == pz.
20
21 Revision 1.4  2001/11/27 13:13:07  morsch
22 Maximum lifetime for long-lived particles to be put on the stack is parameter.
23 It can be set via SetMaximumLifetime(..).
24
25 Revision 1.3  2001/10/16 08:48:56  morsch
26 Common vertex related code moved to base class AliGenerator.
27
28 Revision 1.2  2001/10/15 08:15:51  morsch
29 Event vertex and vertex truncation setting moved into AliMC.
30
31 Revision 1.1  2001/07/13 10:56:00  morsch
32 AliGenMC base class for AliGenParam and AliGenPythia commonalities.
33
34 */
35
36 #include "AliGenMC.h"
37 #include "AliPDG.h"
38 #include <TParticle.h>
39
40  ClassImp(AliGenMC)
41
42 AliGenMC::AliGenMC()
43                  :AliGenerator()
44 {
45 // Default Constructor
46     SetCutOnChild();
47     SetChildMomentumRange();
48     SetChildPtRange();
49     SetChildPhiRange();
50     SetChildThetaRange(); 
51     SetChildYRange(); 
52     SetMaximumLifetime();
53 }
54
55 AliGenMC::AliGenMC(Int_t npart)
56                  :AliGenerator(npart)
57 {
58 //  Constructor
59     SetCutOnChild();
60     SetChildMomentumRange();
61     SetChildPtRange();
62     SetChildPhiRange();
63     SetChildThetaRange();
64     SetChildYRange(); 
65 // 
66     fParentSelect.Set(8);
67     fChildSelect.Set(8);
68     for (Int_t i=0; i<8; i++) fParentSelect[i]=fChildSelect[i]=0;
69     SetMaximumLifetime();
70 }
71
72 AliGenMC::AliGenMC(const AliGenMC & mc)
73 {
74 // copy constructor
75 }
76
77 AliGenMC::~AliGenMC()
78 {
79 // Destructor
80 }
81
82 void AliGenMC::Init()
83 {
84 //
85 //  Initialization
86     switch (fForceDecay) {
87     case kSemiElectronic:
88     case kDiElectron:
89     case kBJpsiDiElectron:
90     case kBPsiPrimeDiElectron:
91         fChildSelect[0] = kElectron;    
92         break;
93     case kSemiMuonic:
94     case kDiMuon:
95     case kBJpsiDiMuon:
96     case kBPsiPrimeDiMuon:
97     case kPiToMu:
98     case kKaToMu:
99         fChildSelect[0]=kMuonMinus;
100         break;
101     case kHadronicD:
102         fChildSelect[0]=kPiPlus;
103         fChildSelect[1]=kKPlus;
104         break;
105     case kAll:
106     case kNoDecay:
107         break;
108     }
109 }
110
111
112 Bool_t AliGenMC::ParentSelected(Int_t ip)
113 {
114 // True if particle is in list of parent particles to be selected
115     for (Int_t i=0; i<8; i++)
116     {
117         if (fParentSelect[i]==ip) return kTRUE;
118     }
119     return kFALSE;
120 }
121
122 Bool_t AliGenMC::ChildSelected(Int_t ip)
123 {
124 // True if particle is in list of decay products to be selected
125     for (Int_t i=0; i<5; i++)
126     {
127         if (fChildSelect[i]==ip) return kTRUE;
128     }
129     return kFALSE;
130 }
131
132 Bool_t AliGenMC::KinematicSelection(TParticle *particle, Int_t flag)
133 {
134 // Perform kinematic selection
135     Float_t px    = particle->Px();
136     Float_t py    = particle->Py();
137     Float_t pz    = particle->Pz();
138     Float_t  e    = particle->Energy();
139     Float_t pt    = particle->Pt();
140     Float_t p     = particle->P();
141     Float_t theta = particle->Theta();
142     Float_t mass  = particle->GetMass();
143     Float_t mt2   = pt * pt + mass * mass;
144     
145     Float_t phi   = Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
146     Double_t y, y0;
147
148     if (TMath::Abs(pz) <  e) {
149         y = 0.5*TMath::Log((e+pz)/(e-pz));
150     } else {
151         y = 1.e10;
152     }
153     
154     if (mt2) {
155         y0 = 0.5*TMath::Log((e+TMath::Abs(pz))*(e+TMath::Abs(pz))/mt2);
156     } else {
157         if (TMath::Abs(y) < 1.e10) {
158             y0 = y;
159         } else {
160             y0 = 1.e10;
161         }
162     }
163       
164     y = (pz < 0) ? -y0 : y0;
165     
166     if (flag == 0) {
167 //
168 // Primary particle cuts
169 //
170 //  transverse momentum cut    
171         if (pt > fPtMax || pt < fPtMin) {
172 //          printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
173             return kFALSE;
174         }
175 //
176 // momentum cut
177         if (p > fPMax || p < fPMin) {
178 //          printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
179             return kFALSE;
180         }
181 //
182 // theta cut
183         if (theta > fThetaMax || theta < fThetaMin) {
184 //          printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
185             return kFALSE;
186         }
187 //
188 // rapidity cut
189         if (y > fYMax || y < fYMin) {
190 //          printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
191             return kFALSE;
192         }
193 //
194 // phi cut
195         if (phi > fPhiMax || phi < fPhiMin) {
196 //          printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
197             return kFALSE;
198         }
199     } else {
200 //
201 // Decay product cuts
202 //
203 //  transverse momentum cut    
204         if (pt > fChildPtMax || pt < fChildPtMin) {
205 //          printf("\n failed pt cut %f %f %f \n",pt,fChildPtMin,fChildPtMax);
206             return kFALSE;
207         }
208 //
209 // momentum cut
210         if (p > fChildPMax || p < fChildPMin) {
211 //          printf("\n failed p cut %f %f %f \n",p,fChildPMin,fChildPMax);
212             return kFALSE;
213         }
214 //
215 // theta cut
216         if (theta > fChildThetaMax || theta < fChildThetaMin) {
217 //          printf("\n failed theta cut %f %f %f \n",theta,fChildThetaMin,fChildThetaMax);
218             return kFALSE;
219         }
220 //
221 // rapidity cut
222         if (y > fChildYMax || y < fChildYMin) {
223 //          printf("\n failed y cut %f %f %f \n",y,fChildYMin,fChildYMax);
224             return kFALSE;
225         }
226 //
227 // phi cut
228         if (phi > fChildPhiMax || phi < fChildPhiMin) {
229 //          printf("\n failed phi cut %f %f %f \n",phi,fChildPhiMin,fChildPhiMax);
230             return kFALSE;
231         }
232     }
233     
234     
235
236     return kTRUE;
237 }
238
239 Int_t AliGenMC::CheckPDGCode(Int_t pdgcode)
240 {
241 //
242 //  If the particle is in a diffractive state, then take action accordingly
243   switch (pdgcode) {
244   case 91:
245     return 92;
246   case 110:
247     //rho_diff0 -- difficult to translate, return rho0
248     return 113;
249   case 210:
250     //pi_diffr+ -- change to pi+
251     return 211;
252   case 220:
253     //omega_di0 -- change to omega0
254     return 223;
255   case 330:
256     //phi_diff0 -- return phi0
257     return 333;
258   case 440:
259     //J/psi_di0 -- return J/psi
260     return 443;
261   case 2110:
262     //n_diffr -- return neutron
263     return 2112;
264   case 2210:
265     //p_diffr+ -- return proton
266     return 2212;
267   }
268   //non diffractive state -- return code unchanged
269   return pdgcode;
270 }
271           
272 AliGenMC& AliGenMC::operator=(const  AliGenMC& rhs)
273 {
274 // Assignment operator
275     return *this;
276 }
277