Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KineticTrack.hh
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // $Id: G4KineticTrack.hh,v 1.0 1998/05/20
29 // -----------------------------------------------------------------------------
30 // GEANT 4 class header file
31 //
32 // History: first implementation, A. Feliciello, 20th May 1998
33 // -----------------------------------------------------------------------------
34 
35 #ifndef G4KineticTrack_h
36 #define G4KineticTrack_h 1
37 
39 
40 #include "globals.hh"
41 #include "G4ios.hh"
42 
43 
44 #include "Randomize.hh"
45 #include "G4ThreeVector.hh"
46 #include "G4LorentzVector.hh"
47 #include "G4VKineticNucleon.hh"
48 #include "G4Nucleon.hh"
49 #include "G4ParticleDefinition.hh"
50 #include "G4VDecayChannel.hh"
51 #include "G4Log.hh"
52 
53 // #include "G4Allocator.hh"
54 
56 
58 {
59  public:
60 
62 
64 
65  G4KineticTrack(const G4ParticleDefinition* aDefinition,
66  G4double aFormationTime,
67  const G4ThreeVector& aPosition,
68  const G4LorentzVector& a4Momentum);
70  const G4ThreeVector& aPosition,
71  const G4LorentzVector& a4Momentum);
72 
74 
76 
77  G4int operator==(const G4KineticTrack& right) const;
78 
79  G4int operator!=(const G4KineticTrack& right) const;
80 /*
81  inline void *operator new(size_t);
82  inline void operator delete(void *aTrack);
83 */
84  const G4ParticleDefinition* GetDefinition() const;
85  void SetDefinition(const G4ParticleDefinition* aDefinition);
86 
87  G4double GetFormationTime() const;
88  void SetFormationTime(G4double aFormationTime);
89 
90  const G4ThreeVector& GetPosition() const;
91  void SetPosition(const G4ThreeVector aPosition);
92 
93  const G4LorentzVector& Get4Momentum() const;
94  void Set4Momentum(const G4LorentzVector& a4Momentum);
95  void Update4Momentum(G4double aEnergy); // update E and p, not changing mass
96  void Update4Momentum(const G4ThreeVector & aMomentum); // idem
97  void SetTrackingMomentum(const G4LorentzVector& a4Momentum);
98  void UpdateTrackingMomentum(G4double aEnergy); // update E and p, not changing mass
99  void UpdateTrackingMomentum(const G4ThreeVector & aMomentum); // idem
100 
101  const G4LorentzVector& GetTrackingMomentum() const;
102 
104 
105  void Hit();
106  void SetNucleon(G4Nucleon * aN) {theNucleon = aN;}
107 
108  G4bool IsParticipant() const;
109 
111 
112  // LB move to public (before was private) LB
113  G4double* GetActualWidth() const;
114 
115  G4double GetActualMass() const;
116  G4int GetnChannels() const;
117 
118 // position relativ to nucleus "state"
121 
122  CascadeState SetState(const CascadeState new_state);
123  CascadeState GetState() const;
124  void SetProjectilePotential(const G4double aPotential);
126 
127 
128  private:
129 
130 
131  void SetnChannels(const G4int aChannel);
132 
133  void SetActualWidth(G4double* anActualWidth);
134 
135  G4double EvaluateTotalActualWidth();
136 
137  G4double EvaluateCMMomentum (const G4double mass,
138  const G4double* m_ij) const;
139 
140  G4double IntegrateCMMomentum(const G4double lowerLimit) const;
141 
142  G4double IntegrateCMMomentum(const G4double lowerLimit ,const G4double polemass) const;
143 
144  G4double IntegrateCMMomentum2() const;
145 
146  public:
147 
148  G4double BrWig(const G4double Gamma,
149  const G4double rmass,
150  const G4double mass) const;
151 
152 private:
153  G4double IntegrandFunction1 (G4double xmass) const;
154  G4double IntegrandFunction2 (G4double xmass) const;
155  G4double IntegrandFunction3 (G4double xmass) const;
156  G4double IntegrandFunction4 (G4double xmass) const;
157 public:
158  // friend G4double IntegrandFunction3 (G4double xmass);
159 
160  // friend G4double IntegrandFunction4 (G4double xmass);
161 
162  private:
163 
164  const G4ParticleDefinition* theDefinition;
165 
166  G4double theFormationTime;
167 
168  G4ThreeVector thePosition;
169 
170  G4LorentzVector the4Momentum;
171  G4LorentzVector theFermi3Momentum;
172  G4LorentzVector theTotal4Momentum;
173 
174  G4Nucleon * theNucleon;
175 
176  G4int nChannels;
177 
178  G4double theActualMass;
179 
180  G4double* theActualWidth;
181 
182  // Temporary storage for daughter masses and widths
183  // (needed because Integrand Function cannot take > 1 argument)
184  G4double* theDaughterMass;
185  G4double* theDaughterWidth;
186 
187  CascadeState theStateToNucleus;
188 
189  G4double theProjectilePotential;
190 };
191 
192 // extern G4Allocator<G4KineticTrack> theKTAllocator;
193 
194 
195 // Class G4KineticTrack
196 /*
197 inline void * G4KineticTrack::operator new(size_t)
198 {
199  void * aT;
200  aT = (void *) theKTAllocator.MallocSingle();
201  return aT;
202 }
203 
204 inline void G4KineticTrack::operator delete(void * aT)
205 {
206  theKTAllocator.FreeSingle((G4KineticTrack *) aT);
207 }
208 */
209 
211 {
212  return theDefinition;
213 }
214 
216 {
217  theDefinition = aDefinition;
218 }
219 
221 {
222  return theFormationTime;
223 }
224 
225 inline void G4KineticTrack::SetFormationTime(G4double aFormationTime)
226 {
227  theFormationTime = aFormationTime;
228 }
229 
231 {
232  return thePosition;
233 }
234 
235 inline void G4KineticTrack::SetPosition(const G4ThreeVector aPosition)
236 {
237  thePosition = aPosition;
238 }
239 
241 {
242  return theTotal4Momentum;
243 }
244 
246 {
247  return the4Momentum;
248 }
249 
250 inline void G4KineticTrack::Set4Momentum(const G4LorentzVector& a4Momentum)
251 {
252 // set the4Momentum and update theTotal4Momentum
253 
254  theTotal4Momentum=a4Momentum;
255  the4Momentum = theTotal4Momentum;
256  theFermi3Momentum=G4LorentzVector(0);
257 }
258 
260 {
261 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
262 // updates theTotal4Momentum as well.
263  G4double newP(0);
264  G4double mass2=theTotal4Momentum.mag2();
265  if ( sqr(aEnergy) > mass2 )
266  {
267  newP = std::sqrt(sqr(aEnergy) - mass2 );
268  } else
269  {
270  aEnergy=std::sqrt(mass2);
271  }
272  Set4Momentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
273 }
274 
275 inline void G4KineticTrack::Update4Momentum(const G4ThreeVector & aMomentum)
276 {
277 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
278 // updates theTotal4Momentum as well.
279  G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
280  Set4Momentum(G4LorentzVector(aMomentum, newE));
281 }
282 
284 {
285 // set the4Momentum and update theTotal4Momentum, keep the mass of aMomentum
286 
287  the4Momentum = aMomentum;
288  theTotal4Momentum=the4Momentum+theFermi3Momentum;
289 // keep mass of aMomentum for the total momentum
290  G4double mass2 = aMomentum.mag2();
291  G4double p2=theTotal4Momentum.vect().mag2();
292  theTotal4Momentum.setE(std::sqrt(mass2+p2));
293 }
294 
296 {
297 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
298 // updates theTotal4Momentum as well.
299  G4double newP(0);
300  G4double mass2=theTotal4Momentum.mag2();
301  if ( sqr(aEnergy) > mass2 )
302  {
303  newP = std::sqrt(sqr(aEnergy) - mass2 );
304  } else
305  {
306  aEnergy=std::sqrt(mass2);
307  }
308  SetTrackingMomentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
309 }
310 
312 {
313 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
314 // updates theTotal4Momentum as well.
315  G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
316  SetTrackingMomentum(G4LorentzVector(aMomentum, newE));
317 }
318 
320 {
321  return std::sqrt(std::abs(the4Momentum.mag2()));
322 }
323 
325 {
326  return nChannels;
327 }
328 
329 inline void G4KineticTrack::SetnChannels(const G4int numberOfChannels)
330 {
331  nChannels = numberOfChannels;
332 }
333 
335 {
336  return theActualWidth;
337 }
338 
339 inline void G4KineticTrack::SetActualWidth(G4double* anActualWidth)
340 {
341  theActualWidth = anActualWidth;
342 }
343 
344 
345 
346 inline G4double G4KineticTrack::EvaluateTotalActualWidth()
347 {
348  G4int index;
349  G4double theTotalActualWidth = 0.0;
350  for (index = nChannels - 1; index >= 0; index--)
351  {
352  theTotalActualWidth += theActualWidth[index];
353  }
354  return theTotalActualWidth;
355 }
356 
358 {
359  G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
360  G4double tau = CLHEP::hbar_Planck * (-1.0 / theTotalActualWidth);
361  G4double theResidualLifetime = tau * G4Log(G4UniformRand());
362  return theResidualLifetime*the4Momentum.gamma();
363 }
364 
365 inline G4double G4KineticTrack::EvaluateCMMomentum(const G4double mass,
366  const G4double* m_ij) const
367 {
368  G4double theCMMomentum;
369  if((m_ij[0]+m_ij[1])<mass)
370  theCMMomentum = 1 / (2 * mass) *
371  std::sqrt (((mass * mass) - (m_ij[0] + m_ij[1]) * (m_ij[0] + m_ij[1])) *
372  ((mass * mass) - (m_ij[0] - m_ij[1]) * (m_ij[0] - m_ij[1])));
373  else
374  theCMMomentum=0.;
375 
376  return theCMMomentum;
377 }
378 
379 inline G4double G4KineticTrack::BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
380 {
381  G4double Norm = CLHEP::twopi;
382  return (Gamma/((mass-rmass)*(mass-rmass)+Gamma*Gamma/4.))/Norm;
383 }
384 
385 inline
387 {
388  if(theNucleon)
389  {
390  theNucleon->Hit(1);
391  }
392 }
393 
394 inline
396 {
397  if(!theNucleon) return true;
398  return theNucleon->AreYouHit();
399 }
400 
401 inline
403 {
404  return theStateToNucleus;
405 }
406 
407 inline
409 {
410  CascadeState old_state=theStateToNucleus;
411  theStateToNucleus=new_state;
412  return old_state;
413 }
414 
415 inline
417 {
418  theProjectilePotential = aPotential;
419 }
420 inline
422 {
423  return theProjectilePotential;
424 }
425 
426 #endif
427 
428 
429 
void SetDefinition(const G4ParticleDefinition *aDefinition)
G4bool IsParticipant() const
void Update4Momentum(G4double aEnergy)
G4int operator==(const G4KineticTrack &right) const
CascadeState GetState() const
G4int GetnChannels() const
const G4ThreeVector & GetPosition() const
G4KineticTrackVector * Decay()
void SetFormationTime(G4double aFormationTime)
void UpdateTrackingMomentum(G4double aEnergy)
G4double GetActualMass() const
void SetProjectilePotential(const G4double aPotential)
int G4int
Definition: G4Types.hh:78
G4double * GetActualWidth() const
G4int operator!=(const G4KineticTrack &right) const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
CascadeState SetState(const CascadeState new_state)
G4bool nucleon(G4int ityp)
G4double GetFormationTime() const
bool G4bool
Definition: G4Types.hh:79
void SetNucleon(G4Nucleon *aN)
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
void SetPosition(const G4ThreeVector aPosition)
G4KineticTrack & operator=(const G4KineticTrack &right)
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
const G4LorentzVector & GetTrackingMomentum() const
Hep3Vector unit() const
double mag2() const
double mag2() const
void SetTrackingMomentum(const G4LorentzVector &a4Momentum)
void Hit(G4VSplitableHadron *aHit)
Definition: G4Nucleon.hh:90
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
const G4LorentzVector & Get4Momentum() const
G4double SampleResidualLifetime()
static constexpr double hbar_Planck
const G4ParticleDefinition * GetDefinition() const
G4double BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4double GetProjectilePotential() const
CLHEP::HepLorentzVector G4LorentzVector