ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KineticTrack.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4KineticTrack.hh
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 // -----------------------------------------------------------------------------
29 // GEANT 4 class header file
30 //
31 // History: first implementation, A. Feliciello, 20th May 1998
32 // -----------------------------------------------------------------------------
33 
34 #ifndef G4KineticTrack_h
35 #define G4KineticTrack_h 1
36 
37 #include <CLHEP/Units/PhysicalConstants.h>
38 
39 #include "globals.hh"
40 #include "G4ios.hh"
41 
42 
43 #include "Randomize.hh"
44 #include "G4ThreeVector.hh"
45 #include "G4LorentzVector.hh"
46 #include "G4VKineticNucleon.hh"
47 #include "G4Nucleon.hh"
48 #include "G4ParticleDefinition.hh"
49 #include "G4VDecayChannel.hh"
50 #include "G4Log.hh"
51 
52 // #include "G4Allocator.hh"
53 
54 class G4KineticTrackVector;
55 
56 class G4KineticTrack : public G4VKineticNucleon
57 {
58  public:
59 
60  G4KineticTrack();
61 
62  G4KineticTrack(const G4KineticTrack& right);
63 
64  G4KineticTrack(const G4ParticleDefinition* aDefinition,
65  G4double aFormationTime,
66  const G4ThreeVector& aPosition,
67  const G4LorentzVector& a4Momentum);
68  G4KineticTrack(G4Nucleon * nucleon,
69  const G4ThreeVector& aPosition,
70  const G4LorentzVector& a4Momentum);
71 
72  ~G4KineticTrack();
73 
74  G4KineticTrack& operator=(const G4KineticTrack& right);
75 
76  G4bool operator==(const G4KineticTrack& right) const;
77 
78  G4bool operator!=(const G4KineticTrack& right) const;
79 /*
80  inline void *operator new(size_t);
81  inline void operator delete(void *aTrack);
82 */
83  const G4ParticleDefinition* GetDefinition() const;
84  void SetDefinition(const G4ParticleDefinition* aDefinition);
85 
86  G4double GetFormationTime() const;
87  void SetFormationTime(G4double aFormationTime);
88 
89  const G4ThreeVector& GetPosition() const;
90  void SetPosition(const G4ThreeVector aPosition);
91 
92  const G4LorentzVector& Get4Momentum() const;
93  void Set4Momentum(const G4LorentzVector& a4Momentum);
94  void Update4Momentum(G4double aEnergy); // update E and p, not changing mass
95  void Update4Momentum(const G4ThreeVector & aMomentum); // idem
96  void SetTrackingMomentum(const G4LorentzVector& a4Momentum);
97  void UpdateTrackingMomentum(G4double aEnergy); // update E and p, not changing mass
98  void UpdateTrackingMomentum(const G4ThreeVector & aMomentum); // idem
99 
100  const G4LorentzVector& GetTrackingMomentum() const;
101 
102  G4double SampleResidualLifetime();
103 
104  void Hit();
105  void SetNucleon(G4Nucleon * aN) {theNucleon = aN;}
106 
107  G4bool IsParticipant() const;
108 
109  G4KineticTrackVector* Decay();
110 
111  // LB move to public (before was private) LB
112  G4double* GetActualWidth() const;
113 
114  G4double GetActualMass() const;
115  G4int GetnChannels() const;
116 
117 // position relativ to nucleus "state"
118  enum CascadeState {undefined, outside, going_in, inside,
119  going_out, gone_out, captured, miss_nucleus };
120 
121  CascadeState SetState(const CascadeState new_state);
122  CascadeState GetState() const;
123  void SetProjectilePotential(const G4double aPotential);
124  G4double GetProjectilePotential() const;
125 
126 
127  private:
128 
129 
130  void SetnChannels(const G4int aChannel);
131 
132  void SetActualWidth(G4double* anActualWidth);
133 
134  G4double EvaluateTotalActualWidth();
135 
136  G4double EvaluateCMMomentum (const G4double mass,
137  const G4double* m_ij) const;
138 
139  G4double IntegrateCMMomentum(const G4double lowerLimit) const;
140 
141  G4double IntegrateCMMomentum(const G4double lowerLimit ,const G4double polemass) const;
142 
143  G4double IntegrateCMMomentum2() const;
144 
145  public:
146 
147  G4double BrWig(const G4double Gamma,
148  const G4double rmass,
149  const G4double mass) const;
150 
151 private:
152  G4double IntegrandFunction1 (G4double xmass) const;
153  G4double IntegrandFunction2 (G4double xmass) const;
154  G4double IntegrandFunction3 (G4double xmass) const;
155  G4double IntegrandFunction4 (G4double xmass) const;
156 public:
157  // friend G4double IntegrandFunction3 (G4double xmass);
158 
159  // friend G4double IntegrandFunction4 (G4double xmass);
160 
161  private:
162 
163  const G4ParticleDefinition* theDefinition;
164 
165  G4double theFormationTime;
166 
167  G4ThreeVector thePosition;
168 
169  G4LorentzVector the4Momentum;
170  G4LorentzVector theFermi3Momentum;
171  G4LorentzVector theTotal4Momentum;
172 
173  G4Nucleon * theNucleon;
174 
175  G4int nChannels;
176 
177  G4double theActualMass;
178 
179  G4double* theActualWidth;
180 
181  // Temporary storage for daughter masses and widths
182  // (needed because Integrand Function cannot take > 1 argument)
183  G4double* theDaughterMass;
184  G4double* theDaughterWidth;
185 
186  CascadeState theStateToNucleus;
187 
188  G4double theProjectilePotential;
189 };
190 
191 // extern G4Allocator<G4KineticTrack> theKTAllocator;
192 
193 
194 // Class G4KineticTrack
195 /*
196 inline void * G4KineticTrack::operator new(size_t)
197 {
198  void * aT;
199  aT = (void *) theKTAllocator.MallocSingle();
200  return aT;
201 }
202 
203 inline void G4KineticTrack::operator delete(void * aT)
204 {
205  theKTAllocator.FreeSingle((G4KineticTrack *) aT);
206 }
207 */
208 
209 inline const G4ParticleDefinition* G4KineticTrack::GetDefinition() const
210 {
211  return theDefinition;
212 }
213 
214 inline void G4KineticTrack::SetDefinition(const G4ParticleDefinition* aDefinition)
215 {
216  theDefinition = aDefinition;
217 }
218 
219 inline G4double G4KineticTrack::GetFormationTime() const
220 {
221  return theFormationTime;
222 }
223 
224 inline void G4KineticTrack::SetFormationTime(G4double aFormationTime)
225 {
226  theFormationTime = aFormationTime;
227 }
228 
229 inline const G4ThreeVector& G4KineticTrack::GetPosition() const
230 {
231  return thePosition;
232 }
233 
234 inline void G4KineticTrack::SetPosition(const G4ThreeVector aPosition)
235 {
236  thePosition = aPosition;
237 }
238 
239 inline const G4LorentzVector& G4KineticTrack::Get4Momentum() const
240 {
241  return theTotal4Momentum;
242 }
243 
244 inline const G4LorentzVector& G4KineticTrack::GetTrackingMomentum() const
245 {
246  return the4Momentum;
247 }
248 
249 inline void G4KineticTrack::Set4Momentum(const G4LorentzVector& a4Momentum)
250 {
251 // set the4Momentum and update theTotal4Momentum
252 
253  theTotal4Momentum=a4Momentum;
254  the4Momentum = theTotal4Momentum;
255  theFermi3Momentum=G4LorentzVector(0);
256 }
257 
258 inline void G4KineticTrack::Update4Momentum(G4double aEnergy)
259 {
260 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
261 // updates theTotal4Momentum as well.
262  G4double newP(0);
263  G4double mass2=theTotal4Momentum.mag2();
264  if ( sqr(aEnergy) > mass2 )
265  {
266  newP = std::sqrt(sqr(aEnergy) - mass2 );
267  } else
268  {
269  aEnergy=std::sqrt(mass2);
270  }
271  Set4Momentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
272 }
273 
274 inline void G4KineticTrack::Update4Momentum(const G4ThreeVector & aMomentum)
275 {
276 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
277 // updates theTotal4Momentum as well.
278  G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
279  Set4Momentum(G4LorentzVector(aMomentum, newE));
280 }
281 
282 inline void G4KineticTrack::SetTrackingMomentum(const G4LorentzVector& aMomentum)
283 {
284 // set the4Momentum and update theTotal4Momentum, keep the mass of aMomentum
285 
286  the4Momentum = aMomentum;
287  theTotal4Momentum=the4Momentum+theFermi3Momentum;
288 // keep mass of aMomentum for the total momentum
289  G4double mass2 = aMomentum.mag2();
290  G4double p2=theTotal4Momentum.vect().mag2();
291  theTotal4Momentum.setE(std::sqrt(mass2+p2));
292 }
293 
294 inline void G4KineticTrack::UpdateTrackingMomentum(G4double aEnergy)
295 {
296 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
297 // updates theTotal4Momentum as well.
298  G4double newP(0);
299  G4double mass2=theTotal4Momentum.mag2();
300  if ( sqr(aEnergy) > mass2 )
301  {
302  newP = std::sqrt(sqr(aEnergy) - mass2 );
303  } else
304  {
305  aEnergy=std::sqrt(mass2);
306  }
307  SetTrackingMomentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
308 }
309 
310 inline void G4KineticTrack::UpdateTrackingMomentum(const G4ThreeVector & aMomentum)
311 {
312 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
313 // updates theTotal4Momentum as well.
314  G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
315  SetTrackingMomentum(G4LorentzVector(aMomentum, newE));
316 }
317 
318 inline G4double G4KineticTrack::GetActualMass() const
319 {
320  return std::sqrt(std::abs(the4Momentum.mag2()));
321 }
322 
323 inline G4int G4KineticTrack::GetnChannels() const
324 {
325  return nChannels;
326 }
327 
328 inline void G4KineticTrack::SetnChannels(const G4int numberOfChannels)
329 {
330  nChannels = numberOfChannels;
331 }
332 
333 inline G4double* G4KineticTrack::GetActualWidth() const
334 {
335  return theActualWidth;
336 }
337 
338 inline void G4KineticTrack::SetActualWidth(G4double* anActualWidth)
339 {
340  theActualWidth = anActualWidth;
341 }
342 
343 
344 
345 inline G4double G4KineticTrack::EvaluateTotalActualWidth()
346 {
347  G4int index;
348  G4double theTotalActualWidth = 0.0;
349  for (index = nChannels - 1; index >= 0; index--)
350  {
351  theTotalActualWidth += theActualWidth[index];
352  }
353  return theTotalActualWidth;
354 }
355 
356 inline G4double G4KineticTrack::SampleResidualLifetime()
357 {
358  G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
359  G4double tau = CLHEP::hbar_Planck * (-1.0 / theTotalActualWidth);
360  G4double theResidualLifetime = tau * G4Log(G4UniformRand());
361  return theResidualLifetime*the4Momentum.gamma();
362 }
363 
364 inline G4double G4KineticTrack::EvaluateCMMomentum(const G4double mass,
365  const G4double* m_ij) const
366 {
367  G4double theCMMomentum;
368  if((m_ij[0]+m_ij[1])<mass)
369  theCMMomentum = 1 / (2 * mass) *
370  std::sqrt (((mass * mass) - (m_ij[0] + m_ij[1]) * (m_ij[0] + m_ij[1])) *
371  ((mass * mass) - (m_ij[0] - m_ij[1]) * (m_ij[0] - m_ij[1])));
372  else
373  theCMMomentum=0.;
374 
375  return theCMMomentum;
376 }
377 
378 inline G4double G4KineticTrack::BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
379 {
380  G4double Norm = CLHEP::twopi;
381  return (Gamma/((mass-rmass)*(mass-rmass)+Gamma*Gamma/4.))/Norm;
382 }
383 
384 inline
385 void G4KineticTrack::Hit()
386 {
387  if(theNucleon)
388  {
389  theNucleon->Hit(1);
390  }
391 }
392 
393 inline
394 G4bool G4KineticTrack::IsParticipant() const
395 {
396  if(!theNucleon) return true;
397  return theNucleon->AreYouHit();
398 }
399 
400 inline
401 G4KineticTrack::CascadeState G4KineticTrack::GetState() const
402 {
403  return theStateToNucleus;
404 }
405 
406 inline
407 G4KineticTrack::CascadeState G4KineticTrack::SetState(const CascadeState new_state)
408 {
409  CascadeState old_state=theStateToNucleus;
410  theStateToNucleus=new_state;
411  return old_state;
412 }
413 
414 inline
415 void G4KineticTrack::SetProjectilePotential(G4double aPotential)
416 {
417  theProjectilePotential = aPotential;
418 }
419 inline
420 G4double G4KineticTrack::GetProjectilePotential() const
421 {
422  return theProjectilePotential;
423 }
424 
425 #endif
426 
427 
428