| Geant4 Cross Reference |
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 void SetCreatorModelID(G4int id);
127 G4int GetCreatorModelID() const;
128
129 const G4ParticleDefinition* GetParentResonanceDef() const;
130 void SetParentResonanceDef(const G4ParticleDefinition* parent);
131 G4int GetParentResonanceID() const;
132 void SetParentResonanceID(const G4int parentID);
133
134 private:
135
136
137 void SetnChannels(const G4int aChannel);
138
139 void SetActualWidth(G4double* anActualWidth);
140
141 G4double EvaluateTotalActualWidth();
142
143 G4double EvaluateCMMomentum (const G4double mass,
144 const G4double* m_ij) const;
145
146 G4double IntegrateCMMomentum(const G4double lowerLimit) const;
147
148 G4double IntegrateCMMomentum(const G4double lowerLimit ,const G4double polemass) const;
149
150 G4double IntegrateCMMomentum2() const;
151
152 public:
153
154 G4double BrWig(const G4double Gamma,
155 const G4double rmass,
156 const G4double mass) const;
157
158 private:
159 G4double IntegrandFunction1 (G4double xmass) const;
160 G4double IntegrandFunction2 (G4double xmass) const;
161 G4double IntegrandFunction3 (G4double xmass) const;
162 G4double IntegrandFunction4 (G4double xmass) const;
163 public:
164 // friend G4double IntegrandFunction3 (G4double xmass);
165
166 // friend G4double IntegrandFunction4 (G4double xmass);
167
168 private:
169
170 const G4ParticleDefinition* theDefinition;
171
172 G4double theFormationTime;
173
174 G4ThreeVector thePosition;
175
176 G4LorentzVector the4Momentum;
177 G4LorentzVector theFermi3Momentum;
178 G4LorentzVector theTotal4Momentum;
179
180 G4Nucleon * theNucleon;
181
182 G4int nChannels;
183
184 G4double theActualMass;
185
186 G4double* theActualWidth;
187
188 // Temporary storage for daughter masses and widths
189 // (needed because Integrand Function cannot take > 1 argument)
190 G4double* theDaughterMass;
191 G4double* theDaughterWidth;
192
193 CascadeState theStateToNucleus;
194
195 G4double theProjectilePotential;
196
197 G4int theCreatorModel;
198
199 const G4ParticleDefinition* theParentResonanceDef = nullptr;
200 G4int theParentResonanceID;
201 };
202
203 // extern G4Allocator<G4KineticTrack> theKTAllocator;
204
205
206 // Class G4KineticTrack
207 /*
208 inline void * G4KineticTrack::operator new(size_t)
209 {
210 void * aT;
211 aT = (void *) theKTAllocator.MallocSingle();
212 return aT;
213 }
214
215 inline void G4KineticTrack::operator delete(void * aT)
216 {
217 theKTAllocator.FreeSingle((G4KineticTrack *) aT);
218 }
219 */
220
221 inline const G4ParticleDefinition* G4KineticTrack::GetDefinition() const
222 {
223 return theDefinition;
224 }
225
226 inline void G4KineticTrack::SetDefinition(const G4ParticleDefinition* aDefinition)
227 {
228 theDefinition = aDefinition;
229 }
230
231 inline G4double G4KineticTrack::GetFormationTime() const
232 {
233 return theFormationTime;
234 }
235
236 inline void G4KineticTrack::SetFormationTime(G4double aFormationTime)
237 {
238 theFormationTime = aFormationTime;
239 }
240
241 inline const G4ThreeVector& G4KineticTrack::GetPosition() const
242 {
243 return thePosition;
244 }
245
246 inline void G4KineticTrack::SetPosition(const G4ThreeVector aPosition)
247 {
248 thePosition = aPosition;
249 }
250
251 inline const G4LorentzVector& G4KineticTrack::Get4Momentum() const
252 {
253 return theTotal4Momentum;
254 }
255
256 inline const G4LorentzVector& G4KineticTrack::GetTrackingMomentum() const
257 {
258 return the4Momentum;
259 }
260
261 inline void G4KineticTrack::Set4Momentum(const G4LorentzVector& a4Momentum)
262 {
263 // set the4Momentum and update theTotal4Momentum
264
265 theTotal4Momentum=a4Momentum;
266 the4Momentum = theTotal4Momentum;
267 theFermi3Momentum=G4LorentzVector(0);
268 }
269
270 inline void G4KineticTrack::Update4Momentum(G4double aEnergy)
271 {
272 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
273 // updates theTotal4Momentum as well.
274 G4double newP(0);
275 G4double mass2=theTotal4Momentum.mag2();
276 if ( sqr(aEnergy) > mass2 )
277 {
278 newP = std::sqrt(sqr(aEnergy) - mass2 );
279 } else
280 {
281 aEnergy=std::sqrt(mass2);
282 }
283 Set4Momentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
284 }
285
286 inline void G4KineticTrack::Update4Momentum(const G4ThreeVector & aMomentum)
287 {
288 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
289 // updates theTotal4Momentum as well.
290 G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
291 Set4Momentum(G4LorentzVector(aMomentum, newE));
292 }
293
294 inline void G4KineticTrack::SetTrackingMomentum(const G4LorentzVector& aMomentum)
295 {
296 // set the4Momentum and update theTotal4Momentum, keep the mass of aMomentum
297
298 the4Momentum = aMomentum;
299 theTotal4Momentum=the4Momentum+theFermi3Momentum;
300 // keep mass of aMomentum for the total momentum
301 G4double mass2 = aMomentum.mag2();
302 G4double p2=theTotal4Momentum.vect().mag2();
303 theTotal4Momentum.setE(std::sqrt(mass2+p2));
304 }
305
306 inline void G4KineticTrack::UpdateTrackingMomentum(G4double aEnergy)
307 {
308 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
309 // updates theTotal4Momentum as well.
310 G4double newP(0);
311 G4double mass2=theTotal4Momentum.mag2();
312 if ( sqr(aEnergy) > mass2 )
313 {
314 newP = std::sqrt(sqr(aEnergy) - mass2 );
315 } else
316 {
317 aEnergy=std::sqrt(mass2);
318 }
319 SetTrackingMomentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
320 }
321
322 inline void G4KineticTrack::UpdateTrackingMomentum(const G4ThreeVector & aMomentum)
323 {
324 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
325 // updates theTotal4Momentum as well.
326 G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
327 SetTrackingMomentum(G4LorentzVector(aMomentum, newE));
328 }
329
330 inline G4double G4KineticTrack::GetActualMass() const
331 {
332 return std::sqrt(std::abs(the4Momentum.mag2()));
333 }
334
335 inline G4int G4KineticTrack::GetnChannels() const
336 {
337 return nChannels;
338 }
339
340 inline void G4KineticTrack::SetnChannels(const G4int numberOfChannels)
341 {
342 nChannels = numberOfChannels;
343 }
344
345 inline G4double* G4KineticTrack::GetActualWidth() const
346 {
347 return theActualWidth;
348 }
349
350 inline void G4KineticTrack::SetActualWidth(G4double* anActualWidth)
351 {
352 theActualWidth = anActualWidth;
353 }
354
355
356
357 inline G4double G4KineticTrack::EvaluateTotalActualWidth()
358 {
359 G4int index;
360 G4double theTotalActualWidth = 0.0;
361 for (index = nChannels - 1; index >= 0; index--)
362 {
363 theTotalActualWidth += theActualWidth[index];
364 }
365 return theTotalActualWidth;
366 }
367
368 inline G4double G4KineticTrack::SampleResidualLifetime()
369 {
370 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
371 G4double tau = CLHEP::hbar_Planck * (-1.0 / theTotalActualWidth);
372 G4double theResidualLifetime = tau * G4Log(G4UniformRand());
373 return theResidualLifetime*the4Momentum.gamma();
374 }
375
376 inline G4double G4KineticTrack::EvaluateCMMomentum(const G4double mass,
377 const G4double* m_ij) const
378 {
379 G4double theCMMomentum;
380 if((m_ij[0]+m_ij[1])<mass)
381 theCMMomentum = 1 / (2 * mass) *
382 std::sqrt (((mass * mass) - (m_ij[0] + m_ij[1]) * (m_ij[0] + m_ij[1])) *
383 ((mass * mass) - (m_ij[0] - m_ij[1]) * (m_ij[0] - m_ij[1])));
384 else
385 theCMMomentum=0.;
386
387 return theCMMomentum;
388 }
389
390 inline G4double G4KineticTrack::BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
391 {
392 G4double Norm = CLHEP::twopi;
393 return (Gamma/((mass-rmass)*(mass-rmass)+Gamma*Gamma/4.))/Norm;
394 }
395
396 inline
397 void G4KineticTrack::Hit()
398 {
399 if(theNucleon)
400 {
401 theNucleon->Hit(1);
402 }
403 }
404
405 inline
406 G4bool G4KineticTrack::IsParticipant() const
407 {
408 if(!theNucleon) return true;
409 return theNucleon->AreYouHit();
410 }
411
412 inline
413 G4KineticTrack::CascadeState G4KineticTrack::GetState() const
414 {
415 return theStateToNucleus;
416 }
417
418 inline
419 G4KineticTrack::CascadeState G4KineticTrack::SetState(const CascadeState new_state)
420 {
421 CascadeState old_state=theStateToNucleus;
422 theStateToNucleus=new_state;
423 return old_state;
424 }
425
426 inline
427 void G4KineticTrack::SetProjectilePotential(G4double aPotential)
428 {
429 theProjectilePotential = aPotential;
430 }
431 inline
432 G4double G4KineticTrack::GetProjectilePotential() const
433 {
434 return theProjectilePotential;
435 }
436
437 inline
438 void G4KineticTrack::SetCreatorModelID(G4int id)
439 {
440 theCreatorModel = id;
441 }
442 inline
443 G4int G4KineticTrack::GetCreatorModelID() const
444 {
445 return theCreatorModel;
446 }
447
448 inline
449 const G4ParticleDefinition* G4KineticTrack::GetParentResonanceDef() const
450 {
451 return theParentResonanceDef;
452 }
453
454 inline
455 void G4KineticTrack::SetParentResonanceDef(const G4ParticleDefinition* parentDef)
456 {
457 theParentResonanceDef = parentDef;
458 }
459
460 inline
461 G4int G4KineticTrack::GetParentResonanceID() const
462 {
463 return theParentResonanceID;
464 }
465
466 inline
467 void G4KineticTrack::SetParentResonanceID(const G4int parentID)
468 {
469 theParentResonanceID = parentID;
470 }
471
472 #endif
473