ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FastStep.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4FastStep.cc
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 "G4SystemOfUnits.hh"
47 #include "G4UnitsTable.hh"
48 #include "G4Track.hh"
49 #include "G4Step.hh"
50 #include "G4TrackFastVector.hh"
51 #include "G4DynamicParticle.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();
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:
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 //----------------------------------------
88 {
91 }
92 
93 //--------------------
94 //
95 //--------------------
96 void
99  G4bool localCoordinates)
100 {
101  // Compute the position coordinate in global
102  // reference system if needed ...
103  G4ThreeVector globalPosition = position;
104  if (localCoordinates)
105  globalPosition = fFastTrack->GetInverseAffineTransformation()->
106  TransformPoint(position);
107  // ...and feed the globalPosition:
108  thePositionChange = globalPosition;
109 }
110 
111 void
114  G4bool localCoordinates)
115 {
116  ProposePrimaryTrackFinalPosition(position, localCoordinates);
117 }
118 
119 //--------------------
120 //
121 //--------------------
122 void
125  G4bool localCoordinates)
126 {
127  // Compute the momentum in global reference
128  // system if needed ...
129  G4ThreeVector globalMomentum = momentum;
130  if (localCoordinates)
131  globalMomentum = fFastTrack->GetInverseAffineTransformation()->
132  TransformAxis(momentum);
133  // ...and feed the globalMomentum (ensuring unitarity)
134  SetMomentumChange(globalMomentum.unit());
135 }
136 
137 void
140  G4bool localCoordinates)
141 {
142  ProposePrimaryTrackFinalMomentumDirection(momentum, localCoordinates);
143 }
144 
145 //--------------------
146 //
147 //--------------------
148 void
151  const G4ThreeVector &direction,
152  G4bool localCoordinates)
153 {
154  // Compute global direction if needed...
155  G4ThreeVector globalDirection = direction;
156  if (localCoordinates)
157  globalDirection =fFastTrack->GetInverseAffineTransformation()->
158  TransformAxis(direction);
159  // ...and feed the globalMomentum (ensuring unitarity)
160  SetMomentumChange(globalDirection.unit());
161  SetPrimaryTrackFinalKineticEnergy(kineticEnergy);
162 }
163 
164 void
167  const G4ThreeVector &direction,
168  G4bool localCoordinates)
169 {
170  ProposePrimaryTrackFinalKineticEnergyAndDirection(kineticEnergy, direction, localCoordinates);
171 }
172 
173 //--------------------
174 //
175 //--------------------
176 void
179  G4bool localCoordinates)
180 {
181  // Compute polarization in global system if needed:
182  G4ThreeVector globalPolarization(polarization);
183  if (localCoordinates)
184  globalPolarization = fFastTrack->GetInverseAffineTransformation()->
185  TransformAxis(globalPolarization);
186  // Feed the particle globalPolarization:
187  thePolarizationChange = globalPolarization;
188 }
189 
190 void
193  G4bool localCoordinates)
194 {
195  ProposePrimaryTrackFinalPolarization(polarization, localCoordinates);
196 }
197 
198 //--------------------
199 //
200 //--------------------
203  G4ThreeVector polarization,
205  G4double time,
206  G4bool localCoordinates )
207 {
208  G4DynamicParticle dummyDynamics(dynamics);
209 
210  // ------------------------------------------
211  // Add the polarization to the dummyDynamics:
212  // ------------------------------------------
213  dummyDynamics.SetPolarization(polarization.x(),
214  polarization.y(),
215  polarization.z());
216 
217  return CreateSecondaryTrack(dummyDynamics, position, time, localCoordinates);
218 }
219 
220 //--------------------
221 //
222 //--------------------
226  G4double time,
227  G4bool localCoordinates )
228 {
229  // ----------------------------------------
230  // Quantities in global coordinates system.
231  //
232  // The allocated globalDynamics is deleted
233  // by the destructor of the G4Track.
234  // ----------------------------------------
235  G4DynamicParticle* globalDynamics =
236  new G4DynamicParticle(dynamics);
237  G4ThreeVector globalPosition(position);
238 
239  // -----------------------------------
240  // Convert to global system if needed:
241  // -----------------------------------
242  if (localCoordinates)
243  {
244  // -- Momentum Direction:
245  globalDynamics->SetMomentumDirection(fFastTrack->
246  GetInverseAffineTransformation()->
247  TransformAxis(globalDynamics->
248  GetMomentumDirection()));
249  // -- Polarization:
250  G4ThreeVector globalPolarization;
251  globalPolarization = fFastTrack->GetInverseAffineTransformation()->
252  TransformAxis(globalDynamics->GetPolarization());
253  globalDynamics->SetPolarization(
254  globalPolarization.x(),
255  globalPolarization.y(),
256  globalPolarization.z()
257  );
258 
259  // -- Position:
260  globalPosition = fFastTrack->GetInverseAffineTransformation()->
261  TransformPoint(globalPosition);
262  }
263 
264  //-------------------------------------
265  // Create the G4Track of the secondary:
266  //-------------------------------------
267  G4Track* secondary = new G4Track(
268  globalDynamics,
269  time,
270  globalPosition
271  );
272 
273  //-------------------------------
274  // and feed the changes:
275  //-------------------------------
276  AddSecondary(secondary);
277 
278  //--------------------------------------
279  // returns the pointer on the secondary:
280  //--------------------------------------
281  return secondary;
282 }
283 
284 // G4FastStep should never be Initialized in this way
285 // but we must define it to avoid warnings.
287 {
288  G4ExceptionDescription tellWhatIsWrong;
289  tellWhatIsWrong << "G4FastStep can be initialised only through G4FastTrack."
290  << G4endl;
291  G4Exception("G4FastStep::Initialize(const G4Track&)",
292  "FastSim005",
294  tellWhatIsWrong);
295 }
296 
298  : G4VParticleChange(),
299  theEnergyChange ( 0.0 ),
300  theTimeChange ( 0.0 ),
301  theProperTimeChange( 0.0 ),
302  fFastTrack ( nullptr ),
303  theWeightChange ( 0.0 )
304 {
305  if (verboseLevel>2)
306  {
307  G4cerr << "G4FastStep::G4FastStep()" << G4endl;
308  }
309 }
310 
312 {
313  if (verboseLevel>2)
314  {
315  G4cerr << "G4FastStep::~G4FastStep()" << G4endl;
316  }
317 }
318 
319 // copy and assignment operators are implemented as "shallow copy"
322 {
323  *this = right;
324 }
325 
326 
328 {
329  if (this != &right)
330  {
346  fFastTrack = right.fFastTrack;
347  }
348  return *this;
349 }
350 
351 
352 
353 
354 
356 {
357  return ((G4VParticleChange *)this == (G4VParticleChange *) &right);
358 }
359 
361 {
362  return ((G4VParticleChange *)this != (G4VParticleChange *) &right);
363 }
364 
365 //----------------------------------------------------------------
366 // methods for updating G4Step
367 //
368 
370 {
371  // A physics process always calculates the final state of the particle
372 
373  // Take note that the return type of GetMomentumChange is a
374  // pointer to G4ParticleMometum. Also it is a normalized
375  // momentum vector.
376 
377  // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
378  G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
379  G4Track* aTrack = pStep->GetTrack();
380  // G4double mass = aTrack->GetDynamicParticle()->GetMass();
381 
382  // update kinetic energy and momentum direction
383  pPostStepPoint->SetMomentumDirection(theMomentumChange);
384  pPostStepPoint->SetKineticEnergy( theEnergyChange );
385 
386  // update polarization
387  pPostStepPoint->SetPolarization( thePolarizationChange );
388 
389  // update position and time
390  pPostStepPoint->SetPosition( thePositionChange );
391  pPostStepPoint->SetGlobalTime( theTimeChange );
392  pPostStepPoint->AddLocalTime( theTimeChange
393  - aTrack->GetGlobalTime());
394  pPostStepPoint->SetProperTime( theProperTimeChange );
395 
396  // update weight
397  pPostStepPoint->SetWeight( theWeightChange );
398 
399  if (debugFlag) CheckIt(*aTrack);
400 
401 
402  // Update the G4Step specific attributes
403  return UpdateStepInfo(pStep);
404 }
405 
407 {
408  // A physics process always calculates the final state of the particle
409 
410  // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
411  G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
412  G4Track* aTrack = pStep->GetTrack();
413  // G4double mass = aTrack->GetDynamicParticle()->GetMass();
414 
415  // update kinetic energy and momentum direction
416  pPostStepPoint->SetMomentumDirection(theMomentumChange);
417  pPostStepPoint->SetKineticEnergy( theEnergyChange );
418 
419  // update polarization
420  pPostStepPoint->SetPolarization( thePolarizationChange );
421 
422  // update position and time
423  pPostStepPoint->SetPosition( thePositionChange );
424  pPostStepPoint->SetGlobalTime( theTimeChange );
425  pPostStepPoint->AddLocalTime( theTimeChange
426  - aTrack->GetGlobalTime());
427  pPostStepPoint->SetProperTime( theProperTimeChange );
428 
429  // update weight
430  pPostStepPoint->SetWeight( theWeightChange );
431 
432  if (debugFlag) CheckIt(*aTrack);
433 
434  // Update the G4Step specific attributes
435  return UpdateStepInfo(pStep);
436 }
437 
438 //----------------------------------------------------------------
439 // methods for printing messages
440 //
441 
443 {
444 // use base-class DumpInfo
446 
447  G4cout << " Position - x (mm) : " << G4BestUnit( thePositionChange.x(), "Length" ) << G4endl;
448  G4cout << " Position - y (mm) : " << G4BestUnit( thePositionChange.y(), "Length" ) << G4endl;
449  G4cout << " Position - z (mm) : " << G4BestUnit( thePositionChange.z(), "Length" ) << G4endl;
450  G4cout << " Time (ns) : " << G4BestUnit( theTimeChange, "Time" ) << G4endl;
451  G4cout << " Proper Time (ns) : " << G4BestUnit( theProperTimeChange, "Time" ) << G4endl;
452  G4int olprc = G4cout.precision(3);
453  G4cout << " Momentum Direct - x : " << std::setw(20) << theMomentumChange.x() << G4endl;
454  G4cout << " Momentum Direct - y : " << std::setw(20) << theMomentumChange.y() << G4endl;
455  G4cout << " Momentum Direct - z : " << std::setw(20) << theMomentumChange.z() << G4endl;
456  G4cout.precision(olprc);
457  G4cout << " Kinetic Energy (MeV): " << G4BestUnit( theEnergyChange, "Energy" ) << G4endl;
458  G4cout.precision(3);
459  G4cout << " Polarization - x : " << std::setw(20) << thePolarizationChange.x() << G4endl;
460  G4cout << " Polarization - y : " << std::setw(20) << thePolarizationChange.y() << G4endl;
461  G4cout << " Polarization - z : " << std::setw(20) << thePolarizationChange.z() << G4endl;
462  G4cout.precision(olprc);
463 }
464 
466 {
467  //
468  // In the G4FastStep::CheckIt
469  // We only check a bit
470  //
471  // If the user violates the energy,
472  // We don't care, we agree.
473  //
474  // But for theMomentumDirectionChange,
475  // We do pay attention.
476  // And if too large is its range,
477  // We issue an Exception.
478  //
479  //
480  // It means, the G4FastStep::CheckIt issues an exception only for the
481  // theMomentumDirectionChange which should be an unit vector
482  // and it corrects it because it could cause problems for the ulterior
483  // tracking.For the rest, only warning are issued.
484 
485  G4bool itsOK = true;
486  G4bool exitWithError = false;
487  G4double accuracy;
488 
489  // Energy should not be larger than the initial value
490  accuracy = ( theEnergyChange - aTrack.GetKineticEnergy())/MeV;
491  if (accuracy > GetAccuracyForWarning())
492  {
494  ed << "The energy becomes larger than the initial value, difference = " << accuracy << " MeV" << G4endl;
495  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
496  "FastSim006",
497  JustWarning, ed);
498  itsOK = false;
499  if (accuracy > GetAccuracyForException()) {exitWithError = true;}
500  }
501 
502  G4bool itsOKforMomentum = true;
503  if ( theEnergyChange >0.)
504  {
505  accuracy = std::abs(theMomentumChange.mag2()-1.0);
506  if (accuracy > GetAccuracyForWarning())
507  {
509  ed << "The Momentum Change is not a unit vector, difference = " << accuracy << G4endl;
510  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
511  "FastSim007",
512  JustWarning, ed);
513  itsOK = itsOKforMomentum = false;
514  if (accuracy > GetAccuracyForException()) {exitWithError = true;}
515  }
516  }
517 
518  accuracy = (aTrack.GetGlobalTime()- theTimeChange)/ns;
519  if (accuracy > GetAccuracyForWarning())
520  {
522  ed << "The global time is getting backward, difference = " << accuracy << " ns" << G4endl;
523  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
524  "FastSim008",
525  JustWarning, ed);
526  itsOK = false;
527  }
528 
529  accuracy = (aTrack.GetProperTime() - theProperTimeChange )/ns;
530  if (accuracy > GetAccuracyForWarning())
531  {
533  ed << "The proper time is getting backward, difference = " << accuracy << " ns" << G4endl;
534  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
535  "FastSim009",
536  JustWarning, ed);
537  itsOK = false;
538  }
539 
540  if (!itsOK)
541  {
542  G4cout << "ERROR - G4FastStep::CheckIt() " << G4endl;
543  G4cout << " Pointer : " << this << G4endl ;
544  DumpInfo();
545  }
546 
547  // Exit with error
548  if (exitWithError)
549  {
551  ed << "An inaccuracy in G4FastStep is beyond tolerance." << G4endl;
552  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
553  "FastSim010",
554  FatalException, ed);
555  }
556 
557  //correction for Momentum only.
558  if (!itsOKforMomentum) {
559  G4double vmag = theMomentumChange.mag();
561  }
562 
563  itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
564  return itsOK;
565 }