| 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.21 2001/12/24 05:14:40 kurasige Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-04-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" <<
38 46
39 // ------------------------------------------- << 47
40 G4ParticleChange::G4ParticleChange() << 48 G4ParticleChange::G4ParticleChange():G4VParticleChange()
41 {} << 49 {
42 << 50 }
43 // ------------------------------------------- << 51
44 void G4ParticleChange::AddSecondary(G4DynamicP << 52 G4ParticleChange::~G4ParticleChange()
45 G4bool IsG << 53 {
>> 54 #ifdef G4VERBOSE
>> 55 if (verboseLevel>2) {
>> 56 G4cout << "G4ParticleChange::~G4ParticleChange() " << G4endl;
>> 57 }
>> 58 #endif
>> 59 }
>> 60
>> 61 // copy constructor
>> 62 G4ParticleChange::G4ParticleChange(const G4ParticleChange &right): G4VParticleChange(right)
>> 63 {
>> 64 if (verboseLevel>1) {
>> 65 G4cout << "G4ParticleChange:: copy constructor is called " << G4endl;
>> 66 }
>> 67 theMomentumDirectionChange = right.theMomentumDirectionChange;
>> 68 thePolarizationChange = right.thePolarizationChange;
>> 69 thePositionChange = right.thePositionChange;
>> 70 theTimeChange = right.theTimeChange;
>> 71 theEnergyChange = right.theEnergyChange;
>> 72 theMassChange = right.theMassChange;
>> 73 theChargeChange = right.theChargeChange;
>> 74 theWeightChange = right.theWeightChange;
>> 75 theProperTimeChange = right.theProperTimeChange;
>> 76 }
>> 77
>> 78 // assignemnt operator
>> 79 G4ParticleChange & G4ParticleChange::operator=(const G4ParticleChange &right)
>> 80 {
>> 81 if (verboseLevel>1) {
>> 82 G4cout << "G4ParticleChange:: assignment operator is called " << G4endl;
>> 83 }
>> 84 if (this != &right)
>> 85 {
>> 86 theListOfSecondaries = right.theListOfSecondaries;
>> 87 theSizeOftheListOfSecondaries = right.theSizeOftheListOfSecondaries;
>> 88 theNumberOfSecondaries = right.theNumberOfSecondaries;
>> 89 theStatusChange = right.theStatusChange;
>> 90
>> 91 theMomentumDirectionChange = right.theMomentumDirectionChange;
>> 92 thePolarizationChange = right.thePolarizationChange;
>> 93 thePositionChange = right.thePositionChange;
>> 94 theTimeChange = right.theTimeChange;
>> 95 theEnergyChange = right.theEnergyChange;
>> 96 theMassChange = right.theMassChange;
>> 97 theChargeChange = right.theChargeChange;
>> 98 theWeightChange = right.theWeightChange;
>> 99
>> 100 theTrueStepLength = right.theTrueStepLength;
>> 101 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
>> 102 theSteppingControlFlag = right.theSteppingControlFlag;
>> 103 }
>> 104 return *this;
>> 105 }
>> 106
>> 107 G4bool G4ParticleChange::operator==(const G4ParticleChange &right) const
>> 108 {
>> 109 return ((G4VParticleChange *)this == (G4VParticleChange *) &right);
>> 110 }
>> 111
>> 112 G4bool G4ParticleChange::operator!=(const G4ParticleChange &right) const
>> 113 {
>> 114 return ((G4VParticleChange *)this != (G4VParticleChange *) &right);
>> 115 }
>> 116
>> 117
>> 118 //----------------------------------------------------------------
>> 119 // methods for handling secondaries
>> 120 //
>> 121
>> 122 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
>> 123 G4bool IsGoodForTracking )
46 { 124 {
47 // create track << 125 // create track
48 G4Track* aTrack = new G4Track(aParticle, Get << 126 G4Track* aTrack = new G4Track(aParticle, theTimeChange, thePositionChange);
49 127
50 // set IsGoodGorTrackingFlag 128 // set IsGoodGorTrackingFlag
51 if(IsGoodForTracking) << 129 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
52 aTrack->SetGoodForTrackingFlag(); <<
53 130
54 // touchable handle is copied to keep the po << 131 // Touchable handle is copied to keep the pointer
55 aTrack->SetTouchableHandle(theCurrentTrack-> 132 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
56 << 133
57 // add a secondary << 134 // add a secondary
58 G4VParticleChange::AddSecondary(aTrack); 135 G4VParticleChange::AddSecondary(aTrack);
59 } 136 }
60 137
61 // ------------------------------------------- << 138 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
62 void G4ParticleChange::AddSecondary(G4DynamicP << 139 G4ThreeVector newPosition,
63 G4ThreeVec << 140 G4bool IsGoodForTracking )
64 G4bool IsG <<
65 { 141 {
66 // create track << 142 // create track
67 G4Track* aTrack = new G4Track(aParticle, Get << 143 G4Track* aTrack = new G4Track(aParticle, theTimeChange, newPosition);
68 144
69 // set IsGoodGorTrackingFlag 145 // set IsGoodGorTrackingFlag
70 if(IsGoodForTracking) << 146 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
71 aTrack->SetGoodForTrackingFlag(); <<
72 147
73 // touchable is a temporary object, so you c << 148 // Touchable is a temporary object, so you cannot keep the pointer
74 aTrack->SetTouchableHandle(nullptr); << 149 aTrack->SetTouchableHandle((G4VTouchable*)0);
75 150
76 // add a secondary << 151 // add a secondary
77 G4VParticleChange::AddSecondary(aTrack); 152 G4VParticleChange::AddSecondary(aTrack);
78 } 153 }
79 154
80 // ------------------------------------------- << 155 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
81 void G4ParticleChange::AddSecondary(G4DynamicP << 156 G4double newTime,
82 G4double n << 157 G4bool IsGoodForTracking )
83 { 158 {
84 // create track << 159 // create track
85 G4Track* aTrack = new G4Track(aParticle, new << 160 G4Track* aTrack = new G4Track(aParticle, newTime, thePositionChange);
86 161
87 // set IsGoodGorTrackingFlag 162 // set IsGoodGorTrackingFlag
88 if(IsGoodForTracking) << 163 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
89 aTrack->SetGoodForTrackingFlag(); << 164
90 << 165 // Touchable handle is copied to keep the pointer
91 // touchable handle is copied to keep the po <<
92 aTrack->SetTouchableHandle(theCurrentTrack-> 166 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
93 167
94 // add a secondary << 168 // add a secondary
95 G4VParticleChange::AddSecondary(aTrack); 169 G4VParticleChange::AddSecondary(aTrack);
96 } 170 }
97 171
98 // ------------------------------------------- <<
99 <<
100 void G4ParticleChange::AddSecondary(G4Track* a 172 void G4ParticleChange::AddSecondary(G4Track* aTrack)
101 { 173 {
102 // add a secondary << 174 // add a secondary
103 G4VParticleChange::AddSecondary(aTrack); 175 G4VParticleChange::AddSecondary(aTrack);
104 } 176 }
105 177
106 // ------------------------------------------- << 178 //----------------------------------------------------------------
>> 179 // functions for Initialization
>> 180 //
>> 181
107 void G4ParticleChange::Initialize(const G4Trac 182 void G4ParticleChange::Initialize(const G4Track& track)
108 { 183 {
109 // use base class's method at first 184 // use base class's method at first
110 G4VParticleChange::Initialize(track); 185 G4VParticleChange::Initialize(track);
>> 186 theCurrentTrack= &track;
111 187
112 // set Energy/Momentum etc. equal to those o 188 // set Energy/Momentum etc. equal to those of the parent particle
113 const G4DynamicParticle* pParticle = track.G << 189 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
114 theEnergyChange = pPartic << 190 theEnergyChange = pParticle->GetKineticEnergy();
115 theVelocityChange = track.G << 191 theMomentumDirectionChange = pParticle->GetMomentumDirection();
116 isVelocityChanged = false; << 192 thePolarizationChange = pParticle->GetPolarization();
117 theMomentumDirectionChange = pPartic << 193 theProperTimeChange = pParticle->GetProperTime();
118 thePolarizationChange = pPartic << 194
119 theProperTimeChange = pPartic << 195 // Set mass/charge of DynamicParticle
120 << 196 theMassChange = pParticle->GetMass();
121 // set mass/charge/MagneticMoment of Dynamic << 197 theChargeChange = pParticle->GetCharge();
122 theMassChange = pParticle->GetMass << 198
123 theChargeChange = pParticle->GetChar << 199 // set Position/Time etc. equal to those of the parent track
124 theMagneticMomentChange = pParticle->GetMagn << 200 thePositionChange = track.GetPosition();
125 << 201 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 202
132 // set initial global time of the parent tra << 203 theWeightChange = track.GetWeight();
133 theGlobalTime0 = track.GetGlobalTime(); <<
134 } 204 }
135 205
136 // ------------------------------------------- << 206 //----------------------------------------------------------------
>> 207 // methods for updating G4Step
>> 208 //
>> 209
137 G4Step* G4ParticleChange::UpdateStepForAlongSt 210 G4Step* G4ParticleChange::UpdateStepForAlongStep(G4Step* pStep)
138 { 211 {
139 // A physics process always calculates the f 212 // A physics process always calculates the final state of the
140 // particle relative to the initial state at 213 // particle relative to the initial state at the beginning
141 // of the Step, i.e., based on information o 214 // of the Step, i.e., based on information of G4Track (or
142 // equivalently the PreStepPoint). << 215 // equivalently the PreStepPoint).
143 // So, the differences (delta) between these 216 // So, the differences (delta) between these two states have to be
144 // calculated and be accumulated in PostStep << 217 // calculated and be accumulated in PostStepPoint.
145 << 218
146 // Take note that the return type of GetMome << 219 // Take note that the return type of GetMomentumDirectionChange is a
147 // pointer to G4ParticleMometum. Also it is << 220 // pointer to G4ParticleMometum. Also it is a normalized
148 // momentum vector << 221 // momentum vector.
>> 222
>> 223 G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
>> 224 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
>> 225 G4Track* aTrack = pStep->GetTrack();
>> 226 G4double mass = theMassChange;
149 227
150 const G4StepPoint* pPreStepPoint = pStep->Ge << 228 // Set Mass/Charge
151 G4StepPoint* pPostStepPoint = pStep->GetPost <<
152 <<
153 // set Mass/Charge/MagneticMoment <<
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 if (debugFlag) CheckIt(*aTrack);
230 #endif 276 #endif
231 277
232 // update the G4Step specific attributes << 278 // Update the G4Step specific attributes
233 return UpdateStepInfo(pStep); 279 return UpdateStepInfo(pStep);
234 } 280 }
235 281
236 // ------------------------------------------- <<
237 G4Step* G4ParticleChange::UpdateStepForPostSte 282 G4Step* G4ParticleChange::UpdateStepForPostStep(G4Step* pStep)
238 { << 283 {
239 // A physics process always calculates the f 284 // A physics process always calculates the final state of the particle
240 285
241 // Take note that the return type of GetMome << 286 // Take note that the return type of GetMomentumChange is a
242 // pointer to G4ParticleMometum. Also it is << 287 // pointer to G4ParticleMometum. Also it is a normalized
243 // momentum vector << 288 // momentum vector.
244 289
245 G4StepPoint* pPostStepPoint = pStep->GetPost << 290 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
246 G4Track* pTrack = pStep->GetTrac << 291 G4Track* aTrack = pStep->GetTrack();
247 292
248 // set Mass/Charge << 293 // Set Mass/Charge
249 pPostStepPoint->SetMass(theMassChange); 294 pPostStepPoint->SetMass(theMassChange);
250 pPostStepPoint->SetCharge(theChargeChange); << 295 pPostStepPoint->SetCharge(theChargeChange);
251 pPostStepPoint->SetMagneticMoment(theMagneti << 296
252 <<
253 // update kinetic energy and momentum direct 297 // update kinetic energy and momentum direction
254 pPostStepPoint->SetMomentumDirection(theMome 298 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
>> 299 pPostStepPoint->SetKineticEnergy( theEnergyChange );
255 300
256 // calculate velocity << 301 // update polarization
257 if(theEnergyChange > 0.0) << 302 pPostStepPoint->SetPolarization( thePolarizationChange );
258 { << 303
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 304 // update position and time
277 pPostStepPoint->SetPosition(thePositionChang << 305 pPostStepPoint->SetPosition( thePositionChange );
278 pPostStepPoint->AddGlobalTime(theTimeChange << 306 pPostStepPoint->SetGlobalTime( theTimeChange );
279 pPostStepPoint->SetLocalTime(theTimeChange); << 307 pPostStepPoint->AddLocalTime( theTimeChange
280 pPostStepPoint->SetProperTime(theProperTimeC << 308 - aTrack->GetGlobalTime());
281 << 309 pPostStepPoint->SetProperTime( theProperTimeChange );
282 if(isParentWeightProposed) << 310
283 { << 311 // update weight
284 pPostStepPoint->SetWeight(theParentWeight) << 312 pPostStepPoint->SetWeight( theWeightChange );
285 } <<
286 313
287 #ifdef G4VERBOSE 314 #ifdef G4VERBOSE
288 if(debugFlag) { CheckIt( *theCurrentTrack ); << 315 if (debugFlag) CheckIt(*aTrack);
289 #endif 316 #endif
290 317
291 // update the G4Step specific attributes << 318 // Update the G4Step specific attributes
292 return UpdateStepInfo(pStep); 319 return UpdateStepInfo(pStep);
293 } 320 }
294 321
295 // ------------------------------------------- << 322
296 G4Step* G4ParticleChange::UpdateStepForAtRest( 323 G4Step* G4ParticleChange::UpdateStepForAtRest(G4Step* pStep)
297 { << 324 {
298 // A physics process always calculates the f 325 // A physics process always calculates the final state of the particle
299 326
300 G4StepPoint* pPostStepPoint = pStep->GetPost << 327 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
>> 328 G4Track* aTrack = pStep->GetTrack();
301 329
302 // set Mass/Charge << 330 // Set Mass/Charge
303 pPostStepPoint->SetMass(theMassChange); 331 pPostStepPoint->SetMass(theMassChange);
304 pPostStepPoint->SetCharge(theChargeChange); << 332 pPostStepPoint->SetCharge(theChargeChange);
305 pPostStepPoint->SetMagneticMoment(theMagneti << 333
306 <<
307 // update kinetic energy and momentum direct 334 // update kinetic energy and momentum direction
308 pPostStepPoint->SetMomentumDirection(theMome 335 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
309 pPostStepPoint->SetKineticEnergy(theEnergyCh << 336 pPostStepPoint->SetKineticEnergy( theEnergyChange );
310 if(!isVelocityChanged) <<
311 theVelocityChange = theCurrentTrack->Calcu <<
312 pPostStepPoint->SetVelocity(theVelocityChang <<
313 337
314 // update polarization 338 // update polarization
315 pPostStepPoint->SetPolarization(thePolarizat << 339 pPostStepPoint->SetPolarization( thePolarizationChange );
316 << 340
317 // update position and time 341 // update position and time
318 pPostStepPoint->SetPosition(thePositionChang << 342 pPostStepPoint->SetPosition( thePositionChange );
319 pPostStepPoint->AddGlobalTime(theTimeChange << 343 pPostStepPoint->SetGlobalTime( theTimeChange );
320 pPostStepPoint->SetLocalTime(theTimeChange); << 344 pPostStepPoint->AddLocalTime( theTimeChange
321 pPostStepPoint->SetProperTime(theProperTimeC << 345 - aTrack->GetGlobalTime());
322 << 346 pPostStepPoint->SetProperTime( theProperTimeChange );
323 if(isParentWeightProposed) << 347
324 { << 348 // update weight
325 pPostStepPoint->SetWeight(theParentWeight) << 349 pPostStepPoint->SetWeight( theWeightChange );
326 } <<
327 350
328 #ifdef G4VERBOSE 351 #ifdef G4VERBOSE
329 if(debugFlag) { CheckIt( *theCurrentTrack ); << 352 if (debugFlag) CheckIt(*aTrack);
330 #endif 353 #endif
331 354
332 // update the G4Step specific attributes << 355 // Update the G4Step specific attributes
333 return UpdateStepInfo(pStep); 356 return UpdateStepInfo(pStep);
334 } 357 }
335 358
336 // ------------------------------------------- << 359 //----------------------------------------------------------------
>> 360 // methods for printing messages
>> 361 //
>> 362
337 void G4ParticleChange::DumpInfo() const 363 void G4ParticleChange::DumpInfo() const
338 { 364 {
339 // use base-class DumpInfo << 365 // use base-class DumpInfo
340 G4VParticleChange::DumpInfo(); 366 G4VParticleChange::DumpInfo();
341 367
342 G4long oldprc = G4cout.precision(8); << 368 G4cout.precision(3);
>> 369
>> 370 G4cout << " Mass (GeV) : "
>> 371 << G4std::setw(20) << theMassChange/GeV
>> 372 << G4endl;
>> 373 G4cout << " Charge (eplus) : "
>> 374 << G4std::setw(20) << theChargeChange/eplus
>> 375 << G4endl;
>> 376 G4cout << " Position - x (mm) : "
>> 377 << G4std::setw(20) << thePositionChange.x()/mm
>> 378 << G4endl;
>> 379 G4cout << " Position - y (mm) : "
>> 380 << G4std::setw(20) << thePositionChange.y()/mm
>> 381 << G4endl;
>> 382 G4cout << " Position - z (mm) : "
>> 383 << G4std::setw(20) << thePositionChange.z()/mm
>> 384 << G4endl;
>> 385 G4cout << " Time (ns) : "
>> 386 << G4std::setw(20) << theTimeChange/ns
>> 387 << G4endl;
>> 388 G4cout << " Proper Time (ns) : "
>> 389 << G4std::setw(20) << theProperTimeChange/ns
>> 390 << G4endl;
>> 391 G4cout << " Momentum Direct - x : "
>> 392 << G4std::setw(20) << theMomentumDirectionChange.x()
>> 393 << G4endl;
>> 394 G4cout << " Momentum Direct - y : "
>> 395 << G4std::setw(20) << theMomentumDirectionChange.y()
>> 396 << G4endl;
>> 397 G4cout << " Momentum Direct - z : "
>> 398 << G4std::setw(20) << theMomentumDirectionChange.z()
>> 399 << G4endl;
>> 400 G4cout << " Kinetic Energy (MeV): "
>> 401 << G4std::setw(20) << theEnergyChange/MeV
>> 402 << G4endl;
>> 403 G4cout << " Polarization - x : "
>> 404 << G4std::setw(20) << thePolarizationChange.x()
>> 405 << G4endl;
>> 406 G4cout << " Polarization - y : "
>> 407 << G4std::setw(20) << thePolarizationChange.y()
>> 408 << G4endl;
>> 409 G4cout << " Polarization - z : "
>> 410 << G4std::setw(20) << thePolarizationChange.z()
>> 411 << G4endl;
>> 412 G4cout << " Track Weight : "
>> 413 << G4std::setw(20) << theWeightChange
>> 414 << G4endl;
>> 415 }
>> 416
>> 417 G4bool G4ParticleChange::CheckIt(const G4Track& aTrack)
>> 418 {
>> 419 G4bool exitWithError = false;
>> 420 G4double accuracy;
>> 421
>> 422 // No check in case of "fStopAndKill"
>> 423 if (GetStatusChange() == fStopAndKill ) {
>> 424 return G4VParticleChange::CheckIt(aTrack);
>> 425 }
>> 426
>> 427 // MomentumDirection should be unit vector
>> 428 G4bool itsOKforMomentum = true;
>> 429 if ( theEnergyChange >0.) {
>> 430 accuracy = abs(theMomentumDirectionChange.mag2()-1.0);
>> 431 if (accuracy > accuracyForWarning) {
>> 432 G4cout << " G4ParticleChange::CheckIt : ";
>> 433 G4cout << "the Momentum Change is not unit vector !!" << G4endl;
>> 434 G4cout << " Difference: " << accuracy << G4endl;
>> 435 itsOKforMomentum = false;
>> 436 if (accuracy > accuracyForException) exitWithError = true;
>> 437 }
>> 438 }
>> 439
>> 440 // Both global and proper time should not go back
>> 441 G4bool itsOKforGlobalTime = true;
>> 442 accuracy = (aTrack.GetGlobalTime()- theTimeChange)/ns;
>> 443 if (accuracy > accuracyForWarning) {
>> 444 G4cout << " G4ParticleChange::CheckIt : ";
>> 445 G4cout << "the global time goes back !!" << G4endl;
>> 446 G4cout << " Difference: " << accuracy << "[ns] " <<G4endl;
>> 447 itsOKforGlobalTime = false;
>> 448 if (accuracy > accuracyForException) exitWithError = true;
>> 449 }
>> 450
>> 451 G4bool itsOKforProperTime = true;
>> 452 accuracy = (aTrack.GetProperTime() - theProperTimeChange )/ns;
>> 453 if (accuracy > accuracyForWarning) {
>> 454 G4cout << " G4ParticleChange::CheckIt : ";
>> 455 G4cout << "the proper time goes back !!" << G4endl;
>> 456 G4cout << " Difference: " << accuracy << "[ns] " <<G4endl;
>> 457 itsOKforProperTime = false;
>> 458 if (accuracy > accuracyForException) exitWithError = true;
>> 459 }
>> 460
>> 461 // Kinetic Energy should not be negative
>> 462 G4bool itsOKforEnergy = true;
>> 463 accuracy = -1.0*theEnergyChange/MeV;
>> 464 if (accuracy > accuracyForWarning) {
>> 465 G4cout << " G4ParticleChange::CheckIt : ";
>> 466 G4cout << "the kinetic energy is negative !!" << G4endl;
>> 467 G4cout << " Difference: " << accuracy << "[MeV] " <<G4endl;
>> 468 itsOKforEnergy = false;
>> 469 if (accuracy > accuracyForException) exitWithError = true;
>> 470 }
>> 471
>> 472 G4bool itsOK = itsOKforMomentum && itsOKforEnergy && itsOKforProperTime && itsOKforGlobalTime;
>> 473 // dump out information of this particle change
>> 474 if (!itsOK) {
>> 475 G4cout << " G4ParticleChange::CheckIt " <<G4endl;
>> 476 DumpInfo();
>> 477 }
>> 478
>> 479 // Exit with error
>> 480 if (exitWithError) G4Exception("G4ParticleChange::CheckIt");
343 481
344 G4cout << " Mass (GeV) : " < << 482 //correction
345 << G4endl; << 483 if (!itsOKforMomentum) {
346 G4cout << " Charge (eplus) : " < << 484 G4double vmag = theMomentumDirectionChange.mag();
347 << theChargeChange / eplus << G4endl; << 485 theMomentumDirectionChange = (1./vmag)*theMomentumDirectionChange;
348 G4cout << " MagneticMoment : " < << 486 }
349 << theMagneticMomentChange << G4endl; << 487 if (!itsOKforGlobalTime) {
350 G4cout << " = : " < << 488 theTimeChange = aTrack.GetGlobalTime();
351 << theMagneticMomentChange << 489 }
352 * 2. * theMassChange / c_squared / << 490 if (!itsOKforProperTime) {
353 << "*[e hbar]/[2 m]" << G4endl; << 491 theProperTimeChange = aTrack.GetProperTime();
354 G4cout << " Position - x (mm) : " < << 492 }
355 << thePositionChange.x() / mm << G4en << 493 if (!itsOKforEnergy) {
356 G4cout << " Position - y (mm) : " < << 494 theEnergyChange = 0.0;
357 << thePositionChange.y() / mm << G4en << 495 }
358 G4cout << " Position - z (mm) : " < << 496
359 << thePositionChange.z() / mm << G4en << 497 itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
360 G4cout << " Time (ns) : " < << 498 return itsOK;
361 << theTimeChange / ns << G4endl; <<
362 G4cout << " Proper Time (ns) : " < <<
363 << theProperTimeChange / ns << G4endl <<
364 G4cout << " Momentum Direct - x : " < <<
365 << theMomentumDirectionChange.x() << <<
366 G4cout << " Momentum Direct - y : " < <<
367 << theMomentumDirectionChange.y() << <<
368 G4cout << " Momentum Direct - z : " < <<
369 << theMomentumDirectionChange.z() << <<
370 G4cout << " Kinetic Energy (MeV): " < <<
371 << theEnergyChange / MeV << G4endl; <<
372 G4cout << " Velocity (/c) : " < <<
373 << theVelocityChange / c_light << G4e <<
374 G4cout << " Polarization - x : " < <<
375 << thePolarizationChange.x() << G4end <<
376 G4cout << " Polarization - y : " < <<
377 << thePolarizationChange.y() << G4end <<
378 G4cout << " Polarization - z : " < <<
379 << thePolarizationChange.z() << G4end <<
380 G4cout.precision(oldprc); <<
381 } 499 }
>> 500
>> 501
>> 502
382 503