| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4ParticleChange class implementation 26 // G4ParticleChange class implementation
27 // 27 //
28 // Author: Hisaya Kurashige, 23 March 1998 28 // Author: Hisaya Kurashige, 23 March 1998
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4ParticleChange.hh" 31 #include "G4ParticleChange.hh"
32 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
34 #include "G4Track.hh" 34 #include "G4Track.hh"
35 #include "G4Step.hh" 35 #include "G4Step.hh"
>> 36 #include "G4TrackFastVector.hh"
36 #include "G4DynamicParticle.hh" 37 #include "G4DynamicParticle.hh"
37 #include "G4ExceptionSeverity.hh" 38 #include "G4ExceptionSeverity.hh"
38 39
39 // ------------------------------------------- 40 // --------------------------------------------------------------------
40 G4ParticleChange::G4ParticleChange() 41 G4ParticleChange::G4ParticleChange()
41 {} << 42 : G4VParticleChange()
>> 43 {
>> 44 }
>> 45
>> 46 // --------------------------------------------------------------------
>> 47 G4ParticleChange::~G4ParticleChange()
>> 48 {
>> 49 }
>> 50
>> 51 // --------------------------------------------------------------------
>> 52 G4ParticleChange::G4ParticleChange(const G4ParticleChange& right)
>> 53 : G4VParticleChange(right)
>> 54 {
>> 55 theCurrentTrack = right.theCurrentTrack;
>> 56
>> 57 theMomentumDirectionChange = right.theMomentumDirectionChange;
>> 58 thePolarizationChange = right.thePolarizationChange;
>> 59 thePositionChange = right.thePositionChange;
>> 60 theGlobalTime0 = right.theGlobalTime0;
>> 61 theLocalTime0 = right.theLocalTime0;
>> 62 theTimeChange = right.theTimeChange;
>> 63 theProperTimeChange = right.theProperTimeChange;
>> 64 theEnergyChange = right.theEnergyChange;
>> 65 theVelocityChange = right.theVelocityChange;
>> 66 isVelocityChanged = true;
>> 67 theMassChange = right.theMassChange;
>> 68 theChargeChange = right.theChargeChange;
>> 69 theMagneticMomentChange = right.theMagneticMomentChange;
>> 70 }
>> 71
>> 72 // --------------------------------------------------------------------
>> 73 G4ParticleChange& G4ParticleChange::operator=(const G4ParticleChange& right)
>> 74 {
>> 75 if(this != &right)
>> 76 {
>> 77 if(theNumberOfSecondaries > 0)
>> 78 {
>> 79 for(G4int index = 0; index < theNumberOfSecondaries; ++index)
>> 80 {
>> 81 if((*theListOfSecondaries)[index])
>> 82 delete(*theListOfSecondaries)[index];
>> 83 }
>> 84 }
>> 85 delete theListOfSecondaries;
>> 86
>> 87 theListOfSecondaries = new G4TrackFastVector();
>> 88 theNumberOfSecondaries = right.theNumberOfSecondaries;
>> 89 for(G4int index = 0; index < theNumberOfSecondaries; ++index)
>> 90 {
>> 91 G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index]));
>> 92 theListOfSecondaries->SetElement(index, newTrack);
>> 93 }
>> 94
>> 95 theStatusChange = right.theStatusChange;
>> 96 theCurrentTrack = right.theCurrentTrack;
>> 97
>> 98 theMomentumDirectionChange = right.theMomentumDirectionChange;
>> 99 thePolarizationChange = right.thePolarizationChange;
>> 100 thePositionChange = right.thePositionChange;
>> 101 theGlobalTime0 = right.theGlobalTime0;
>> 102 theLocalTime0 = right.theLocalTime0;
>> 103 theTimeChange = right.theTimeChange;
>> 104 theProperTimeChange = right.theProperTimeChange;
>> 105 theEnergyChange = right.theEnergyChange;
>> 106 theVelocityChange = right.theVelocityChange;
>> 107 isVelocityChanged = true;
>> 108 theMassChange = right.theMassChange;
>> 109 theChargeChange = right.theChargeChange;
>> 110 theMagneticMomentChange = right.theMagneticMomentChange;
>> 111
>> 112 theTrueStepLength = right.theTrueStepLength;
>> 113 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
>> 114 theSteppingControlFlag = right.theSteppingControlFlag;
>> 115 }
>> 116 return *this;
>> 117 }
>> 118
>> 119 // --------------------------------------------------------------------
>> 120 G4bool G4ParticleChange::operator==(const G4ParticleChange& right) const
>> 121 {
>> 122 return ((G4VParticleChange*) this == (G4VParticleChange*) &right);
>> 123 }
>> 124
>> 125 // --------------------------------------------------------------------
>> 126 G4bool G4ParticleChange::operator!=(const G4ParticleChange& right) const
>> 127 {
>> 128 return ((G4VParticleChange*) this != (G4VParticleChange*) &right);
>> 129 }
42 130
43 // ------------------------------------------- 131 // --------------------------------------------------------------------
44 void G4ParticleChange::AddSecondary(G4DynamicP 132 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
45 G4bool IsG 133 G4bool IsGoodForTracking)
46 { 134 {
47 // create track 135 // create track
48 G4Track* aTrack = new G4Track(aParticle, Get 136 G4Track* aTrack = new G4Track(aParticle, GetGlobalTime(), thePositionChange);
49 137
50 // set IsGoodGorTrackingFlag 138 // set IsGoodGorTrackingFlag
51 if(IsGoodForTracking) 139 if(IsGoodForTracking)
52 aTrack->SetGoodForTrackingFlag(); 140 aTrack->SetGoodForTrackingFlag();
53 141
54 // touchable handle is copied to keep the po 142 // touchable handle is copied to keep the pointer
55 aTrack->SetTouchableHandle(theCurrentTrack-> 143 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
56 144
57 // add a secondary 145 // add a secondary
58 G4VParticleChange::AddSecondary(aTrack); 146 G4VParticleChange::AddSecondary(aTrack);
59 } 147 }
60 148
61 // ------------------------------------------- 149 // --------------------------------------------------------------------
62 void G4ParticleChange::AddSecondary(G4DynamicP 150 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
63 G4ThreeVec 151 G4ThreeVector newPosition,
64 G4bool IsG 152 G4bool IsGoodForTracking)
65 { 153 {
66 // create track 154 // create track
67 G4Track* aTrack = new G4Track(aParticle, Get 155 G4Track* aTrack = new G4Track(aParticle, GetGlobalTime(), newPosition);
68 156
69 // set IsGoodGorTrackingFlag 157 // set IsGoodGorTrackingFlag
70 if(IsGoodForTracking) 158 if(IsGoodForTracking)
71 aTrack->SetGoodForTrackingFlag(); 159 aTrack->SetGoodForTrackingFlag();
72 160
73 // touchable is a temporary object, so you c 161 // touchable is a temporary object, so you cannot keep the pointer
74 aTrack->SetTouchableHandle(nullptr); 162 aTrack->SetTouchableHandle(nullptr);
75 163
76 // add a secondary 164 // add a secondary
77 G4VParticleChange::AddSecondary(aTrack); 165 G4VParticleChange::AddSecondary(aTrack);
78 } 166 }
79 167
80 // ------------------------------------------- 168 // --------------------------------------------------------------------
81 void G4ParticleChange::AddSecondary(G4DynamicP 169 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
82 G4double n 170 G4double newTime, G4bool IsGoodForTracking)
83 { 171 {
84 // create track 172 // create track
85 G4Track* aTrack = new G4Track(aParticle, new 173 G4Track* aTrack = new G4Track(aParticle, newTime, thePositionChange);
86 174
87 // set IsGoodGorTrackingFlag 175 // set IsGoodGorTrackingFlag
88 if(IsGoodForTracking) 176 if(IsGoodForTracking)
89 aTrack->SetGoodForTrackingFlag(); 177 aTrack->SetGoodForTrackingFlag();
90 178
91 // touchable handle is copied to keep the po 179 // touchable handle is copied to keep the pointer
92 aTrack->SetTouchableHandle(theCurrentTrack-> 180 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
93 181
94 // add a secondary 182 // add a secondary
95 G4VParticleChange::AddSecondary(aTrack); 183 G4VParticleChange::AddSecondary(aTrack);
96 } 184 }
97 185
98 // ------------------------------------------- 186 // --------------------------------------------------------------------
99 <<
100 void G4ParticleChange::AddSecondary(G4Track* a 187 void G4ParticleChange::AddSecondary(G4Track* aTrack)
101 { 188 {
102 // add a secondary 189 // add a secondary
103 G4VParticleChange::AddSecondary(aTrack); 190 G4VParticleChange::AddSecondary(aTrack);
104 } 191 }
105 192
106 // ------------------------------------------- 193 // --------------------------------------------------------------------
107 void G4ParticleChange::Initialize(const G4Trac 194 void G4ParticleChange::Initialize(const G4Track& track)
108 { 195 {
109 // use base class's method at first 196 // use base class's method at first
110 G4VParticleChange::Initialize(track); 197 G4VParticleChange::Initialize(track);
>> 198 theCurrentTrack = &track;
111 199
112 // set Energy/Momentum etc. equal to those o 200 // set Energy/Momentum etc. equal to those of the parent particle
113 const G4DynamicParticle* pParticle = track.G 201 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
114 theEnergyChange = pPartic 202 theEnergyChange = pParticle->GetKineticEnergy();
115 theVelocityChange = track.G 203 theVelocityChange = track.GetVelocity();
116 isVelocityChanged = false; 204 isVelocityChanged = false;
117 theMomentumDirectionChange = pPartic 205 theMomentumDirectionChange = pParticle->GetMomentumDirection();
118 thePolarizationChange = pPartic 206 thePolarizationChange = pParticle->GetPolarization();
119 theProperTimeChange = pPartic 207 theProperTimeChange = pParticle->GetProperTime();
120 208
121 // set mass/charge/MagneticMoment of Dynamic 209 // set mass/charge/MagneticMoment of DynamicParticle
122 theMassChange = pParticle->GetMass 210 theMassChange = pParticle->GetMass();
123 theChargeChange = pParticle->GetChar 211 theChargeChange = pParticle->GetCharge();
124 theMagneticMomentChange = pParticle->GetMagn 212 theMagneticMomentChange = pParticle->GetMagneticMoment();
125 213
126 // set Position equal to those of the parent 214 // set Position equal to those of the parent track
127 thePositionChange = track.GetPosition(); 215 thePositionChange = track.GetPosition();
128 216
129 // set TimeChange equal to local time of the 217 // set TimeChange equal to local time of the parent track
130 theTimeChange = theLocalTime0 = track.GetLoc << 218 theTimeChange = track.GetLocalTime();
131 219
132 // set initial global time of the parent tra << 220 // set initial Local/Global time of the parent track
>> 221 theLocalTime0 = track.GetLocalTime();
133 theGlobalTime0 = track.GetGlobalTime(); 222 theGlobalTime0 = track.GetGlobalTime();
134 } 223 }
135 224
136 // ------------------------------------------- 225 // --------------------------------------------------------------------
137 G4Step* G4ParticleChange::UpdateStepForAlongSt 226 G4Step* G4ParticleChange::UpdateStepForAlongStep(G4Step* pStep)
138 { 227 {
139 // A physics process always calculates the f 228 // A physics process always calculates the final state of the
140 // particle relative to the initial state at 229 // particle relative to the initial state at the beginning
141 // of the Step, i.e., based on information o 230 // of the Step, i.e., based on information of G4Track (or
142 // equivalently the PreStepPoint). 231 // equivalently the PreStepPoint).
143 // So, the differences (delta) between these 232 // So, the differences (delta) between these two states have to be
144 // calculated and be accumulated in PostStep 233 // calculated and be accumulated in PostStepPoint
145 234
146 // Take note that the return type of GetMome 235 // Take note that the return type of GetMomentumDirectionChange() is a
147 // pointer to G4ParticleMometum. Also it is 236 // pointer to G4ParticleMometum. Also it is a normalized
148 // momentum vector 237 // momentum vector
149 238
150 const G4StepPoint* pPreStepPoint = pStep->Ge << 239 G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
151 G4StepPoint* pPostStepPoint = pStep->GetPost 240 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
>> 241 G4Track* pTrack = pStep->GetTrack();
>> 242 G4double mass = theMassChange;
152 243
153 // set Mass/Charge/MagneticMoment 244 // set Mass/Charge/MagneticMoment
154 pPostStepPoint->SetMass(theMassChange); 245 pPostStepPoint->SetMass(theMassChange);
155 pPostStepPoint->SetCharge(theChargeChange); 246 pPostStepPoint->SetCharge(theChargeChange);
156 pPostStepPoint->SetMagneticMoment(theMagneti 247 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
157 248
158 // calculate new kinetic energy 249 // calculate new kinetic energy
159 G4double preEnergy = pPreStepPoint->GetKinet 250 G4double preEnergy = pPreStepPoint->GetKineticEnergy();
160 G4double energy = 251 G4double energy =
161 pPostStepPoint->GetKineticEnergy() + (theE 252 pPostStepPoint->GetKineticEnergy() + (theEnergyChange - preEnergy);
162 253
163 // update kinetic energy and momentum direct 254 // update kinetic energy and momentum direction
164 if(energy > 0.0) 255 if(energy > 0.0)
165 { 256 {
166 // calculate new momentum 257 // calculate new momentum
167 G4ThreeVector pMomentum = pPostStepPoint-> 258 G4ThreeVector pMomentum = pPostStepPoint->GetMomentum()
168 + (CalcMomentum(theEnergyChange, theMome << 259 + (CalcMomentum(theEnergyChange, theMomentumDirectionChange, mass)
169 theMassChange) - pPreStepPoint->GetM << 260 - pPreStepPoint->GetMomentum());
170 G4double tMomentum2 = pMomentum.mag2(); << 261 G4double tMomentum = pMomentum.mag();
171 G4ThreeVector direction(1.0, 0.0, 0.0); 262 G4ThreeVector direction(1.0, 0.0, 0.0);
172 if(tMomentum2 > 0.) << 263 if(tMomentum > 0.)
173 { 264 {
174 direction = pMomentum / std::sqrt(tMomen << 265 G4double inv_Momentum = 1.0 / tMomentum;
>> 266 direction = pMomentum * inv_Momentum;
175 } 267 }
176 pPostStepPoint->SetMomentumDirection(direc 268 pPostStepPoint->SetMomentumDirection(direction);
177 pPostStepPoint->SetKineticEnergy(energy); 269 pPostStepPoint->SetKineticEnergy(energy);
178 <<
179 // if velocity is not set it should be rec <<
180 // <<
181 if(!isVelocityChanged) <<
182 { <<
183 if (theMassChange > 0.0) <<
184 { <<
185 theVelocityChange = CLHEP::c_light * <<
186 std::sqrt(energy*(energy + 2*theMassChange <<
187 } <<
188 else <<
189 { <<
190 // zero mass particle <<
191 theVelocityChange = CLHEP::c_light; <<
192 // optical photon case <<
193 if(theCurrentTrack->GetParticleDefinition()- <<
194 { <<
195 G4Track* pTrack = pStep->GetTrack(); <<
196 G4double e = pTrack->GetKineticEnerg <<
197 pTrack->SetKineticEnergy(energy); <<
198 theVelocityChange = pTrack->Calculat <<
199 pTrack->SetKineticEnergy(e); <<
200 } <<
201 } <<
202 } <<
203 pPostStepPoint->SetVelocity(theVelocityCha <<
204 } 270 }
205 else 271 else
206 { 272 {
207 // stop case 273 // stop case
>> 274 // pPostStepPoint->SetMomentumDirection(G4ThreeVector(1., 0., 0.));
208 pPostStepPoint->SetKineticEnergy(0.0); 275 pPostStepPoint->SetKineticEnergy(0.0);
209 pPostStepPoint->SetVelocity(0.0); <<
210 } 276 }
>> 277 // calculate velocity
>> 278 if(!isVelocityChanged)
>> 279 {
>> 280 if(energy > 0.0)
>> 281 {
>> 282 pTrack->SetKineticEnergy(energy);
>> 283 theVelocityChange = pTrack->CalculateVelocity();
>> 284 pTrack->SetKineticEnergy(preEnergy);
>> 285 }
>> 286 else if(theMassChange > 0.0)
>> 287 {
>> 288 theVelocityChange = 0.0;
>> 289 }
>> 290 }
>> 291 pPostStepPoint->SetVelocity(theVelocityChange);
211 292
212 // update polarization 293 // update polarization
213 pPostStepPoint->AddPolarization(thePolarizat 294 pPostStepPoint->AddPolarization(thePolarizationChange -
214 pPreStepPoin 295 pPreStepPoint->GetPolarization());
215 296
216 // update position and time 297 // update position and time
217 pPostStepPoint->AddPosition(thePositionChang 298 pPostStepPoint->AddPosition(thePositionChange - pPreStepPoint->GetPosition());
218 pPostStepPoint->AddGlobalTime(theTimeChange 299 pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
219 pPostStepPoint->AddLocalTime(theTimeChange - 300 pPostStepPoint->AddLocalTime(theTimeChange - theLocalTime0);
220 pPostStepPoint->AddProperTime(theProperTimeC 301 pPostStepPoint->AddProperTime(theProperTimeChange -
221 pPreStepPoint- 302 pPreStepPoint->GetProperTime());
222 303
223 if(isParentWeightProposed) 304 if(isParentWeightProposed)
224 { 305 {
225 pPostStepPoint->SetWeight(theParentWeight) 306 pPostStepPoint->SetWeight(theParentWeight);
226 } 307 }
227 308
228 #ifdef G4VERBOSE 309 #ifdef G4VERBOSE
229 if(debugFlag) { CheckIt( *theCurrentTrack ); << 310 G4Track* aTrack = pStep->GetTrack();
>> 311 if(debugFlag) { CheckIt(*aTrack); }
230 #endif 312 #endif
231 313
232 // update the G4Step specific attributes 314 // update the G4Step specific attributes
233 return UpdateStepInfo(pStep); 315 return UpdateStepInfo(pStep);
234 } 316 }
235 317
236 // ------------------------------------------- 318 // --------------------------------------------------------------------
237 G4Step* G4ParticleChange::UpdateStepForPostSte 319 G4Step* G4ParticleChange::UpdateStepForPostStep(G4Step* pStep)
238 { 320 {
239 // A physics process always calculates the f 321 // A physics process always calculates the final state of the particle
240 322
241 // Take note that the return type of GetMome 323 // Take note that the return type of GetMomentumChange() is a
242 // pointer to G4ParticleMometum. Also it is 324 // pointer to G4ParticleMometum. Also it is a normalized
243 // momentum vector 325 // momentum vector
244 326
245 G4StepPoint* pPostStepPoint = pStep->GetPost 327 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
246 G4Track* pTrack = pStep->GetTrac 328 G4Track* pTrack = pStep->GetTrack();
247 329
248 // set Mass/Charge 330 // set Mass/Charge
249 pPostStepPoint->SetMass(theMassChange); 331 pPostStepPoint->SetMass(theMassChange);
250 pPostStepPoint->SetCharge(theChargeChange); 332 pPostStepPoint->SetCharge(theChargeChange);
251 pPostStepPoint->SetMagneticMoment(theMagneti 333 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
252 334
253 // update kinetic energy and momentum direct 335 // update kinetic energy and momentum direction
254 pPostStepPoint->SetMomentumDirection(theMome 336 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
>> 337 pPostStepPoint->SetKineticEnergy(theEnergyChange);
255 338
256 // calculate velocity 339 // calculate velocity
257 if(theEnergyChange > 0.0) << 340 pTrack->SetKineticEnergy(theEnergyChange);
>> 341 if(!isVelocityChanged)
258 { 342 {
259 pPostStepPoint->SetKineticEnergy(theEnergy << 343 if(theEnergyChange > 0.0)
260 pTrack->SetKineticEnergy(theEnergyChange); <<
261 if(!isVelocityChanged) <<
262 { 344 {
263 theVelocityChange = pTrack->CalculateVel 345 theVelocityChange = pTrack->CalculateVelocity();
264 } 346 }
265 pPostStepPoint->SetVelocity(theVelocityCha << 347 else if(theMassChange > 0.0)
266 } << 348 {
267 else << 349 theVelocityChange = 0.0;
268 { << 350 }
269 pPostStepPoint->SetKineticEnergy(0.0); <<
270 pPostStepPoint->SetVelocity(0.0); <<
271 } 351 }
>> 352 pPostStepPoint->SetVelocity(theVelocityChange);
272 353
273 // update polarization 354 // update polarization
274 pPostStepPoint->SetPolarization(thePolarizat 355 pPostStepPoint->SetPolarization(thePolarizationChange);
275 356
276 // update position and time 357 // update position and time
277 pPostStepPoint->SetPosition(thePositionChang 358 pPostStepPoint->SetPosition(thePositionChange);
278 pPostStepPoint->AddGlobalTime(theTimeChange 359 pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
279 pPostStepPoint->SetLocalTime(theTimeChange); 360 pPostStepPoint->SetLocalTime(theTimeChange);
280 pPostStepPoint->SetProperTime(theProperTimeC 361 pPostStepPoint->SetProperTime(theProperTimeChange);
281 362
282 if(isParentWeightProposed) 363 if(isParentWeightProposed)
283 { 364 {
284 pPostStepPoint->SetWeight(theParentWeight) 365 pPostStepPoint->SetWeight(theParentWeight);
285 } 366 }
286 367
287 #ifdef G4VERBOSE 368 #ifdef G4VERBOSE
288 if(debugFlag) { CheckIt( *theCurrentTrack ); << 369 G4Track* aTrack = pStep->GetTrack();
>> 370 if(debugFlag) { CheckIt(*aTrack); }
289 #endif 371 #endif
290 372
291 // update the G4Step specific attributes 373 // update the G4Step specific attributes
292 return UpdateStepInfo(pStep); 374 return UpdateStepInfo(pStep);
293 } 375 }
294 376
295 // ------------------------------------------- 377 // --------------------------------------------------------------------
296 G4Step* G4ParticleChange::UpdateStepForAtRest( 378 G4Step* G4ParticleChange::UpdateStepForAtRest(G4Step* pStep)
297 { 379 {
298 // A physics process always calculates the f 380 // A physics process always calculates the final state of the particle
299 381
300 G4StepPoint* pPostStepPoint = pStep->GetPost 382 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
301 383
302 // set Mass/Charge 384 // set Mass/Charge
303 pPostStepPoint->SetMass(theMassChange); 385 pPostStepPoint->SetMass(theMassChange);
304 pPostStepPoint->SetCharge(theChargeChange); 386 pPostStepPoint->SetCharge(theChargeChange);
305 pPostStepPoint->SetMagneticMoment(theMagneti 387 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
306 388
307 // update kinetic energy and momentum direct 389 // update kinetic energy and momentum direction
308 pPostStepPoint->SetMomentumDirection(theMome 390 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
309 pPostStepPoint->SetKineticEnergy(theEnergyCh 391 pPostStepPoint->SetKineticEnergy(theEnergyChange);
310 if(!isVelocityChanged) 392 if(!isVelocityChanged)
311 theVelocityChange = theCurrentTrack->Calcu << 393 theVelocityChange = pStep->GetTrack()->CalculateVelocity();
312 pPostStepPoint->SetVelocity(theVelocityChang 394 pPostStepPoint->SetVelocity(theVelocityChange);
313 395
314 // update polarization 396 // update polarization
315 pPostStepPoint->SetPolarization(thePolarizat 397 pPostStepPoint->SetPolarization(thePolarizationChange);
316 398
317 // update position and time 399 // update position and time
318 pPostStepPoint->SetPosition(thePositionChang 400 pPostStepPoint->SetPosition(thePositionChange);
319 pPostStepPoint->AddGlobalTime(theTimeChange 401 pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
320 pPostStepPoint->SetLocalTime(theTimeChange); 402 pPostStepPoint->SetLocalTime(theTimeChange);
321 pPostStepPoint->SetProperTime(theProperTimeC 403 pPostStepPoint->SetProperTime(theProperTimeChange);
322 404
323 if(isParentWeightProposed) 405 if(isParentWeightProposed)
324 { 406 {
325 pPostStepPoint->SetWeight(theParentWeight) 407 pPostStepPoint->SetWeight(theParentWeight);
326 } 408 }
327 409
328 #ifdef G4VERBOSE 410 #ifdef G4VERBOSE
329 if(debugFlag) { CheckIt( *theCurrentTrack ); << 411 G4Track* aTrack = pStep->GetTrack();
>> 412 if(debugFlag) { CheckIt(*aTrack); }
330 #endif 413 #endif
331 414
332 // update the G4Step specific attributes 415 // update the G4Step specific attributes
333 return UpdateStepInfo(pStep); 416 return UpdateStepInfo(pStep);
334 } 417 }
335 418
336 // ------------------------------------------- 419 // --------------------------------------------------------------------
337 void G4ParticleChange::DumpInfo() const 420 void G4ParticleChange::DumpInfo() const
338 { 421 {
339 // use base-class DumpInfo 422 // use base-class DumpInfo
340 G4VParticleChange::DumpInfo(); 423 G4VParticleChange::DumpInfo();
341 424
342 G4long oldprc = G4cout.precision(8); << 425 G4int oldprc = G4cout.precision(3);
343 426
344 G4cout << " Mass (GeV) : " < << 427 G4cout << " Mass (GeV) : " << std::setw(20) << theMassChange / GeV
345 << G4endl; 428 << G4endl;
346 G4cout << " Charge (eplus) : " < << 429 G4cout << " Charge (eplus) : " << std::setw(20)
347 << theChargeChange / eplus << G4endl; 430 << theChargeChange / eplus << G4endl;
348 G4cout << " MagneticMoment : " < << 431 G4cout << " MagneticMoment : " << std::setw(20)
349 << theMagneticMomentChange << G4endl; 432 << theMagneticMomentChange << G4endl;
350 G4cout << " = : " < << 433 G4cout << " : = " << std::setw(20)
351 << theMagneticMomentChange 434 << theMagneticMomentChange
352 * 2. * theMassChange / c_squared / 435 * 2. * theMassChange / c_squared / eplus / hbar_Planck
353 << "*[e hbar]/[2 m]" << G4endl; 436 << "*[e hbar]/[2 m]" << G4endl;
354 G4cout << " Position - x (mm) : " < 437 G4cout << " Position - x (mm) : " << std::setw(20)
355 << thePositionChange.x() / mm << G4en 438 << thePositionChange.x() / mm << G4endl;
356 G4cout << " Position - y (mm) : " < 439 G4cout << " Position - y (mm) : " << std::setw(20)
357 << thePositionChange.y() / mm << G4en 440 << thePositionChange.y() / mm << G4endl;
358 G4cout << " Position - z (mm) : " < 441 G4cout << " Position - z (mm) : " << std::setw(20)
359 << thePositionChange.z() / mm << G4en 442 << thePositionChange.z() / mm << G4endl;
360 G4cout << " Time (ns) : " < 443 G4cout << " Time (ns) : " << std::setw(20)
361 << theTimeChange / ns << G4endl; 444 << theTimeChange / ns << G4endl;
362 G4cout << " Proper Time (ns) : " < 445 G4cout << " Proper Time (ns) : " << std::setw(20)
363 << theProperTimeChange / ns << G4endl 446 << theProperTimeChange / ns << G4endl;
364 G4cout << " Momentum Direct - x : " < 447 G4cout << " Momentum Direct - x : " << std::setw(20)
365 << theMomentumDirectionChange.x() << 448 << theMomentumDirectionChange.x() << G4endl;
366 G4cout << " Momentum Direct - y : " < 449 G4cout << " Momentum Direct - y : " << std::setw(20)
367 << theMomentumDirectionChange.y() << 450 << theMomentumDirectionChange.y() << G4endl;
368 G4cout << " Momentum Direct - z : " < 451 G4cout << " Momentum Direct - z : " << std::setw(20)
369 << theMomentumDirectionChange.z() << 452 << theMomentumDirectionChange.z() << G4endl;
370 G4cout << " Kinetic Energy (MeV): " < 453 G4cout << " Kinetic Energy (MeV): " << std::setw(20)
371 << theEnergyChange / MeV << G4endl; 454 << theEnergyChange / MeV << G4endl;
372 G4cout << " Velocity (/c) : " < << 455 G4cout << " Velocity (/c): " << std::setw(20)
373 << theVelocityChange / c_light << G4e 456 << theVelocityChange / c_light << G4endl;
374 G4cout << " Polarization - x : " < 457 G4cout << " Polarization - x : " << std::setw(20)
375 << thePolarizationChange.x() << G4end 458 << thePolarizationChange.x() << G4endl;
376 G4cout << " Polarization - y : " < 459 G4cout << " Polarization - y : " << std::setw(20)
377 << thePolarizationChange.y() << G4end 460 << thePolarizationChange.y() << G4endl;
378 G4cout << " Polarization - z : " < 461 G4cout << " Polarization - z : " << std::setw(20)
379 << thePolarizationChange.z() << G4end 462 << thePolarizationChange.z() << G4endl;
380 G4cout.precision(oldprc); 463 G4cout.precision(oldprc);
>> 464 }
>> 465
>> 466 // --------------------------------------------------------------------
>> 467 G4bool G4ParticleChange::CheckIt(const G4Track& aTrack)
>> 468 {
>> 469 G4bool exitWithError = false;
>> 470 G4double accuracy;
>> 471 static G4ThreadLocal G4int nError = 0;
>> 472 #ifdef G4VERBOSE
>> 473 const G4int maxError = 30;
>> 474 #endif
>> 475
>> 476 // no check in case of "fStopAndKill"
>> 477 if(GetTrackStatus() == fStopAndKill)
>> 478 return G4VParticleChange::CheckIt(aTrack);
>> 479
>> 480 // MomentumDirection should be unit vector
>> 481 G4bool itsOKforMomentum = true;
>> 482 if(theEnergyChange > 0.)
>> 483 {
>> 484 accuracy = std::fabs(theMomentumDirectionChange.mag2() - 1.0);
>> 485 if(accuracy > accuracyForWarning)
>> 486 {
>> 487 itsOKforMomentum = false;
>> 488 nError += 1;
>> 489 exitWithError = exitWithError || (accuracy > accuracyForException);
>> 490 #ifdef G4VERBOSE
>> 491 if(nError < maxError)
>> 492 {
>> 493 G4cout << " G4ParticleChange::CheckIt : ";
>> 494 G4cout << "the Momentum Change is not unit vector !!"
>> 495 << " Difference: " << accuracy << G4endl;
>> 496 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 497 << " E=" << aTrack.GetKineticEnergy() / MeV
>> 498 << " pos=" << aTrack.GetPosition().x() / m << ", "
>> 499 << aTrack.GetPosition().y() / m << ", "
>> 500 << aTrack.GetPosition().z() / m << G4endl;
>> 501 }
>> 502 #endif
>> 503 }
>> 504 }
>> 505
>> 506 // both global and proper time should not go back
>> 507 G4bool itsOKforGlobalTime = true;
>> 508 accuracy = (aTrack.GetLocalTime() - theTimeChange) / ns;
>> 509 if(accuracy > accuracyForWarning)
>> 510 {
>> 511 itsOKforGlobalTime = false;
>> 512 nError += 1;
>> 513 exitWithError = exitWithError || (accuracy > accuracyForException);
>> 514 #ifdef G4VERBOSE
>> 515 if(nError < maxError)
>> 516 {
>> 517 G4cout << " G4ParticleChange::CheckIt : ";
>> 518 G4cout << "the local time goes back !!"
>> 519 << " Difference: " << accuracy << "[ns] " << G4endl;
>> 520 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 521 << " E=" << aTrack.GetKineticEnergy() / MeV
>> 522 << " pos=" << aTrack.GetPosition().x() / m << ", "
>> 523 << aTrack.GetPosition().y() / m << ", "
>> 524 << aTrack.GetPosition().z() / m
>> 525 << " global time=" << aTrack.GetGlobalTime() / ns
>> 526 << " local time=" << aTrack.GetLocalTime() / ns
>> 527 << " proper time=" << aTrack.GetProperTime() / ns << G4endl;
>> 528 }
>> 529 #endif
>> 530 }
>> 531
>> 532 G4bool itsOKforProperTime = true;
>> 533 accuracy = (aTrack.GetProperTime() - theProperTimeChange) / ns;
>> 534 if(accuracy > accuracyForWarning)
>> 535 {
>> 536 itsOKforProperTime = false;
>> 537 nError += 1;
>> 538 exitWithError = exitWithError || (accuracy > accuracyForException);
>> 539 #ifdef G4VERBOSE
>> 540 if(nError < maxError)
>> 541 {
>> 542 G4cout << " G4ParticleChange::CheckIt : ";
>> 543 G4cout << "the proper time goes back !!"
>> 544 << " Difference: " << accuracy << "[ns] " << G4endl;
>> 545 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 546 << " E=" << aTrack.GetKineticEnergy() / MeV
>> 547 << " pos=" << aTrack.GetPosition().x() / m << ", "
>> 548 << aTrack.GetPosition().y() / m << ", "
>> 549 << aTrack.GetPosition().z() / m
>> 550 << " global time=" << aTrack.GetGlobalTime() / ns
>> 551 << " local time=" << aTrack.GetLocalTime() / ns
>> 552 << " proper time=" << aTrack.GetProperTime() / ns << G4endl;
>> 553 }
>> 554 #endif
>> 555 }
>> 556
>> 557 // kinetic Energy should not be negative
>> 558 G4bool itsOKforEnergy = true;
>> 559 accuracy = -1.0 * theEnergyChange / MeV;
>> 560 if(accuracy > accuracyForWarning)
>> 561 {
>> 562 itsOKforEnergy = false;
>> 563 nError += 1;
>> 564 exitWithError = exitWithError || (accuracy > accuracyForException);
>> 565 #ifdef G4VERBOSE
>> 566 if(nError < maxError)
>> 567 {
>> 568 G4cout << " G4ParticleChange::CheckIt : ";
>> 569 G4cout << "the kinetic energy is negative !!"
>> 570 << " Difference: " << accuracy << "[MeV] " << G4endl;
>> 571 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 572 << " E=" << aTrack.GetKineticEnergy() / MeV
>> 573 << " pos=" << aTrack.GetPosition().x() / m << ", "
>> 574 << aTrack.GetPosition().y() / m << ", "
>> 575 << aTrack.GetPosition().z() / m << G4endl;
>> 576 }
>> 577 #endif
>> 578 }
>> 579
>> 580 // velocity should not be less than c_light
>> 581 G4bool itsOKforVelocity = true;
>> 582 if(theVelocityChange < 0.)
>> 583 {
>> 584 itsOKforVelocity = false;
>> 585 nError += 1;
>> 586 exitWithError = true;
>> 587 #ifdef G4VERBOSE
>> 588 if(nError < maxError)
>> 589 {
>> 590 G4cout << " G4ParticleChange::CheckIt : ";
>> 591 G4cout << "the velocity is negative !!"
>> 592 << " Velocity: " << theVelocityChange / c_light << G4endl;
>> 593 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 594 << " E=" << aTrack.GetKineticEnergy() / MeV
>> 595 << " pos=" << aTrack.GetPosition().x() / m << ", "
>> 596 << aTrack.GetPosition().y() / m << ", "
>> 597 << aTrack.GetPosition().z() / m << G4endl;
>> 598 }
>> 599 #endif
>> 600 }
>> 601
>> 602 accuracy = theVelocityChange / c_light - 1.0;
>> 603 if(accuracy > accuracyForWarning)
>> 604 {
>> 605 itsOKforVelocity = false;
>> 606 nError += 1;
>> 607 exitWithError = exitWithError || (accuracy > accuracyForException);
>> 608 #ifdef G4VERBOSE
>> 609 if(nError < maxError)
>> 610 {
>> 611 G4cout << " G4ParticleChange::CheckIt : ";
>> 612 G4cout << "the velocity is greater than c_light !!" << G4endl;
>> 613 G4cout << " Velocity: " << theVelocityChange / c_light << G4endl;
>> 614 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 615 << " E=" << aTrack.GetKineticEnergy() / MeV
>> 616 << " pos=" << aTrack.GetPosition().x() / m << ", "
>> 617 << aTrack.GetPosition().y() / m << ", "
>> 618 << aTrack.GetPosition().z() / m << G4endl;
>> 619 }
>> 620 #endif
>> 621 }
>> 622
>> 623 G4bool itsOK = itsOKforMomentum && itsOKforEnergy && itsOKforVelocity &&
>> 624 itsOKforProperTime && itsOKforGlobalTime;
>> 625 // dump out information of this particle change
>> 626 #ifdef G4VERBOSE
>> 627 if(!itsOK)
>> 628 {
>> 629 DumpInfo();
>> 630 }
>> 631 #endif
>> 632
>> 633 // Exit with error
>> 634 if(exitWithError)
>> 635 {
>> 636 G4Exception("G4ParticleChange::CheckIt()", "TRACK003", EventMustBeAborted,
>> 637 "momentum, energy, and/or time was illegal");
>> 638 }
>> 639 // correction
>> 640 if(!itsOKforMomentum)
>> 641 {
>> 642 G4double vmag = theMomentumDirectionChange.mag();
>> 643 theMomentumDirectionChange = (1. / vmag) * theMomentumDirectionChange;
>> 644 }
>> 645 if(!itsOKforGlobalTime)
>> 646 {
>> 647 theTimeChange = aTrack.GetLocalTime();
>> 648 }
>> 649 if(!itsOKforProperTime)
>> 650 {
>> 651 theProperTimeChange = aTrack.GetProperTime();
>> 652 }
>> 653 if(!itsOKforEnergy)
>> 654 {
>> 655 theEnergyChange = 0.0;
>> 656 }
>> 657 if(!itsOKforVelocity)
>> 658 {
>> 659 theVelocityChange = c_light;
>> 660 }
>> 661
>> 662 itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
>> 663 return itsOK;
381 } 664 }
382 665