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