| 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 // G4ParticleChange class implementation
27 //
28 // Author: Hisaya Kurashige, 23 March 1998
29 // --------------------------------------------------------------------
30
31 #include "G4ParticleChange.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4Track.hh"
35 #include "G4Step.hh"
36 #include "G4DynamicParticle.hh"
37 #include "G4ExceptionSeverity.hh"
38
39 // --------------------------------------------------------------------
40 G4ParticleChange::G4ParticleChange()
41 {}
42
43 // --------------------------------------------------------------------
44 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
45 G4bool IsGoodForTracking)
46 {
47 // create track
48 G4Track* aTrack = new G4Track(aParticle, GetGlobalTime(), thePositionChange);
49
50 // set IsGoodGorTrackingFlag
51 if(IsGoodForTracking)
52 aTrack->SetGoodForTrackingFlag();
53
54 // touchable handle is copied to keep the pointer
55 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
56
57 // add a secondary
58 G4VParticleChange::AddSecondary(aTrack);
59 }
60
61 // --------------------------------------------------------------------
62 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
63 G4ThreeVector newPosition,
64 G4bool IsGoodForTracking)
65 {
66 // create track
67 G4Track* aTrack = new G4Track(aParticle, GetGlobalTime(), newPosition);
68
69 // set IsGoodGorTrackingFlag
70 if(IsGoodForTracking)
71 aTrack->SetGoodForTrackingFlag();
72
73 // touchable is a temporary object, so you cannot keep the pointer
74 aTrack->SetTouchableHandle(nullptr);
75
76 // add a secondary
77 G4VParticleChange::AddSecondary(aTrack);
78 }
79
80 // --------------------------------------------------------------------
81 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
82 G4double newTime, G4bool IsGoodForTracking)
83 {
84 // create track
85 G4Track* aTrack = new G4Track(aParticle, newTime, thePositionChange);
86
87 // set IsGoodGorTrackingFlag
88 if(IsGoodForTracking)
89 aTrack->SetGoodForTrackingFlag();
90
91 // touchable handle is copied to keep the pointer
92 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
93
94 // add a secondary
95 G4VParticleChange::AddSecondary(aTrack);
96 }
97
98 // --------------------------------------------------------------------
99
100 void G4ParticleChange::AddSecondary(G4Track* aTrack)
101 {
102 // add a secondary
103 G4VParticleChange::AddSecondary(aTrack);
104 }
105
106 // --------------------------------------------------------------------
107 void G4ParticleChange::Initialize(const G4Track& track)
108 {
109 // use base class's method at first
110 G4VParticleChange::Initialize(track);
111
112 // set Energy/Momentum etc. equal to those of the parent particle
113 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
114 theEnergyChange = pParticle->GetKineticEnergy();
115 theVelocityChange = track.GetVelocity();
116 isVelocityChanged = false;
117 theMomentumDirectionChange = pParticle->GetMomentumDirection();
118 thePolarizationChange = pParticle->GetPolarization();
119 theProperTimeChange = pParticle->GetProperTime();
120
121 // set mass/charge/MagneticMoment of DynamicParticle
122 theMassChange = pParticle->GetMass();
123 theChargeChange = pParticle->GetCharge();
124 theMagneticMomentChange = pParticle->GetMagneticMoment();
125
126 // set Position equal to those of the parent track
127 thePositionChange = track.GetPosition();
128
129 // set TimeChange equal to local time of the parent track
130 theTimeChange = theLocalTime0 = track.GetLocalTime();
131
132 // set initial global time of the parent track
133 theGlobalTime0 = track.GetGlobalTime();
134 }
135
136 // --------------------------------------------------------------------
137 G4Step* G4ParticleChange::UpdateStepForAlongStep(G4Step* pStep)
138 {
139 // A physics process always calculates the final state of the
140 // particle relative to the initial state at the beginning
141 // of the Step, i.e., based on information of G4Track (or
142 // equivalently the PreStepPoint).
143 // So, the differences (delta) between these two states have to be
144 // calculated and be accumulated in PostStepPoint
145
146 // Take note that the return type of GetMomentumDirectionChange() is a
147 // pointer to G4ParticleMometum. Also it is a normalized
148 // momentum vector
149
150 const G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
151 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
152
153 // set Mass/Charge/MagneticMoment
154 pPostStepPoint->SetMass(theMassChange);
155 pPostStepPoint->SetCharge(theChargeChange);
156 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
157
158 // calculate new kinetic energy
159 G4double preEnergy = pPreStepPoint->GetKineticEnergy();
160 G4double energy =
161 pPostStepPoint->GetKineticEnergy() + (theEnergyChange - preEnergy);
162
163 // update kinetic energy and momentum direction
164 if(energy > 0.0)
165 {
166 // calculate new momentum
167 G4ThreeVector pMomentum = pPostStepPoint->GetMomentum()
168 + (CalcMomentum(theEnergyChange, theMomentumDirectionChange,
169 theMassChange) - pPreStepPoint->GetMomentum());
170 G4double tMomentum2 = pMomentum.mag2();
171 G4ThreeVector direction(1.0, 0.0, 0.0);
172 if(tMomentum2 > 0.)
173 {
174 direction = pMomentum / std::sqrt(tMomentum2);
175 }
176 pPostStepPoint->SetMomentumDirection(direction);
177 pPostStepPoint->SetKineticEnergy(energy);
178
179 // if velocity is not set it should be recomputed
180 //
181 if(!isVelocityChanged)
182 {
183 if (theMassChange > 0.0)
184 {
185 theVelocityChange = CLHEP::c_light *
186 std::sqrt(energy*(energy + 2*theMassChange))/(energy + theMassChange);
187 }
188 else
189 {
190 // zero mass particle
191 theVelocityChange = CLHEP::c_light;
192 // optical photon case
193 if(theCurrentTrack->GetParticleDefinition()->GetPDGEncoding() == -22)
194 {
195 G4Track* pTrack = pStep->GetTrack();
196 G4double e = pTrack->GetKineticEnergy();
197 pTrack->SetKineticEnergy(energy);
198 theVelocityChange = pTrack->CalculateVelocityForOpticalPhoton();
199 pTrack->SetKineticEnergy(e);
200 }
201 }
202 }
203 pPostStepPoint->SetVelocity(theVelocityChange);
204 }
205 else
206 {
207 // stop case
208 pPostStepPoint->SetKineticEnergy(0.0);
209 pPostStepPoint->SetVelocity(0.0);
210 }
211
212 // update polarization
213 pPostStepPoint->AddPolarization(thePolarizationChange -
214 pPreStepPoint->GetPolarization());
215
216 // update position and time
217 pPostStepPoint->AddPosition(thePositionChange - pPreStepPoint->GetPosition());
218 pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
219 pPostStepPoint->AddLocalTime(theTimeChange - theLocalTime0);
220 pPostStepPoint->AddProperTime(theProperTimeChange -
221 pPreStepPoint->GetProperTime());
222
223 if(isParentWeightProposed)
224 {
225 pPostStepPoint->SetWeight(theParentWeight);
226 }
227
228 #ifdef G4VERBOSE
229 if(debugFlag) { CheckIt( *theCurrentTrack ); }
230 #endif
231
232 // update the G4Step specific attributes
233 return UpdateStepInfo(pStep);
234 }
235
236 // --------------------------------------------------------------------
237 G4Step* G4ParticleChange::UpdateStepForPostStep(G4Step* pStep)
238 {
239 // A physics process always calculates the final state of the particle
240
241 // Take note that the return type of GetMomentumChange() is a
242 // pointer to G4ParticleMometum. Also it is a normalized
243 // momentum vector
244
245 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
246 G4Track* pTrack = pStep->GetTrack();
247
248 // set Mass/Charge
249 pPostStepPoint->SetMass(theMassChange);
250 pPostStepPoint->SetCharge(theChargeChange);
251 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
252
253 // update kinetic energy and momentum direction
254 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
255
256 // calculate velocity
257 if(theEnergyChange > 0.0)
258 {
259 pPostStepPoint->SetKineticEnergy(theEnergyChange);
260 pTrack->SetKineticEnergy(theEnergyChange);
261 if(!isVelocityChanged)
262 {
263 theVelocityChange = pTrack->CalculateVelocity();
264 }
265 pPostStepPoint->SetVelocity(theVelocityChange);
266 }
267 else
268 {
269 pPostStepPoint->SetKineticEnergy(0.0);
270 pPostStepPoint->SetVelocity(0.0);
271 }
272
273 // update polarization
274 pPostStepPoint->SetPolarization(thePolarizationChange);
275
276 // update position and time
277 pPostStepPoint->SetPosition(thePositionChange);
278 pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
279 pPostStepPoint->SetLocalTime(theTimeChange);
280 pPostStepPoint->SetProperTime(theProperTimeChange);
281
282 if(isParentWeightProposed)
283 {
284 pPostStepPoint->SetWeight(theParentWeight);
285 }
286
287 #ifdef G4VERBOSE
288 if(debugFlag) { CheckIt( *theCurrentTrack ); }
289 #endif
290
291 // update the G4Step specific attributes
292 return UpdateStepInfo(pStep);
293 }
294
295 // --------------------------------------------------------------------
296 G4Step* G4ParticleChange::UpdateStepForAtRest(G4Step* pStep)
297 {
298 // A physics process always calculates the final state of the particle
299
300 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
301
302 // set Mass/Charge
303 pPostStepPoint->SetMass(theMassChange);
304 pPostStepPoint->SetCharge(theChargeChange);
305 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
306
307 // update kinetic energy and momentum direction
308 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
309 pPostStepPoint->SetKineticEnergy(theEnergyChange);
310 if(!isVelocityChanged)
311 theVelocityChange = theCurrentTrack->CalculateVelocity();
312 pPostStepPoint->SetVelocity(theVelocityChange);
313
314 // update polarization
315 pPostStepPoint->SetPolarization(thePolarizationChange);
316
317 // update position and time
318 pPostStepPoint->SetPosition(thePositionChange);
319 pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
320 pPostStepPoint->SetLocalTime(theTimeChange);
321 pPostStepPoint->SetProperTime(theProperTimeChange);
322
323 if(isParentWeightProposed)
324 {
325 pPostStepPoint->SetWeight(theParentWeight);
326 }
327
328 #ifdef G4VERBOSE
329 if(debugFlag) { CheckIt( *theCurrentTrack ); }
330 #endif
331
332 // update the G4Step specific attributes
333 return UpdateStepInfo(pStep);
334 }
335
336 // --------------------------------------------------------------------
337 void G4ParticleChange::DumpInfo() const
338 {
339 // use base-class DumpInfo
340 G4VParticleChange::DumpInfo();
341
342 G4long oldprc = G4cout.precision(8);
343
344 G4cout << " Mass (GeV) : " << std::setw(20) << theMassChange / GeV
345 << G4endl;
346 G4cout << " Charge (eplus) : " << std::setw(20)
347 << theChargeChange / eplus << G4endl;
348 G4cout << " MagneticMoment : " << std::setw(20)
349 << theMagneticMomentChange << G4endl;
350 G4cout << " = : " << std::setw(20)
351 << theMagneticMomentChange
352 * 2. * theMassChange / c_squared / eplus / hbar_Planck
353 << "*[e hbar]/[2 m]" << G4endl;
354 G4cout << " Position - x (mm) : " << std::setw(20)
355 << thePositionChange.x() / mm << G4endl;
356 G4cout << " Position - y (mm) : " << std::setw(20)
357 << thePositionChange.y() / mm << G4endl;
358 G4cout << " Position - z (mm) : " << std::setw(20)
359 << thePositionChange.z() / mm << G4endl;
360 G4cout << " Time (ns) : " << std::setw(20)
361 << theTimeChange / ns << G4endl;
362 G4cout << " Proper Time (ns) : " << std::setw(20)
363 << theProperTimeChange / ns << G4endl;
364 G4cout << " Momentum Direct - x : " << std::setw(20)
365 << theMomentumDirectionChange.x() << G4endl;
366 G4cout << " Momentum Direct - y : " << std::setw(20)
367 << theMomentumDirectionChange.y() << G4endl;
368 G4cout << " Momentum Direct - z : " << std::setw(20)
369 << theMomentumDirectionChange.z() << G4endl;
370 G4cout << " Kinetic Energy (MeV): " << std::setw(20)
371 << theEnergyChange / MeV << G4endl;
372 G4cout << " Velocity (/c) : " << std::setw(20)
373 << theVelocityChange / c_light << G4endl;
374 G4cout << " Polarization - x : " << std::setw(20)
375 << thePolarizationChange.x() << G4endl;
376 G4cout << " Polarization - y : " << std::setw(20)
377 << thePolarizationChange.y() << G4endl;
378 G4cout << " Polarization - z : " << std::setw(20)
379 << thePolarizationChange.z() << G4endl;
380 G4cout.precision(oldprc);
381 }
382