| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4ParticleChange class implementation <<
27 // 23 //
28 // Author: Hisaya Kurashige, 23 March 1998 << 24 // $Id: G4ParticleChange.cc,v 1.27 2004/12/02 06:38:02 kurasige Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-08-00-patch-01 $
>> 26 //
>> 27 //
>> 28 // --------------------------------------------------------------
>> 29 // GEANT 4 class implementation file
>> 30 //
>> 31 //
>> 32 //
>> 33 // ------------------------------------------------------------
>> 34 // Implemented for the new scheme 23 Mar. 1998 H.Kurahige
>> 35 // Change default debug flag to false 10 May. 1998 H.Kurahige
>> 36 // Add Track weight 12 Nov. 1998 H.Kurashige
>> 37 // Activate CheckIt method for VERBOSE mode 14 Dec. 1998 H.Kurashige
>> 38 // Modified CheckIt method for time 9 Feb. 1999 H.Kurashige
>> 39 // --------------------------------------------------------------
30 40
31 #include "G4ParticleChange.hh" 41 #include "G4ParticleChange.hh"
32 #include "G4PhysicalConstants.hh" <<
33 #include "G4SystemOfUnits.hh" <<
34 #include "G4Track.hh" 42 #include "G4Track.hh"
35 #include "G4Step.hh" 43 #include "G4Step.hh"
>> 44 #include "G4TrackFastVector.hh"
36 #include "G4DynamicParticle.hh" 45 #include "G4DynamicParticle.hh"
37 #include "G4ExceptionSeverity.hh" 46 #include "G4ExceptionSeverity.hh"
38 47
39 // ------------------------------------------- << 48
40 G4ParticleChange::G4ParticleChange() << 49 G4ParticleChange::G4ParticleChange():G4VParticleChange()
41 {} <<
42 <<
43 // ------------------------------------------- <<
44 void G4ParticleChange::AddSecondary(G4DynamicP <<
45 G4bool IsG <<
46 { 50 {
47 // create track << 51 }
48 G4Track* aTrack = new G4Track(aParticle, Get << 52
>> 53 G4ParticleChange::~G4ParticleChange()
>> 54 {
>> 55 #ifdef G4VERBOSE
>> 56 if (verboseLevel>2) {
>> 57 G4cout << "G4ParticleChange::~G4ParticleChange() " << G4endl;
>> 58 }
>> 59 #endif
>> 60 }
>> 61
>> 62 // copy constructor
>> 63 G4ParticleChange::G4ParticleChange(const G4ParticleChange &right): G4VParticleChange(right)
>> 64 {
>> 65 if (verboseLevel>1) {
>> 66 G4cout << "G4ParticleChange:: copy constructor is called " << G4endl;
>> 67 }
>> 68 theMomentumDirectionChange = right.theMomentumDirectionChange;
>> 69 thePolarizationChange = right.thePolarizationChange;
>> 70 thePositionChange = right.thePositionChange;
>> 71 theTimeChange = right.theTimeChange;
>> 72 theEnergyChange = right.theEnergyChange;
>> 73 theMassChange = right.theMassChange;
>> 74 theChargeChange = right.theChargeChange;
>> 75 theWeightChange = right.theWeightChange;
>> 76 theProperTimeChange = right.theProperTimeChange;
>> 77 }
>> 78
>> 79 // assignemnt operator
>> 80 G4ParticleChange & G4ParticleChange::operator=(const G4ParticleChange &right)
>> 81 {
>> 82 if (verboseLevel>1) {
>> 83 G4cout << "G4ParticleChange:: assignment operator is called " << G4endl;
>> 84 }
>> 85 if (this != &right)
>> 86 {
>> 87 theListOfSecondaries = right.theListOfSecondaries;
>> 88 theSizeOftheListOfSecondaries = right.theSizeOftheListOfSecondaries;
>> 89 theNumberOfSecondaries = right.theNumberOfSecondaries;
>> 90 theStatusChange = right.theStatusChange;
>> 91
>> 92 theMomentumDirectionChange = right.theMomentumDirectionChange;
>> 93 thePolarizationChange = right.thePolarizationChange;
>> 94 thePositionChange = right.thePositionChange;
>> 95 theTimeChange = right.theTimeChange;
>> 96 theEnergyChange = right.theEnergyChange;
>> 97 theMassChange = right.theMassChange;
>> 98 theChargeChange = right.theChargeChange;
>> 99 theWeightChange = right.theWeightChange;
>> 100
>> 101 theTrueStepLength = right.theTrueStepLength;
>> 102 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
>> 103 theSteppingControlFlag = right.theSteppingControlFlag;
>> 104 }
>> 105 return *this;
>> 106 }
>> 107
>> 108 G4bool G4ParticleChange::operator==(const G4ParticleChange &right) const
>> 109 {
>> 110 return ((G4VParticleChange *)this == (G4VParticleChange *) &right);
>> 111 }
>> 112
>> 113 G4bool G4ParticleChange::operator!=(const G4ParticleChange &right) const
>> 114 {
>> 115 return ((G4VParticleChange *)this != (G4VParticleChange *) &right);
>> 116 }
>> 117
>> 118
>> 119 //----------------------------------------------------------------
>> 120 // methods for handling secondaries
>> 121 //
>> 122
>> 123 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
>> 124 G4bool IsGoodForTracking )
>> 125 {
>> 126 // create track
>> 127 G4Track* aTrack = new G4Track(aParticle, theTimeChange, thePositionChange);
49 128
50 // set IsGoodGorTrackingFlag 129 // set IsGoodGorTrackingFlag
51 if(IsGoodForTracking) << 130 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
52 aTrack->SetGoodForTrackingFlag(); <<
53 131
54 // touchable handle is copied to keep the po << 132 // Touchable handle is copied to keep the pointer
55 aTrack->SetTouchableHandle(theCurrentTrack-> 133 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
56 << 134
57 // add a secondary << 135 // add a secondary
58 G4VParticleChange::AddSecondary(aTrack); 136 G4VParticleChange::AddSecondary(aTrack);
59 } 137 }
60 138
61 // ------------------------------------------- << 139 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
62 void G4ParticleChange::AddSecondary(G4DynamicP << 140 G4ThreeVector newPosition,
63 G4ThreeVec << 141 G4bool IsGoodForTracking )
64 G4bool IsG <<
65 { 142 {
66 // create track << 143 // create track
67 G4Track* aTrack = new G4Track(aParticle, Get << 144 G4Track* aTrack = new G4Track(aParticle, theTimeChange, newPosition);
68 145
69 // set IsGoodGorTrackingFlag 146 // set IsGoodGorTrackingFlag
70 if(IsGoodForTracking) << 147 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
71 aTrack->SetGoodForTrackingFlag(); <<
72 148
73 // touchable is a temporary object, so you c << 149 // Touchable is a temporary object, so you cannot keep the pointer
74 aTrack->SetTouchableHandle(nullptr); << 150 aTrack->SetTouchableHandle((G4VTouchable*)0);
75 151
76 // add a secondary << 152 // add a secondary
77 G4VParticleChange::AddSecondary(aTrack); 153 G4VParticleChange::AddSecondary(aTrack);
78 } 154 }
79 155
80 // ------------------------------------------- << 156 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
81 void G4ParticleChange::AddSecondary(G4DynamicP << 157 G4double newTime,
82 G4double n << 158 G4bool IsGoodForTracking )
83 { 159 {
84 // create track << 160 // create track
85 G4Track* aTrack = new G4Track(aParticle, new << 161 G4Track* aTrack = new G4Track(aParticle, newTime, thePositionChange);
86 162
87 // set IsGoodGorTrackingFlag 163 // set IsGoodGorTrackingFlag
88 if(IsGoodForTracking) << 164 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
89 aTrack->SetGoodForTrackingFlag(); << 165
90 << 166 // Touchable handle is copied to keep the pointer
91 // touchable handle is copied to keep the po <<
92 aTrack->SetTouchableHandle(theCurrentTrack-> 167 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
93 168
94 // add a secondary << 169 // add a secondary
95 G4VParticleChange::AddSecondary(aTrack); 170 G4VParticleChange::AddSecondary(aTrack);
96 } 171 }
97 172
98 // ------------------------------------------- <<
99 <<
100 void G4ParticleChange::AddSecondary(G4Track* a 173 void G4ParticleChange::AddSecondary(G4Track* aTrack)
101 { 174 {
102 // add a secondary << 175 // add a secondary
103 G4VParticleChange::AddSecondary(aTrack); 176 G4VParticleChange::AddSecondary(aTrack);
104 } 177 }
105 178
106 // ------------------------------------------- << 179 //----------------------------------------------------------------
>> 180 // functions for Initialization
>> 181 //
>> 182
107 void G4ParticleChange::Initialize(const G4Trac 183 void G4ParticleChange::Initialize(const G4Track& track)
108 { 184 {
109 // use base class's method at first 185 // use base class's method at first
110 G4VParticleChange::Initialize(track); 186 G4VParticleChange::Initialize(track);
>> 187 theCurrentTrack= &track;
111 188
112 // set Energy/Momentum etc. equal to those o 189 // set Energy/Momentum etc. equal to those of the parent particle
113 const G4DynamicParticle* pParticle = track.G << 190 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
114 theEnergyChange = pPartic << 191 theEnergyChange = pParticle->GetKineticEnergy();
115 theVelocityChange = track.G << 192 theMomentumDirectionChange = pParticle->GetMomentumDirection();
116 isVelocityChanged = false; << 193 thePolarizationChange = pParticle->GetPolarization();
117 theMomentumDirectionChange = pPartic << 194 theProperTimeChange = pParticle->GetProperTime();
118 thePolarizationChange = pPartic << 195
119 theProperTimeChange = pPartic << 196 // Set mass/charge of DynamicParticle
120 << 197 theMassChange = pParticle->GetMass();
121 // set mass/charge/MagneticMoment of Dynamic << 198 theChargeChange = pParticle->GetCharge();
122 theMassChange = pParticle->GetMass << 199
123 theChargeChange = pParticle->GetChar << 200 // set Position/Time etc. equal to those of the parent track
124 theMagneticMomentChange = pParticle->GetMagn << 201 thePositionChange = track.GetPosition();
125 << 202 theTimeChange = track.GetGlobalTime();
126 // set Position equal to those of the parent <<
127 thePositionChange = track.GetPosition(); <<
128 <<
129 // set TimeChange equal to local time of the <<
130 theTimeChange = theLocalTime0 = track.GetLoc <<
131 203
132 // set initial global time of the parent tra << 204 theWeightChange = track.GetWeight();
133 theGlobalTime0 = track.GetGlobalTime(); <<
134 } 205 }
135 206
136 // ------------------------------------------- << 207 //----------------------------------------------------------------
>> 208 // methods for updating G4Step
>> 209 //
>> 210
137 G4Step* G4ParticleChange::UpdateStepForAlongSt 211 G4Step* G4ParticleChange::UpdateStepForAlongStep(G4Step* pStep)
138 { 212 {
139 // A physics process always calculates the f 213 // A physics process always calculates the final state of the
140 // particle relative to the initial state at 214 // particle relative to the initial state at the beginning
141 // of the Step, i.e., based on information o 215 // of the Step, i.e., based on information of G4Track (or
142 // equivalently the PreStepPoint). << 216 // equivalently the PreStepPoint).
143 // So, the differences (delta) between these 217 // So, the differences (delta) between these two states have to be
144 // calculated and be accumulated in PostStep << 218 // calculated and be accumulated in PostStepPoint.
145 << 219
146 // Take note that the return type of GetMome << 220 // Take note that the return type of GetMomentumDirectionChange is a
147 // pointer to G4ParticleMometum. Also it is << 221 // pointer to G4ParticleMometum. Also it is a normalized
148 // momentum vector << 222 // momentum vector.
149 << 223
150 const G4StepPoint* pPreStepPoint = pStep->Ge << 224 G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
151 G4StepPoint* pPostStepPoint = pStep->GetPost << 225 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
>> 226 G4double mass = theMassChange;
152 227
153 // set Mass/Charge/MagneticMoment << 228 // Set Mass/Charge
154 pPostStepPoint->SetMass(theMassChange); 229 pPostStepPoint->SetMass(theMassChange);
155 pPostStepPoint->SetCharge(theChargeChange); << 230 pPostStepPoint->SetCharge(theChargeChange);
156 pPostStepPoint->SetMagneticMoment(theMagneti << 231
157 <<
158 // calculate new kinetic energy 232 // calculate new kinetic energy
159 G4double preEnergy = pPreStepPoint->GetKinet << 233 G4double energy = pPostStepPoint->GetKineticEnergy()
160 G4double energy = << 234 + (theEnergyChange - pPreStepPoint->GetKineticEnergy());
161 pPostStepPoint->GetKineticEnergy() + (theE <<
162 235
163 // update kinetic energy and momentum direct 236 // update kinetic energy and momentum direction
164 if(energy > 0.0) << 237 if (energy > 0.0) {
165 { <<
166 // calculate new momentum 238 // calculate new momentum
167 G4ThreeVector pMomentum = pPostStepPoint-> << 239 G4ThreeVector pMomentum = pPostStepPoint->GetMomentum()
168 + (CalcMomentum(theEnergyChange, theMome << 240 + ( CalcMomentum(theEnergyChange, theMomentumDirectionChange, mass)
169 theMassChange) - pPreStepPoint->GetM << 241 - pPreStepPoint->GetMomentum());
170 G4double tMomentum2 = pMomentum.mag2(); << 242 G4double tMomentum = pMomentum.mag();
171 G4ThreeVector direction(1.0, 0.0, 0.0); << 243 G4ThreeVector direction(1.0,0.0,0.0);
172 if(tMomentum2 > 0.) << 244 if( tMomentum > 0. ){
173 { << 245 G4double inv_Momentum= 1.0 / tMomentum;
174 direction = pMomentum / std::sqrt(tMomen << 246 direction= pMomentum * inv_Momentum;
175 } 247 }
176 pPostStepPoint->SetMomentumDirection(direc 248 pPostStepPoint->SetMomentumDirection(direction);
177 pPostStepPoint->SetKineticEnergy(energy); << 249 pPostStepPoint->SetKineticEnergy( energy );
178 << 250 } else {
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 } <<
205 else <<
206 { <<
207 // stop case 251 // stop case
>> 252 //pPostStepPoint->SetMomentumDirection(G4ThreeVector(1., 0., 0.));
208 pPostStepPoint->SetKineticEnergy(0.0); 253 pPostStepPoint->SetKineticEnergy(0.0);
209 pPostStepPoint->SetVelocity(0.0); <<
210 } 254 }
211 255
212 // update polarization 256 // update polarization
213 pPostStepPoint->AddPolarization(thePolarizat << 257 pPostStepPoint->AddPolarization( thePolarizationChange
214 pPreStepPoin << 258 - pPreStepPoint->GetPolarization());
215 << 259
216 // update position and time 260 // update position and time
217 pPostStepPoint->AddPosition(thePositionChang << 261 pPostStepPoint->AddPosition( thePositionChange
218 pPostStepPoint->AddGlobalTime(theTimeChange << 262 - pPreStepPoint->GetPosition() );
219 pPostStepPoint->AddLocalTime(theTimeChange - << 263 pPostStepPoint->AddGlobalTime( theTimeChange
220 pPostStepPoint->AddProperTime(theProperTimeC << 264 - pPreStepPoint->GetGlobalTime());
221 pPreStepPoint- << 265 pPostStepPoint->AddLocalTime( theTimeChange
222 << 266 - pPreStepPoint->GetGlobalTime());
223 if(isParentWeightProposed) << 267 pPostStepPoint->AddProperTime( theProperTimeChange
224 { << 268 - pPreStepPoint->GetProperTime());
225 pPostStepPoint->SetWeight(theParentWeight) << 269
226 } << 270 // update weight
>> 271 G4double newWeight= theWeightChange/(pPreStepPoint->GetWeight())*(pPostStepPoint->GetWeight());
>> 272 pPostStepPoint->SetWeight( newWeight );
227 273
228 #ifdef G4VERBOSE 274 #ifdef G4VERBOSE
229 if(debugFlag) { CheckIt( *theCurrentTrack ); << 275 G4Track* aTrack = pStep->GetTrack();
>> 276 if (debugFlag) CheckIt(*aTrack);
230 #endif 277 #endif
231 278
232 // update the G4Step specific attributes << 279 // Update the G4Step specific attributes
233 return UpdateStepInfo(pStep); 280 return UpdateStepInfo(pStep);
234 } 281 }
235 282
236 // ------------------------------------------- <<
237 G4Step* G4ParticleChange::UpdateStepForPostSte 283 G4Step* G4ParticleChange::UpdateStepForPostStep(G4Step* pStep)
238 { << 284 {
239 // A physics process always calculates the f 285 // A physics process always calculates the final state of the particle
240 286
241 // Take note that the return type of GetMome << 287 // Take note that the return type of GetMomentumChange is a
242 // pointer to G4ParticleMometum. Also it is << 288 // pointer to G4ParticleMometum. Also it is a normalized
243 // momentum vector << 289 // momentum vector.
244 290
245 G4StepPoint* pPostStepPoint = pStep->GetPost << 291 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
246 G4Track* pTrack = pStep->GetTrac << 292 G4Track* aTrack = pStep->GetTrack();
247 293
248 // set Mass/Charge << 294 // Set Mass/Charge
249 pPostStepPoint->SetMass(theMassChange); 295 pPostStepPoint->SetMass(theMassChange);
250 pPostStepPoint->SetCharge(theChargeChange); << 296 pPostStepPoint->SetCharge(theChargeChange);
251 pPostStepPoint->SetMagneticMoment(theMagneti << 297
252 <<
253 // update kinetic energy and momentum direct 298 // update kinetic energy and momentum direction
254 pPostStepPoint->SetMomentumDirection(theMome 299 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
>> 300 pPostStepPoint->SetKineticEnergy( theEnergyChange );
255 301
256 // calculate velocity << 302 // update polarization
257 if(theEnergyChange > 0.0) << 303 pPostStepPoint->SetPolarization( thePolarizationChange );
258 { << 304
259 pPostStepPoint->SetKineticEnergy(theEnergy <<
260 pTrack->SetKineticEnergy(theEnergyChange); <<
261 if(!isVelocityChanged) <<
262 { <<
263 theVelocityChange = pTrack->CalculateVel <<
264 } <<
265 pPostStepPoint->SetVelocity(theVelocityCha <<
266 } <<
267 else <<
268 { <<
269 pPostStepPoint->SetKineticEnergy(0.0); <<
270 pPostStepPoint->SetVelocity(0.0); <<
271 } <<
272 <<
273 // update polarization <<
274 pPostStepPoint->SetPolarization(thePolarizat <<
275 <<
276 // update position and time 305 // update position and time
277 pPostStepPoint->SetPosition(thePositionChang << 306 pPostStepPoint->SetPosition( thePositionChange );
278 pPostStepPoint->AddGlobalTime(theTimeChange << 307 pPostStepPoint->SetGlobalTime( theTimeChange );
279 pPostStepPoint->SetLocalTime(theTimeChange); << 308 pPostStepPoint->AddLocalTime( theTimeChange
280 pPostStepPoint->SetProperTime(theProperTimeC << 309 - aTrack->GetGlobalTime());
281 << 310 pPostStepPoint->SetProperTime( theProperTimeChange );
282 if(isParentWeightProposed) << 311
283 { << 312 // update weight
284 pPostStepPoint->SetWeight(theParentWeight) << 313 pPostStepPoint->SetWeight( theWeightChange );
285 } <<
286 314
287 #ifdef G4VERBOSE 315 #ifdef G4VERBOSE
288 if(debugFlag) { CheckIt( *theCurrentTrack ); << 316 if (debugFlag) CheckIt(*aTrack);
289 #endif 317 #endif
290 318
291 // update the G4Step specific attributes << 319 // Update the G4Step specific attributes
292 return UpdateStepInfo(pStep); 320 return UpdateStepInfo(pStep);
293 } 321 }
294 322
295 // ------------------------------------------- << 323
296 G4Step* G4ParticleChange::UpdateStepForAtRest( 324 G4Step* G4ParticleChange::UpdateStepForAtRest(G4Step* pStep)
297 { << 325 {
298 // A physics process always calculates the f 326 // A physics process always calculates the final state of the particle
299 327
300 G4StepPoint* pPostStepPoint = pStep->GetPost << 328 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
>> 329 G4Track* aTrack = pStep->GetTrack();
301 330
302 // set Mass/Charge << 331 // Set Mass/Charge
303 pPostStepPoint->SetMass(theMassChange); 332 pPostStepPoint->SetMass(theMassChange);
304 pPostStepPoint->SetCharge(theChargeChange); << 333 pPostStepPoint->SetCharge(theChargeChange);
305 pPostStepPoint->SetMagneticMoment(theMagneti << 334
306 <<
307 // update kinetic energy and momentum direct 335 // update kinetic energy and momentum direction
308 pPostStepPoint->SetMomentumDirection(theMome 336 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
309 pPostStepPoint->SetKineticEnergy(theEnergyCh << 337 pPostStepPoint->SetKineticEnergy( theEnergyChange );
310 if(!isVelocityChanged) <<
311 theVelocityChange = theCurrentTrack->Calcu <<
312 pPostStepPoint->SetVelocity(theVelocityChang <<
313 338
314 // update polarization 339 // update polarization
315 pPostStepPoint->SetPolarization(thePolarizat << 340 pPostStepPoint->SetPolarization( thePolarizationChange );
316 << 341
317 // update position and time 342 // update position and time
318 pPostStepPoint->SetPosition(thePositionChang << 343 pPostStepPoint->SetPosition( thePositionChange );
319 pPostStepPoint->AddGlobalTime(theTimeChange << 344 pPostStepPoint->SetGlobalTime( theTimeChange );
320 pPostStepPoint->SetLocalTime(theTimeChange); << 345 pPostStepPoint->AddLocalTime( theTimeChange
321 pPostStepPoint->SetProperTime(theProperTimeC << 346 - aTrack->GetGlobalTime());
322 << 347 pPostStepPoint->SetProperTime( theProperTimeChange );
323 if(isParentWeightProposed) << 348
324 { << 349 // update weight
325 pPostStepPoint->SetWeight(theParentWeight) << 350 pPostStepPoint->SetWeight( theWeightChange );
326 } <<
327 351
328 #ifdef G4VERBOSE 352 #ifdef G4VERBOSE
329 if(debugFlag) { CheckIt( *theCurrentTrack ); << 353 if (debugFlag) CheckIt(*aTrack);
330 #endif 354 #endif
331 355
332 // update the G4Step specific attributes << 356 // Update the G4Step specific attributes
333 return UpdateStepInfo(pStep); 357 return UpdateStepInfo(pStep);
334 } 358 }
335 359
336 // ------------------------------------------- << 360 //----------------------------------------------------------------
>> 361 // methods for printing messages
>> 362 //
>> 363
337 void G4ParticleChange::DumpInfo() const 364 void G4ParticleChange::DumpInfo() const
338 { 365 {
339 // use base-class DumpInfo << 366 // use base-class DumpInfo
340 G4VParticleChange::DumpInfo(); 367 G4VParticleChange::DumpInfo();
341 368
342 G4long oldprc = G4cout.precision(8); << 369 G4cout.precision(3);
343 370
344 G4cout << " Mass (GeV) : " < << 371 G4cout << " Mass (GeV) : "
345 << G4endl; << 372 << std::setw(20) << theMassChange/GeV
346 G4cout << " Charge (eplus) : " < << 373 << G4endl;
347 << theChargeChange / eplus << G4endl; << 374 G4cout << " Charge (eplus) : "
348 G4cout << " MagneticMoment : " < << 375 << std::setw(20) << theChargeChange/eplus
349 << theMagneticMomentChange << G4endl; << 376 << G4endl;
350 G4cout << " = : " < << 377 G4cout << " Position - x (mm) : "
351 << theMagneticMomentChange << 378 << std::setw(20) << thePositionChange.x()/mm
352 * 2. * theMassChange / c_squared / << 379 << G4endl;
353 << "*[e hbar]/[2 m]" << G4endl; << 380 G4cout << " Position - y (mm) : "
354 G4cout << " Position - x (mm) : " < << 381 << std::setw(20) << thePositionChange.y()/mm
355 << thePositionChange.x() / mm << G4en << 382 << G4endl;
356 G4cout << " Position - y (mm) : " < << 383 G4cout << " Position - z (mm) : "
357 << thePositionChange.y() / mm << G4en << 384 << std::setw(20) << thePositionChange.z()/mm
358 G4cout << " Position - z (mm) : " < << 385 << G4endl;
359 << thePositionChange.z() / mm << G4en << 386 G4cout << " Time (ns) : "
360 G4cout << " Time (ns) : " < << 387 << std::setw(20) << theTimeChange/ns
361 << theTimeChange / ns << G4endl; << 388 << G4endl;
362 G4cout << " Proper Time (ns) : " < << 389 G4cout << " Proper Time (ns) : "
363 << theProperTimeChange / ns << G4endl << 390 << std::setw(20) << theProperTimeChange/ns
364 G4cout << " Momentum Direct - x : " < << 391 << G4endl;
365 << theMomentumDirectionChange.x() << << 392 G4cout << " Momentum Direct - x : "
366 G4cout << " Momentum Direct - y : " < << 393 << std::setw(20) << theMomentumDirectionChange.x()
367 << theMomentumDirectionChange.y() << << 394 << G4endl;
368 G4cout << " Momentum Direct - z : " < << 395 G4cout << " Momentum Direct - y : "
369 << theMomentumDirectionChange.z() << << 396 << std::setw(20) << theMomentumDirectionChange.y()
370 G4cout << " Kinetic Energy (MeV): " < << 397 << G4endl;
371 << theEnergyChange / MeV << G4endl; << 398 G4cout << " Momentum Direct - z : "
372 G4cout << " Velocity (/c) : " < << 399 << std::setw(20) << theMomentumDirectionChange.z()
373 << theVelocityChange / c_light << G4e << 400 << G4endl;
374 G4cout << " Polarization - x : " < << 401 G4cout << " Kinetic Energy (MeV): "
375 << thePolarizationChange.x() << G4end << 402 << std::setw(20) << theEnergyChange/MeV
376 G4cout << " Polarization - y : " < << 403 << G4endl;
377 << thePolarizationChange.y() << G4end << 404 G4cout << " Polarization - x : "
378 G4cout << " Polarization - z : " < << 405 << std::setw(20) << thePolarizationChange.x()
379 << thePolarizationChange.z() << G4end << 406 << G4endl;
380 G4cout.precision(oldprc); << 407 G4cout << " Polarization - y : "
>> 408 << std::setw(20) << thePolarizationChange.y()
>> 409 << G4endl;
>> 410 G4cout << " Polarization - z : "
>> 411 << std::setw(20) << thePolarizationChange.z()
>> 412 << G4endl;
>> 413 G4cout << " Track Weight : "
>> 414 << std::setw(20) << theWeightChange
>> 415 << G4endl;
381 } 416 }
>> 417
>> 418 G4bool G4ParticleChange::CheckIt(const G4Track& aTrack)
>> 419 {
>> 420 G4bool exitWithError = false;
>> 421 G4double accuracy;
>> 422
>> 423 // No check in case of "fStopAndKill"
>> 424 if (GetTrackStatus() == fStopAndKill ) {
>> 425 return G4VParticleChange::CheckIt(aTrack);
>> 426 }
>> 427
>> 428 // MomentumDirection should be unit vector
>> 429 G4bool itsOKforMomentum = true;
>> 430 if ( theEnergyChange >0.) {
>> 431 accuracy = std::abs(theMomentumDirectionChange.mag2()-1.0);
>> 432 if (accuracy > accuracyForWarning) {
>> 433 #ifdef G4VERBOSE
>> 434 G4cout << " G4ParticleChange::CheckIt : ";
>> 435 G4cout << "the Momentum Change is not unit vector !!" << G4endl;
>> 436 G4cout << " Difference: " << accuracy << G4endl;
>> 437 #endif
>> 438 itsOKforMomentum = false;
>> 439 if (accuracy > accuracyForException) exitWithError = true;
>> 440 }
>> 441 }
>> 442
>> 443 // Both global and proper time should not go back
>> 444 G4bool itsOKforGlobalTime = true;
>> 445 accuracy = (aTrack.GetGlobalTime()- theTimeChange)/ns;
>> 446 if (accuracy > accuracyForWarning) {
>> 447 #ifdef G4VERBOSE
>> 448 G4cout << " G4ParticleChange::CheckIt : ";
>> 449 G4cout << "the global time goes back !!" << G4endl;
>> 450 G4cout << " Difference: " << accuracy << "[ns] " <<G4endl;
>> 451 #endif
>> 452 itsOKforGlobalTime = false;
>> 453 if (accuracy > accuracyForException) exitWithError = true;
>> 454 }
>> 455
>> 456 G4bool itsOKforProperTime = true;
>> 457 accuracy = (aTrack.GetProperTime() - theProperTimeChange )/ns;
>> 458 if (accuracy > accuracyForWarning) {
>> 459 #ifdef G4VERBOSE
>> 460 G4cout << " G4ParticleChange::CheckIt : ";
>> 461 G4cout << "the proper time goes back !!" << G4endl;
>> 462 G4cout << " Difference: " << accuracy << "[ns] " <<G4endl;
>> 463 #endif
>> 464 itsOKforProperTime = false;
>> 465 if (accuracy > accuracyForException) exitWithError = true;
>> 466 }
>> 467
>> 468 // Kinetic Energy should not be negative
>> 469 G4bool itsOKforEnergy = true;
>> 470 accuracy = -1.0*theEnergyChange/MeV;
>> 471 if (accuracy > accuracyForWarning) {
>> 472 #ifdef G4VERBOSE
>> 473 G4cout << " G4ParticleChange::CheckIt : ";
>> 474 G4cout << "the kinetic energy is negative !!" << G4endl;
>> 475 G4cout << " Difference: " << accuracy << "[MeV] " <<G4endl;
>> 476 #endif
>> 477 itsOKforEnergy = false;
>> 478 if (accuracy > accuracyForException) exitWithError = true;
>> 479 }
>> 480
>> 481 G4bool itsOK = itsOKforMomentum && itsOKforEnergy && itsOKforProperTime && itsOKforGlobalTime;
>> 482 // dump out information of this particle change
>> 483 #ifdef G4VERBOSE
>> 484 if (!itsOK) {
>> 485 G4cout << " G4ParticleChange::CheckIt " <<G4endl;
>> 486 DumpInfo();
>> 487 }
>> 488 #endif
>> 489
>> 490 // Exit with error
>> 491 if (exitWithError) {
>> 492 G4Exception("G4ParticleChange::CheckIt",
>> 493 "200",
>> 494 EventMustBeAborted,
>> 495 "momentum, energy, and/or time was illegal");
>> 496 }
>> 497 //correction
>> 498 if (!itsOKforMomentum) {
>> 499 G4double vmag = theMomentumDirectionChange.mag();
>> 500 theMomentumDirectionChange = (1./vmag)*theMomentumDirectionChange;
>> 501 }
>> 502 if (!itsOKforGlobalTime) {
>> 503 theTimeChange = aTrack.GetGlobalTime();
>> 504 }
>> 505 if (!itsOKforProperTime) {
>> 506 theProperTimeChange = aTrack.GetProperTime();
>> 507 }
>> 508 if (!itsOKforEnergy) {
>> 509 theEnergyChange = 0.0;
>> 510 }
>> 511
>> 512 itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
>> 513 return itsOK;
>> 514 }
>> 515
>> 516
>> 517
382 518