| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 //---------------------------------------------------------------
29 //
30 // G4FastStep.cc
31 //
32 // Description:
33 // Encapsulates a G4ParticleChange and insure friendly interface
34 // methods to manage the primary/secondaries final state for
35 // Fast Simulation Models.
36 //
37 // History:
38 // Oct 97: Verderi && MoraDeFreitas - First Implementation.
39 // Apr 98: MoraDeFreitas - G4FastStep becomes the G4ParticleChange
40 // for the Fast Simulation Process.
41 //
42 //---------------------------------------------------------------
43
44 #include "G4FastStep.hh"
45
46 #include "G4DynamicParticle.hh"
47 #include "G4Step.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4Track.hh"
50 #include "G4TrackFastVector.hh"
51 #include "G4UnitsTable.hh"
52
53 void G4FastStep::Initialize(const G4FastTrack& fastTrack)
54 {
55 // keeps the fastTrack reference
56 fFastTrack = &fastTrack;
57
58 // currentTrack will be used to Initialize the other data members
59 const G4Track& currentTrack = *(fFastTrack->GetPrimaryTrack());
60
61 // use base class's method at first
62 G4VParticleChange::Initialize(currentTrack);
63
64 // set Energy/Momentum etc. equal to those of the parent particle
65 const G4DynamicParticle* pParticle = currentTrack.GetDynamicParticle();
66 theEnergyChange = pParticle->GetKineticEnergy();
67 theMomentumChange = pParticle->GetMomentumDirection();
68 thePolarizationChange = pParticle->GetPolarization();
69 theProperTimeChange = pParticle->GetProperTime();
70
71 // set Position/Time etc. equal to those of the parent track
72 thePositionChange = currentTrack.GetPosition();
73 theTimeChange = currentTrack.GetGlobalTime();
74
75 // switch off stepping hit invokation by default:
76 theSteppingControlFlag = AvoidHitInvocation;
77
78 // event biasing weigth:
79 theWeightChange = currentTrack.GetWeight();
80 }
81
82 //----------------------------------------
83 // -- Set the StopAndKilled signal
84 // -- and put kinetic energy to 0.0. in the
85 // -- G4ParticleChange.
86 //----------------------------------------
87 void G4FastStep::KillPrimaryTrack()
88 {
89 ProposePrimaryTrackFinalKineticEnergy(0.);
90 ProposeTrackStatus(fStopAndKill);
91 }
92
93 //--------------------
94 //
95 //--------------------
96 void G4FastStep::ProposePrimaryTrackFinalPosition(const G4ThreeVector& position,
97 G4bool localCoordinates)
98 {
99 // Compute the position coordinate in global
100 // reference system if needed ...
101 G4ThreeVector globalPosition = position;
102 if (localCoordinates)
103 globalPosition = fFastTrack->GetInverseAffineTransformation()->TransformPoint(position);
104 // ...and feed the globalPosition:
105 thePositionChange = globalPosition;
106 }
107
108 void G4FastStep::SetPrimaryTrackFinalPosition(const G4ThreeVector& position,
109 G4bool localCoordinates)
110 {
111 ProposePrimaryTrackFinalPosition(position, localCoordinates);
112 }
113
114 //--------------------
115 //
116 //--------------------
117 void G4FastStep::ProposePrimaryTrackFinalMomentumDirection(const G4ThreeVector& momentum,
118 G4bool localCoordinates)
119 {
120 // Compute the momentum in global reference
121 // system if needed ...
122 G4ThreeVector globalMomentum = momentum;
123 if (localCoordinates)
124 globalMomentum = fFastTrack->GetInverseAffineTransformation()->TransformAxis(momentum);
125 // ...and feed the globalMomentum (ensuring unitarity)
126 SetMomentumChange(globalMomentum.unit());
127 }
128
129 void G4FastStep::SetPrimaryTrackFinalMomentum(const G4ThreeVector& momentum,
130 G4bool localCoordinates)
131 {
132 ProposePrimaryTrackFinalMomentumDirection(momentum, localCoordinates);
133 }
134
135 //--------------------
136 //
137 //--------------------
138 void G4FastStep::ProposePrimaryTrackFinalKineticEnergyAndDirection(G4double kineticEnergy,
139 const G4ThreeVector& direction,
140 G4bool localCoordinates)
141 {
142 // Compute global direction if needed...
143 G4ThreeVector globalDirection = direction;
144 if (localCoordinates)
145 globalDirection = fFastTrack->GetInverseAffineTransformation()->TransformAxis(direction);
146 // ...and feed the globalMomentum (ensuring unitarity)
147 SetMomentumChange(globalDirection.unit());
148 ProposePrimaryTrackFinalKineticEnergy(kineticEnergy);
149 }
150
151 void G4FastStep::SetPrimaryTrackFinalKineticEnergyAndDirection(G4double kineticEnergy,
152 const G4ThreeVector& direction,
153 G4bool localCoordinates)
154 {
155 ProposePrimaryTrackFinalKineticEnergyAndDirection(kineticEnergy, direction, localCoordinates);
156 }
157
158 //--------------------
159 //
160 //--------------------
161 void G4FastStep::ProposePrimaryTrackFinalPolarization(const G4ThreeVector& polarization,
162 G4bool localCoordinates)
163 {
164 // Compute polarization in global system if needed:
165 G4ThreeVector globalPolarization(polarization);
166 if (localCoordinates)
167 globalPolarization =
168 fFastTrack->GetInverseAffineTransformation()->TransformAxis(globalPolarization);
169 // Feed the particle globalPolarization:
170 thePolarizationChange = globalPolarization;
171 }
172
173 void G4FastStep::SetPrimaryTrackFinalPolarization(const G4ThreeVector& polarization,
174 G4bool localCoordinates)
175 {
176 ProposePrimaryTrackFinalPolarization(polarization, localCoordinates);
177 }
178
179 //--------------------
180 //
181 //--------------------
182 G4Track* G4FastStep::CreateSecondaryTrack(const G4DynamicParticle& dynamics,
183 G4ThreeVector polarization, G4ThreeVector position,
184 G4double time, G4bool localCoordinates)
185 {
186 G4DynamicParticle dummyDynamics(dynamics);
187
188 // ------------------------------------------
189 // Add the polarization to the dummyDynamics:
190 // ------------------------------------------
191 dummyDynamics.SetPolarization(polarization.x(), polarization.y(), polarization.z());
192
193 return CreateSecondaryTrack(dummyDynamics, position, time, localCoordinates);
194 }
195
196 //--------------------
197 //
198 //--------------------
199 G4Track* G4FastStep::CreateSecondaryTrack(const G4DynamicParticle& dynamics, G4ThreeVector position,
200 G4double time, G4bool localCoordinates)
201 {
202 // ----------------------------------------
203 // Quantities in global coordinates system.
204 //
205 // The allocated globalDynamics is deleted
206 // by the destructor of the G4Track.
207 // ----------------------------------------
208 auto globalDynamics = new G4DynamicParticle(dynamics);
209 G4ThreeVector globalPosition(position);
210
211 // -----------------------------------
212 // Convert to global system if needed:
213 // -----------------------------------
214 if (localCoordinates) {
215 // -- Momentum Direction:
216 globalDynamics->SetMomentumDirection(
217 fFastTrack->GetInverseAffineTransformation()->TransformAxis(
218 globalDynamics->GetMomentumDirection()));
219 // -- Polarization:
220 G4ThreeVector globalPolarization;
221 globalPolarization = fFastTrack->GetInverseAffineTransformation()->TransformAxis(
222 globalDynamics->GetPolarization());
223 globalDynamics->SetPolarization(globalPolarization.x(), globalPolarization.y(),
224 globalPolarization.z());
225
226 // -- Position:
227 globalPosition = fFastTrack->GetInverseAffineTransformation()->TransformPoint(globalPosition);
228 }
229
230 //-------------------------------------
231 // Create the G4Track of the secondary:
232 //-------------------------------------
233 auto secondary = new G4Track(globalDynamics, time, globalPosition);
234
235 //-------------------------------
236 // and feed the changes:
237 //-------------------------------
238 AddSecondary(secondary);
239
240 //--------------------------------------
241 // returns the pointer on the secondary:
242 //--------------------------------------
243 return secondary;
244 }
245
246 // G4FastStep should never be Initialized in this way
247 // but we must define it to avoid warnings.
248 void G4FastStep::Initialize(const G4Track&)
249 {
250 G4ExceptionDescription tellWhatIsWrong;
251 tellWhatIsWrong << "G4FastStep can be initialised only through G4FastTrack." << G4endl;
252 G4Exception("G4FastStep::Initialize(const G4Track&)", "FastSim005", FatalException,
253 tellWhatIsWrong);
254 }
255
256 //----------------------------------------------------------------
257 // methods for updating G4Step
258 //
259
260 G4Step* G4FastStep::UpdateStepForPostStep(G4Step* pStep)
261 {
262 // A physics process always calculates the final state of the particle
263
264 // Take note that the return type of GetMomentumChange is a
265 // pointer to G4ParticleMometum. Also it is a normalized
266 // momentum vector.
267
268 // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
269 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
270 G4Track* aTrack = pStep->GetTrack();
271 // G4double mass = aTrack->GetDynamicParticle()->GetMass();
272
273 // update kinetic energy and momentum direction
274 pPostStepPoint->SetMomentumDirection(theMomentumChange);
275 pPostStepPoint->SetKineticEnergy(theEnergyChange);
276
277 // update polarization
278 pPostStepPoint->SetPolarization(thePolarizationChange);
279
280 // update position and time
281 pPostStepPoint->SetPosition(thePositionChange);
282 pPostStepPoint->SetGlobalTime(theTimeChange);
283 pPostStepPoint->AddLocalTime(theTimeChange - aTrack->GetGlobalTime());
284 pPostStepPoint->SetProperTime(theProperTimeChange);
285
286 // update weight
287 pPostStepPoint->SetWeight(theWeightChange);
288
289 if (debugFlag) CheckIt(*aTrack);
290
291 // Update the G4Step specific attributes
292 return UpdateStepInfo(pStep);
293 }
294
295 G4Step* G4FastStep::UpdateStepForAtRest(G4Step* pStep)
296 {
297 // A physics process always calculates the final state of the particle
298
299 // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
300 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
301 G4Track* aTrack = pStep->GetTrack();
302 // G4double mass = aTrack->GetDynamicParticle()->GetMass();
303
304 // update kinetic energy and momentum direction
305 pPostStepPoint->SetMomentumDirection(theMomentumChange);
306 pPostStepPoint->SetKineticEnergy(theEnergyChange);
307
308 // update polarization
309 pPostStepPoint->SetPolarization(thePolarizationChange);
310
311 // update position and time
312 pPostStepPoint->SetPosition(thePositionChange);
313 pPostStepPoint->SetGlobalTime(theTimeChange);
314 pPostStepPoint->AddLocalTime(theTimeChange - aTrack->GetGlobalTime());
315 pPostStepPoint->SetProperTime(theProperTimeChange);
316
317 // update weight
318 pPostStepPoint->SetWeight(theWeightChange);
319
320 if (debugFlag) CheckIt(*aTrack);
321
322 // Update the G4Step specific attributes
323 return UpdateStepInfo(pStep);
324 }
325
326 //----------------------------------------------------------------
327 // methods for printing messages
328 //
329
330 void G4FastStep::DumpInfo() const
331 {
332 // use base-class DumpInfo
333 G4VParticleChange::DumpInfo();
334
335 G4cout << " Position - x (mm) : " << G4BestUnit(thePositionChange.x(), "Length")
336 << G4endl;
337 G4cout << " Position - y (mm) : " << G4BestUnit(thePositionChange.y(), "Length")
338 << G4endl;
339 G4cout << " Position - z (mm) : " << G4BestUnit(thePositionChange.z(), "Length")
340 << G4endl;
341 G4cout << " Time (ns) : " << G4BestUnit(theTimeChange, "Time") << G4endl;
342 G4cout << " Proper Time (ns) : " << G4BestUnit(theProperTimeChange, "Time") << G4endl;
343 G4long olprc = G4cout.precision(3);
344 G4cout << " Momentum Direct - x : " << std::setw(20) << theMomentumChange.x() << G4endl;
345 G4cout << " Momentum Direct - y : " << std::setw(20) << theMomentumChange.y() << G4endl;
346 G4cout << " Momentum Direct - z : " << std::setw(20) << theMomentumChange.z() << G4endl;
347 G4cout.precision(olprc);
348 G4cout << " Kinetic Energy (MeV): " << G4BestUnit(theEnergyChange, "Energy") << G4endl;
349 G4cout.precision(3);
350 G4cout << " Polarization - x : " << std::setw(20) << thePolarizationChange.x()
351 << G4endl;
352 G4cout << " Polarization - y : " << std::setw(20) << thePolarizationChange.y()
353 << G4endl;
354 G4cout << " Polarization - z : " << std::setw(20) << thePolarizationChange.z()
355 << G4endl;
356 G4cout.precision(olprc);
357 }
358
359 G4bool G4FastStep::CheckIt(const G4Track& aTrack)
360 {
361 //
362 // In the G4FastStep::CheckIt
363 // We only check a bit
364 //
365 // If the user violates the energy,
366 // We don't care, we agree.
367 //
368 // But for theMomentumDirectionChange,
369 // We do pay attention.
370 // And if too large is its range,
371 // We issue an Exception.
372 //
373 //
374 // It means, the G4FastStep::CheckIt issues an exception only for the
375 // theMomentumDirectionChange which should be an unit vector
376 // and it corrects it because it could cause problems for the ulterior
377 // tracking.For the rest, only warning are issued.
378
379 G4bool itsOK = true;
380 G4bool exitWithError = false;
381 G4double accuracy;
382
383 // Energy should not be larger than the initial value
384 accuracy = (theEnergyChange - aTrack.GetKineticEnergy()) / MeV;
385 if (accuracy > GetAccuracyForWarning()) {
386 G4ExceptionDescription ed;
387 ed << "The energy becomes larger than the initial value, difference = " << accuracy << " MeV"
388 << G4endl;
389 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)", "FastSim006", JustWarning, ed);
390 itsOK = false;
391 if (accuracy > GetAccuracyForException()) {
392 exitWithError = true;
393 }
394 }
395
396 G4bool itsOKforMomentum = true;
397 if (theEnergyChange > 0.) {
398 accuracy = std::abs(theMomentumChange.mag2() - 1.0);
399 if (accuracy > GetAccuracyForWarning()) {
400 G4ExceptionDescription ed;
401 ed << "The Momentum Change is not a unit vector, difference = " << accuracy << G4endl;
402 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)", "FastSim007", JustWarning, ed);
403 itsOK = itsOKforMomentum = false;
404 if (accuracy > GetAccuracyForException()) {
405 exitWithError = true;
406 }
407 }
408 }
409
410 accuracy = (aTrack.GetGlobalTime() - theTimeChange) / ns;
411 if (accuracy > GetAccuracyForWarning()) {
412 G4ExceptionDescription ed;
413 ed << "The global time is getting backward, difference = " << accuracy << " ns" << G4endl;
414 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)", "FastSim008", JustWarning, ed);
415 itsOK = false;
416 }
417
418 accuracy = (aTrack.GetProperTime() - theProperTimeChange) / ns;
419 if (accuracy > GetAccuracyForWarning()) {
420 G4ExceptionDescription ed;
421 ed << "The proper time is getting backward, difference = " << accuracy << " ns" << G4endl;
422 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)", "FastSim009", JustWarning, ed);
423 itsOK = false;
424 }
425
426 if (!itsOK) {
427 G4cout << "ERROR - G4FastStep::CheckIt() " << G4endl;
428 G4cout << " Pointer : " << this << G4endl;
429 DumpInfo();
430 }
431
432 // Exit with error
433 if (exitWithError) {
434 G4ExceptionDescription ed;
435 ed << "An inaccuracy in G4FastStep is beyond tolerance." << G4endl;
436 G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)", "FastSim010", FatalException, ed);
437 }
438
439 // correction for Momentum only.
440 if (!itsOKforMomentum) {
441 G4double vmag = theMomentumChange.mag();
442 theMomentumChange = (1. / vmag) * theMomentumChange;
443 }
444
445 itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
446 return itsOK;
447 }
448