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