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