| 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 // 26 //
>> 27 // $Id$
27 // 28 //
28 //-------------------------------------------- 29 //---------------------------------------------------------------
29 // 30 //
30 // G4FastStep.cc 31 // G4FastStep.cc
31 // 32 //
32 // Description: 33 // Description:
33 // Encapsulates a G4ParticleChange and insu 34 // Encapsulates a G4ParticleChange and insure friendly interface
34 // methods to manage the primary/secondarie << 35 // methods to manage the primary/secondaries final state for
35 // Fast Simulation Models. 36 // Fast Simulation Models.
36 // 37 //
37 // History: 38 // History:
38 // Oct 97: Verderi && MoraDeFreitas - First 39 // Oct 97: Verderi && MoraDeFreitas - First Implementation.
39 // Apr 98: MoraDeFreitas - G4FastStep becom 40 // Apr 98: MoraDeFreitas - G4FastStep becomes the G4ParticleChange
40 // for the Fast Simulatio 41 // for the Fast Simulation Process.
41 // 42 //
42 //-------------------------------------------- 43 //---------------------------------------------------------------
43 44
44 #include "G4FastStep.hh" 45 #include "G4FastStep.hh"
45 46
46 #include "G4DynamicParticle.hh" <<
47 #include "G4Step.hh" <<
48 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
>> 48 #include "G4UnitsTable.hh"
49 #include "G4Track.hh" 49 #include "G4Track.hh"
>> 50 #include "G4Step.hh"
50 #include "G4TrackFastVector.hh" 51 #include "G4TrackFastVector.hh"
51 #include "G4UnitsTable.hh" << 52 #include "G4DynamicParticle.hh"
52 53
53 void G4FastStep::Initialize(const G4FastTrack& 54 void G4FastStep::Initialize(const G4FastTrack& fastTrack)
54 { 55 {
55 // keeps the fastTrack reference 56 // keeps the fastTrack reference
56 fFastTrack = &fastTrack; << 57 fFastTrack=&fastTrack;
57 58
58 // currentTrack will be used to Initialize t 59 // currentTrack will be used to Initialize the other data members
59 const G4Track& currentTrack = *(fFastTrack-> 60 const G4Track& currentTrack = *(fFastTrack->GetPrimaryTrack());
60 61
61 // use base class's method at first 62 // use base class's method at first
62 G4VParticleChange::Initialize(currentTrack); 63 G4VParticleChange::Initialize(currentTrack);
63 64
64 // set Energy/Momentum etc. equal to those o 65 // set Energy/Momentum etc. equal to those of the parent particle
65 const G4DynamicParticle* pParticle = current << 66 const G4DynamicParticle* pParticle = currentTrack.GetDynamicParticle();
66 theEnergyChange = pParticle->GetKineticEnerg << 67 theEnergyChange = pParticle->GetKineticEnergy();
67 theMomentumChange = pParticle->GetMomentumDi << 68 theMomentumChange = pParticle->GetMomentumDirection();
68 thePolarizationChange = pParticle->GetPolari << 69 thePolarizationChange = pParticle->GetPolarization();
69 theProperTimeChange = pParticle->GetProperTi << 70 theProperTimeChange = pParticle->GetProperTime();
70 71
71 // set Position/Time etc. equal to those of 72 // set Position/Time etc. equal to those of the parent track
72 thePositionChange = currentTrack.GetPosition << 73 thePositionChange = currentTrack.GetPosition();
73 theTimeChange = currentTrack.GetGlobalTime() << 74 theTimeChange = currentTrack.GetGlobalTime();
74 75
75 // switch off stepping hit invokation by def 76 // switch off stepping hit invokation by default:
76 theSteppingControlFlag = AvoidHitInvocation; 77 theSteppingControlFlag = AvoidHitInvocation;
77 78
78 // event biasing weigth: 79 // event biasing weigth:
79 theWeightChange = currentTrack.GetWeight(); << 80 theWeightChange = currentTrack.GetWeight();
80 } << 81 }
81 82
82 //---------------------------------------- 83 //----------------------------------------
83 // -- Set the StopAndKilled signal 84 // -- Set the StopAndKilled signal
84 // -- and put kinetic energy to 0.0. in the 85 // -- and put kinetic energy to 0.0. in the
85 // -- G4ParticleChange. 86 // -- G4ParticleChange.
86 //---------------------------------------- 87 //----------------------------------------
87 void G4FastStep::KillPrimaryTrack() 88 void G4FastStep::KillPrimaryTrack()
88 { 89 {
89 ProposePrimaryTrackFinalKineticEnergy(0.); << 90 SetPrimaryTrackFinalKineticEnergy(0.) ;
90 ProposeTrackStatus(fStopAndKill); << 91 ProposeTrackStatus(fStopAndKill) ;
91 } 92 }
92 93
93 //-------------------- 94 //--------------------
94 // 95 //
95 //-------------------- 96 //--------------------
96 void G4FastStep::ProposePrimaryTrackFinalPosit << 97 void
97 << 98 G4FastStep::
>> 99 ProposePrimaryTrackFinalPosition(const G4ThreeVector &position,
>> 100 G4bool localCoordinates)
98 { 101 {
99 // Compute the position coordinate in global 102 // Compute the position coordinate in global
100 // reference system if needed ... 103 // reference system if needed ...
101 G4ThreeVector globalPosition = position; 104 G4ThreeVector globalPosition = position;
102 if (localCoordinates) << 105 if (localCoordinates)
103 globalPosition = fFastTrack->GetInverseAff << 106 globalPosition = fFastTrack->GetInverseAffineTransformation()->
>> 107 TransformPoint(position);
104 // ...and feed the globalPosition: 108 // ...and feed the globalPosition:
105 thePositionChange = globalPosition; 109 thePositionChange = globalPosition;
106 } 110 }
107 111
108 void G4FastStep::SetPrimaryTrackFinalPosition( << 112 void
109 << 113 G4FastStep::
>> 114 SetPrimaryTrackFinalPosition(const G4ThreeVector &position,
>> 115 G4bool localCoordinates)
110 { 116 {
111 ProposePrimaryTrackFinalPosition(position, l 117 ProposePrimaryTrackFinalPosition(position, localCoordinates);
112 } 118 }
113 119
114 //-------------------- 120 //--------------------
115 // 121 //
116 //-------------------- 122 //--------------------
117 void G4FastStep::ProposePrimaryTrackFinalMomen << 123 void
118 << 124 G4FastStep::
>> 125 ProposePrimaryTrackFinalMomentumDirection(const G4ThreeVector &momentum,
>> 126 G4bool localCoordinates)
119 { 127 {
120 // Compute the momentum in global reference 128 // Compute the momentum in global reference
121 // system if needed ... 129 // system if needed ...
122 G4ThreeVector globalMomentum = momentum; 130 G4ThreeVector globalMomentum = momentum;
123 if (localCoordinates) 131 if (localCoordinates)
124 globalMomentum = fFastTrack->GetInverseAff << 132 globalMomentum = fFastTrack->GetInverseAffineTransformation()->
>> 133 TransformAxis(momentum);
125 // ...and feed the globalMomentum (ensuring 134 // ...and feed the globalMomentum (ensuring unitarity)
126 SetMomentumChange(globalMomentum.unit()); 135 SetMomentumChange(globalMomentum.unit());
127 } 136 }
128 137
129 void G4FastStep::SetPrimaryTrackFinalMomentum( << 138 void
130 << 139 G4FastStep::
>> 140 SetPrimaryTrackFinalMomentum(const G4ThreeVector &momentum,
>> 141 G4bool localCoordinates)
131 { 142 {
132 ProposePrimaryTrackFinalMomentumDirection(mo 143 ProposePrimaryTrackFinalMomentumDirection(momentum, localCoordinates);
133 } 144 }
134 145
135 //-------------------- 146 //--------------------
136 // 147 //
137 //-------------------- 148 //--------------------
138 void G4FastStep::ProposePrimaryTrackFinalKinet << 149 void
139 << 150 G4FastStep::
140 << 151 ProposePrimaryTrackFinalKineticEnergyAndDirection(G4double kineticEnergy,
>> 152 const G4ThreeVector &direction,
>> 153 G4bool localCoordinates)
141 { 154 {
142 // Compute global direction if needed... 155 // Compute global direction if needed...
143 G4ThreeVector globalDirection = direction; 156 G4ThreeVector globalDirection = direction;
144 if (localCoordinates) 157 if (localCoordinates)
145 globalDirection = fFastTrack->GetInverseAf << 158 globalDirection =fFastTrack->GetInverseAffineTransformation()->
>> 159 TransformAxis(direction);
146 // ...and feed the globalMomentum (ensuring 160 // ...and feed the globalMomentum (ensuring unitarity)
147 SetMomentumChange(globalDirection.unit()); 161 SetMomentumChange(globalDirection.unit());
148 ProposePrimaryTrackFinalKineticEnergy(kineti << 162 SetPrimaryTrackFinalKineticEnergy(kineticEnergy);
149 } 163 }
150 164
151 void G4FastStep::SetPrimaryTrackFinalKineticEn << 165 void
152 << 166 G4FastStep::
153 << 167 SetPrimaryTrackFinalKineticEnergyAndDirection(G4double kineticEnergy,
>> 168 const G4ThreeVector &direction,
>> 169 G4bool localCoordinates)
154 { 170 {
155 ProposePrimaryTrackFinalKineticEnergyAndDire 171 ProposePrimaryTrackFinalKineticEnergyAndDirection(kineticEnergy, direction, localCoordinates);
156 } 172 }
157 173
158 //-------------------- 174 //--------------------
159 // 175 //
160 //-------------------- 176 //--------------------
161 void G4FastStep::ProposePrimaryTrackFinalPolar << 177 void
162 << 178 G4FastStep::
>> 179 ProposePrimaryTrackFinalPolarization(const G4ThreeVector &polarization,
>> 180 G4bool localCoordinates)
163 { 181 {
164 // Compute polarization in global system if 182 // Compute polarization in global system if needed:
165 G4ThreeVector globalPolarization(polarizatio 183 G4ThreeVector globalPolarization(polarization);
166 if (localCoordinates) 184 if (localCoordinates)
167 globalPolarization = << 185 globalPolarization = fFastTrack->GetInverseAffineTransformation()->
168 fFastTrack->GetInverseAffineTransformati << 186 TransformAxis(globalPolarization);
169 // Feed the particle globalPolarization: 187 // Feed the particle globalPolarization:
170 thePolarizationChange = globalPolarization; 188 thePolarizationChange = globalPolarization;
171 } 189 }
172 190
173 void G4FastStep::SetPrimaryTrackFinalPolarizat << 191 void
174 << 192 G4FastStep::
>> 193 SetPrimaryTrackFinalPolarization(const G4ThreeVector &polarization,
>> 194 G4bool localCoordinates)
175 { 195 {
176 ProposePrimaryTrackFinalPolarization(polariz 196 ProposePrimaryTrackFinalPolarization(polarization, localCoordinates);
177 } 197 }
178 198
179 //-------------------- 199 //--------------------
180 // 200 //
181 //-------------------- 201 //--------------------
182 G4Track* G4FastStep::CreateSecondaryTrack(cons << 202 G4Track* G4FastStep::
183 G4Th << 203 CreateSecondaryTrack(const G4DynamicParticle& dynamics,
184 G4do << 204 G4ThreeVector polarization,
>> 205 G4ThreeVector position,
>> 206 G4double time,
>> 207 G4bool localCoordinates )
185 { 208 {
186 G4DynamicParticle dummyDynamics(dynamics); 209 G4DynamicParticle dummyDynamics(dynamics);
187 << 210
188 // ----------------------------------------- 211 // ------------------------------------------
189 // Add the polarization to the dummyDynamics 212 // Add the polarization to the dummyDynamics:
190 // ----------------------------------------- 213 // ------------------------------------------
191 dummyDynamics.SetPolarization(polarization.x << 214 dummyDynamics.SetPolarization(polarization.x(),
192 << 215 polarization.y(),
>> 216 polarization.z());
>> 217
193 return CreateSecondaryTrack(dummyDynamics, p 218 return CreateSecondaryTrack(dummyDynamics, position, time, localCoordinates);
194 } 219 }
195 220
196 //-------------------- 221 //--------------------
197 // 222 //
198 //-------------------- 223 //--------------------
199 G4Track* G4FastStep::CreateSecondaryTrack(cons << 224 G4Track* G4FastStep::
200 G4do << 225 CreateSecondaryTrack(const G4DynamicParticle& dynamics,
>> 226 G4ThreeVector position,
>> 227 G4double time,
>> 228 G4bool localCoordinates )
201 { 229 {
202 // ---------------------------------------- 230 // ----------------------------------------
203 // Quantities in global coordinates system. 231 // Quantities in global coordinates system.
204 // << 232 //
205 // The allocated globalDynamics is deleted 233 // The allocated globalDynamics is deleted
206 // by the destructor of the G4Track. 234 // by the destructor of the G4Track.
207 // ---------------------------------------- 235 // ----------------------------------------
208 auto globalDynamics = new G4DynamicParticle( << 236 G4DynamicParticle* globalDynamics =
>> 237 new G4DynamicParticle(dynamics);
209 G4ThreeVector globalPosition(position); 238 G4ThreeVector globalPosition(position);
210 << 239
211 // ----------------------------------- 240 // -----------------------------------
212 // Convert to global system if needed: 241 // Convert to global system if needed:
213 // ----------------------------------- 242 // -----------------------------------
214 if (localCoordinates) { << 243 if (localCoordinates)
215 // -- Momentum Direction: << 244 {
216 globalDynamics->SetMomentumDirection( << 245 // -- Momentum Direction:
217 fFastTrack->GetInverseAffineTransformati << 246 globalDynamics->SetMomentumDirection(fFastTrack->
218 globalDynamics->GetMomentumDirection() << 247 GetInverseAffineTransformation()->
219 // -- Polarization: << 248 TransformAxis(globalDynamics->
220 G4ThreeVector globalPolarization; << 249 GetMomentumDirection()));
221 globalPolarization = fFastTrack->GetInvers << 250 // -- Polarization:
222 globalDynamics->GetPolarization()); << 251 G4ThreeVector globalPolarization;
223 globalDynamics->SetPolarization(globalPola << 252 globalPolarization = fFastTrack->GetInverseAffineTransformation()->
224 globalPola << 253 TransformAxis(globalDynamics->GetPolarization());
225 << 254 globalDynamics->SetPolarization(
226 // -- Position: << 255 globalPolarization.x(),
227 globalPosition = fFastTrack->GetInverseAff << 256 globalPolarization.y(),
228 } << 257 globalPolarization.z()
229 << 258 );
>> 259
>> 260 // -- Position:
>> 261 globalPosition = fFastTrack->GetInverseAffineTransformation()->
>> 262 TransformPoint(globalPosition);
>> 263 }
>> 264
230 //------------------------------------- 265 //-------------------------------------
231 // Create the G4Track of the secondary: 266 // Create the G4Track of the secondary:
232 //------------------------------------- 267 //-------------------------------------
233 auto secondary = new G4Track(globalDynamics, << 268 G4Track* secondary = new G4Track(
234 << 269 globalDynamics,
>> 270 time,
>> 271 globalPosition
>> 272 );
>> 273
235 //------------------------------- 274 //-------------------------------
236 // and feed the changes: 275 // and feed the changes:
237 //------------------------------- 276 //-------------------------------
238 AddSecondary(secondary); 277 AddSecondary(secondary);
239 << 278
240 //-------------------------------------- 279 //--------------------------------------
241 // returns the pointer on the secondary: 280 // returns the pointer on the secondary:
242 //-------------------------------------- 281 //--------------------------------------
243 return secondary; 282 return secondary;
244 } 283 }
245 284
246 // G4FastStep should never be Initialized in t 285 // G4FastStep should never be Initialized in this way
247 // but we must define it to avoid warnings. 286 // but we must define it to avoid warnings.
248 void G4FastStep::Initialize(const G4Track&) 287 void G4FastStep::Initialize(const G4Track&)
249 { 288 {
250 G4ExceptionDescription tellWhatIsWrong; << 289 G4ExceptionDescription tellWhatIsWrong;
251 tellWhatIsWrong << "G4FastStep can be initia << 290 tellWhatIsWrong << "G4FastStep can be initialised only through G4FastTrack."
252 G4Exception("G4FastStep::Initialize(const G4 << 291 << G4endl;
253 tellWhatIsWrong); << 292 G4Exception("G4FastStep::Initialize(const G4Track&)",
>> 293 "FastSim005",
>> 294 FatalException,
>> 295 tellWhatIsWrong);
>> 296 }
>> 297
>> 298 G4FastStep::G4FastStep()
>> 299 : G4VParticleChange(),
>> 300 theEnergyChange ( 0.0 ),
>> 301 theTimeChange ( 0.0 ),
>> 302 theProperTimeChange( 0.0 ),
>> 303 fFastTrack ( nullptr ),
>> 304 theWeightChange ( 0.0 )
>> 305 {
>> 306 if (verboseLevel>2)
>> 307 {
>> 308 G4cerr << "G4FastStep::G4FastStep()" << G4endl;
>> 309 }
>> 310 }
>> 311
>> 312 G4FastStep::~G4FastStep()
>> 313 {
>> 314 if (verboseLevel>2)
>> 315 {
>> 316 G4cerr << "G4FastStep::~G4FastStep()" << G4endl;
>> 317 }
>> 318 }
>> 319
>> 320 // copy and assignment operators are implemented as "shallow copy"
>> 321 G4FastStep::G4FastStep(const G4FastStep &right)
>> 322 : G4VParticleChange()
>> 323 {
>> 324 *this = right;
>> 325 }
>> 326
>> 327
>> 328 G4FastStep & G4FastStep::operator=(const G4FastStep &right)
>> 329 {
>> 330 if (this != &right)
>> 331 {
>> 332 G4VParticleChange::operator=(right);
>> 333 theListOfSecondaries = right.theListOfSecondaries;
>> 334 theSizeOftheListOfSecondaries = right.theSizeOftheListOfSecondaries;
>> 335 theNumberOfSecondaries = right.theNumberOfSecondaries;
>> 336 theStatusChange = right.theStatusChange;
>> 337 theMomentumChange = right.theMomentumChange;
>> 338 thePolarizationChange = right.thePolarizationChange;
>> 339 thePositionChange = right.thePositionChange;
>> 340 theProperTimeChange = right.theProperTimeChange;
>> 341 theTimeChange = right.theTimeChange;
>> 342 theEnergyChange = right.theEnergyChange;
>> 343 theTrueStepLength = right.theTrueStepLength;
>> 344 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
>> 345 theSteppingControlFlag = right.theSteppingControlFlag;
>> 346 theWeightChange = right.theWeightChange;
>> 347 fFastTrack = right.fFastTrack;
>> 348 }
>> 349 return *this;
>> 350 }
>> 351
>> 352
>> 353
>> 354
>> 355
>> 356 G4bool G4FastStep::operator==(const G4FastStep &right) const
>> 357 {
>> 358 return ((G4VParticleChange *)this == (G4VParticleChange *) &right);
>> 359 }
>> 360
>> 361 G4bool G4FastStep::operator!=(const G4FastStep &right) const
>> 362 {
>> 363 return ((G4VParticleChange *)this != (G4VParticleChange *) &right);
254 } 364 }
255 365
256 //-------------------------------------------- 366 //----------------------------------------------------------------
257 // methods for updating G4Step << 367 // methods for updating G4Step
258 // 368 //
259 369
260 G4Step* G4FastStep::UpdateStepForPostStep(G4St 370 G4Step* G4FastStep::UpdateStepForPostStep(G4Step* pStep)
261 { << 371 {
262 // A physics process always calculates the f 372 // A physics process always calculates the final state of the particle
263 373
264 // Take note that the return type of GetMome 374 // Take note that the return type of GetMomentumChange is a
265 // pointer to G4ParticleMometum. Also it is << 375 // pointer to G4ParticleMometum. Also it is a normalized
266 // momentum vector. 376 // momentum vector.
267 377
268 // G4StepPoint* pPreStepPoint = pStep->Get << 378 // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
269 G4StepPoint* pPostStepPoint = pStep->GetPost << 379 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
270 G4Track* aTrack = pStep->GetTrack(); << 380 G4Track* aTrack = pStep->GetTrack();
271 // G4double mass = aTrack->GetDynamicPa 381 // G4double mass = aTrack->GetDynamicParticle()->GetMass();
272 << 382
273 // update kinetic energy and momentum direct 383 // update kinetic energy and momentum direction
274 pPostStepPoint->SetMomentumDirection(theMome 384 pPostStepPoint->SetMomentumDirection(theMomentumChange);
275 pPostStepPoint->SetKineticEnergy(theEnergyCh << 385 pPostStepPoint->SetKineticEnergy( theEnergyChange );
276 <<
277 // update polarization <<
278 pPostStepPoint->SetPolarization(thePolarizat <<
279 386
>> 387 // update polarization
>> 388 pPostStepPoint->SetPolarization( thePolarizationChange );
>> 389
280 // update position and time 390 // update position and time
281 pPostStepPoint->SetPosition(thePositionChang << 391 pPostStepPoint->SetPosition( thePositionChange );
282 pPostStepPoint->SetGlobalTime(theTimeChange) << 392 pPostStepPoint->SetGlobalTime( theTimeChange );
283 pPostStepPoint->AddLocalTime(theTimeChange - << 393 pPostStepPoint->AddLocalTime( theTimeChange
284 pPostStepPoint->SetProperTime(theProperTimeC << 394 - aTrack->GetGlobalTime());
>> 395 pPostStepPoint->SetProperTime( theProperTimeChange );
285 396
286 // update weight 397 // update weight
287 pPostStepPoint->SetWeight(theWeightChange); << 398 pPostStepPoint->SetWeight( theWeightChange );
288 399
289 if (debugFlag) CheckIt(*aTrack); 400 if (debugFlag) CheckIt(*aTrack);
290 401
291 // Update the G4Step specific attributes << 402
>> 403 // Update the G4Step specific attributes
292 return UpdateStepInfo(pStep); 404 return UpdateStepInfo(pStep);
293 } 405 }
294 406
295 G4Step* G4FastStep::UpdateStepForAtRest(G4Step 407 G4Step* G4FastStep::UpdateStepForAtRest(G4Step* pStep)
296 { << 408 {
297 // A physics process always calculates the f 409 // A physics process always calculates the final state of the particle
298 410
299 // G4StepPoint* pPreStepPoint = pStep->GetP << 411 // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
300 G4StepPoint* pPostStepPoint = pStep->GetPost << 412 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
301 G4Track* aTrack = pStep->GetTrack(); << 413 G4Track* aTrack = pStep->GetTrack();
302 // G4double mass = aTrack->GetDynamicPar 414 // G4double mass = aTrack->GetDynamicParticle()->GetMass();
303 << 415
304 // update kinetic energy and momentum direct 416 // update kinetic energy and momentum direction
305 pPostStepPoint->SetMomentumDirection(theMome 417 pPostStepPoint->SetMomentumDirection(theMomentumChange);
306 pPostStepPoint->SetKineticEnergy(theEnergyCh << 418 pPostStepPoint->SetKineticEnergy( theEnergyChange );
307 419
308 // update polarization 420 // update polarization
309 pPostStepPoint->SetPolarization(thePolarizat << 421 pPostStepPoint->SetPolarization( thePolarizationChange );
310 << 422
311 // update position and time 423 // update position and time
312 pPostStepPoint->SetPosition(thePositionChang << 424 pPostStepPoint->SetPosition( thePositionChange );
313 pPostStepPoint->SetGlobalTime(theTimeChange) << 425 pPostStepPoint->SetGlobalTime( theTimeChange );
314 pPostStepPoint->AddLocalTime(theTimeChange - << 426 pPostStepPoint->AddLocalTime( theTimeChange
315 pPostStepPoint->SetProperTime(theProperTimeC << 427 - aTrack->GetGlobalTime());
>> 428 pPostStepPoint->SetProperTime( theProperTimeChange );
316 429
317 // update weight 430 // update weight
318 pPostStepPoint->SetWeight(theWeightChange); << 431 pPostStepPoint->SetWeight( theWeightChange );
319 432
320 if (debugFlag) CheckIt(*aTrack); 433 if (debugFlag) CheckIt(*aTrack);
321 434
322 // Update the G4Step specific attributes << 435 // Update the G4Step specific attributes
323 return UpdateStepInfo(pStep); 436 return UpdateStepInfo(pStep);
324 } 437 }
325 438
326 //-------------------------------------------- 439 //----------------------------------------------------------------
327 // methods for printing messages << 440 // methods for printing messages
328 // 441 //
329 442
330 void G4FastStep::DumpInfo() const 443 void G4FastStep::DumpInfo() const
331 { 444 {
332 // use base-class DumpInfo << 445 // use base-class DumpInfo
333 G4VParticleChange::DumpInfo(); 446 G4VParticleChange::DumpInfo();
334 << 447
335 G4cout << " Position - x (mm) : " < << 448 G4cout << " Position - x (mm) : " << G4BestUnit( thePositionChange.x(), "Length" ) << G4endl;
336 << G4endl; << 449 G4cout << " Position - y (mm) : " << G4BestUnit( thePositionChange.y(), "Length" ) << G4endl;
337 G4cout << " Position - y (mm) : " < << 450 G4cout << " Position - z (mm) : " << G4BestUnit( thePositionChange.z(), "Length" ) << G4endl;
338 << G4endl; << 451 G4cout << " Time (ns) : " << G4BestUnit( theTimeChange, "Time" ) << G4endl;
339 G4cout << " Position - z (mm) : " < << 452 G4cout << " Proper Time (ns) : " << G4BestUnit( theProperTimeChange, "Time" ) << G4endl;
340 << G4endl; << 453 G4int olprc = G4cout.precision(3);
341 G4cout << " Time (ns) : " < <<
342 G4cout << " Proper Time (ns) : " < <<
343 G4long olprc = G4cout.precision(3); <<
344 G4cout << " Momentum Direct - x : " < 454 G4cout << " Momentum Direct - x : " << std::setw(20) << theMomentumChange.x() << G4endl;
345 G4cout << " Momentum Direct - y : " < 455 G4cout << " Momentum Direct - y : " << std::setw(20) << theMomentumChange.y() << G4endl;
346 G4cout << " Momentum Direct - z : " < 456 G4cout << " Momentum Direct - z : " << std::setw(20) << theMomentumChange.z() << G4endl;
347 G4cout.precision(olprc); 457 G4cout.precision(olprc);
348 G4cout << " Kinetic Energy (MeV): " < << 458 G4cout << " Kinetic Energy (MeV): " << G4BestUnit( theEnergyChange, "Energy" ) << G4endl;
349 G4cout.precision(3); 459 G4cout.precision(3);
350 G4cout << " Polarization - x : " < << 460 G4cout << " Polarization - x : " << std::setw(20) << thePolarizationChange.x() << G4endl;
351 << G4endl; << 461 G4cout << " Polarization - y : " << std::setw(20) << thePolarizationChange.y() << G4endl;
352 G4cout << " Polarization - y : " < << 462 G4cout << " Polarization - z : " << std::setw(20) << thePolarizationChange.z() << G4endl;
353 << G4endl; <<
354 G4cout << " Polarization - z : " < <<
355 << G4endl; <<
356 G4cout.precision(olprc); 463 G4cout.precision(olprc);
357 } 464 }
358 465
359 G4bool G4FastStep::CheckIt(const G4Track& aTra 466 G4bool G4FastStep::CheckIt(const G4Track& aTrack)
360 { 467 {
361 // 468 //
362 // In the G4FastStep::CheckIt 469 // In the G4FastStep::CheckIt
363 // We only check a bit 470 // We only check a bit
364 // << 471 //
365 // If the user violates the energy, 472 // If the user violates the energy,
366 // We don't care, we agree. 473 // We don't care, we agree.
367 // 474 //
368 // But for theMomentumDirectionChange, 475 // But for theMomentumDirectionChange,
369 // We do pay attention. 476 // We do pay attention.
370 // And if too large is its range, 477 // And if too large is its range,
371 // We issue an Exception. 478 // We issue an Exception.
372 // 479 //
373 // 480 //
374 // It means, the G4FastStep::CheckIt issues 481 // It means, the G4FastStep::CheckIt issues an exception only for the
375 // theMomentumDirectionChange which should b 482 // theMomentumDirectionChange which should be an unit vector
376 // and it corrects it because it could cause 483 // and it corrects it because it could cause problems for the ulterior
377 // tracking.For the rest, only warning are i 484 // tracking.For the rest, only warning are issued.
378 485
379 G4bool itsOK = true; << 486 G4bool itsOK = true;
380 G4bool exitWithError = false; << 487 G4bool exitWithError = false;
381 G4double accuracy; << 488 G4double accuracy;
382 << 489
383 // Energy should not be larger than the init 490 // Energy should not be larger than the initial value
384 accuracy = (theEnergyChange - aTrack.GetKine << 491 accuracy = ( theEnergyChange - aTrack.GetKineticEnergy())/MeV;
385 if (accuracy > GetAccuracyForWarning()) { << 492 if (accuracy > GetAccuracyForWarning())
386 G4ExceptionDescription ed; << 493 {
387 ed << "The energy becomes larger than the << 494 G4ExceptionDescription ed;
388 << G4endl; << 495 ed << "The energy becomes larger than the initial value, difference = " << accuracy << " MeV" << G4endl;
389 G4Exception("G4FastStep::CheckIt(const G4T << 496 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
390 itsOK = false; << 497 "FastSim006",
391 if (accuracy > GetAccuracyForException()) << 498 JustWarning, ed);
392 exitWithError = true; << 499 itsOK = false;
>> 500 if (accuracy > GetAccuracyForException()) {exitWithError = true;}
393 } 501 }
394 } << 502
395 <<
396 G4bool itsOKforMomentum = true; 503 G4bool itsOKforMomentum = true;
397 if (theEnergyChange > 0.) { << 504 if ( theEnergyChange >0.)
398 accuracy = std::abs(theMomentumChange.mag2 << 505 {
399 if (accuracy > GetAccuracyForWarning()) { << 506 accuracy = std::abs(theMomentumChange.mag2()-1.0);
>> 507 if (accuracy > GetAccuracyForWarning())
>> 508 {
>> 509 G4ExceptionDescription ed;
>> 510 ed << "The Momentum Change is not a unit vector, difference = " << accuracy << G4endl;
>> 511 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
>> 512 "FastSim007",
>> 513 JustWarning, ed);
>> 514 itsOK = itsOKforMomentum = false;
>> 515 if (accuracy > GetAccuracyForException()) {exitWithError = true;}
>> 516 }
>> 517 }
>> 518
>> 519 accuracy = (aTrack.GetGlobalTime()- theTimeChange)/ns;
>> 520 if (accuracy > GetAccuracyForWarning())
>> 521 {
400 G4ExceptionDescription ed; 522 G4ExceptionDescription ed;
401 ed << "The Momentum Change is not a unit << 523 ed << "The global time is getting backward, difference = " << accuracy << " ns" << G4endl;
402 G4Exception("G4FastStep::CheckIt(const G << 524 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
403 itsOK = itsOKforMomentum = false; << 525 "FastSim008",
404 if (accuracy > GetAccuracyForException() << 526 JustWarning, ed);
405 exitWithError = true; << 527 itsOK = false;
406 } <<
407 } 528 }
408 } << 529
409 << 530 accuracy = (aTrack.GetProperTime() - theProperTimeChange )/ns;
410 accuracy = (aTrack.GetGlobalTime() - theTime << 531 if (accuracy > GetAccuracyForWarning())
411 if (accuracy > GetAccuracyForWarning()) { << 532 {
412 G4ExceptionDescription ed; << 533 G4ExceptionDescription ed;
413 ed << "The global time is getting backward << 534 ed << "The proper time is getting backward, difference = " << accuracy << " ns" << G4endl;
414 G4Exception("G4FastStep::CheckIt(const G4T << 535 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
415 itsOK = false; << 536 "FastSim009",
416 } << 537 JustWarning, ed);
417 << 538 itsOK = false;
418 accuracy = (aTrack.GetProperTime() - theProp << 539 }
419 if (accuracy > GetAccuracyForWarning()) { << 540
420 G4ExceptionDescription ed; << 541 if (!itsOK)
421 ed << "The proper time is getting backward << 542 {
422 G4Exception("G4FastStep::CheckIt(const G4T << 543 G4cout << "ERROR - G4FastStep::CheckIt() " << G4endl;
423 itsOK = false; << 544 G4cout << " Pointer : " << this << G4endl ;
424 } << 545 DumpInfo();
425 << 546 }
426 if (!itsOK) { << 547
427 G4cout << "ERROR - G4FastStep::CheckIt() " <<
428 G4cout << " Pointer : " << this << <<
429 DumpInfo(); <<
430 } <<
431 <<
432 // Exit with error 548 // Exit with error
433 if (exitWithError) { << 549 if (exitWithError)
434 G4ExceptionDescription ed; << 550 {
435 ed << "An inaccuracy in G4FastStep is beyo << 551 G4ExceptionDescription ed;
436 G4Exception("G4FastStep::CheckIt(const G4T << 552 ed << "An inaccuracy in G4FastStep is beyond tolerance." << G4endl;
437 } << 553 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
438 << 554 "FastSim010",
439 // correction for Momentum only. << 555 FatalException, ed);
>> 556 }
>> 557
>> 558 //correction for Momentum only.
440 if (!itsOKforMomentum) { 559 if (!itsOKforMomentum) {
441 G4double vmag = theMomentumChange.mag(); 560 G4double vmag = theMomentumChange.mag();
442 theMomentumChange = (1. / vmag) * theMomen << 561 theMomentumChange = (1./vmag)*theMomentumChange;
443 } 562 }
444 << 563
445 itsOK = (itsOK) && G4VParticleChange::CheckI << 564 itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
446 return itsOK; 565 return itsOK;
447 } 566 }
448 567