Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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$
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 "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();
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:
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.
286 void G4FastStep::Initialize(const G4Track&)
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 
299 {
300  if (verboseLevel>2)
301  {
302  G4cerr << "G4FastStep::G4FastStep() " << G4endl;
303  }
304 }
305 
307 {
308  if (verboseLevel>2)
309  {
310  G4cerr << "G4FastStep::~G4FastStep() " << G4endl;
311  }
312 }
313 
314 // copy and assignment operators are implemented as "shallow copy"
317 {
318  *this = right;
319 }
320 
321 
323 {
324  if (this != &right)
325  {
331  theMomentumChange = right.theMomentumChange;
332  thePolarizationChange = right.thePolarizationChange;
333  thePositionChange = right.thePositionChange;
334  theTimeChange = right.theTimeChange;
335  theEnergyChange = right.theEnergyChange;
339  theWeightChange = right.theWeightChange;
340  }
341  return *this;
342 }
343 
344 
345 
346 
347 
349 {
350  return ((G4VParticleChange *)this == (G4VParticleChange *) &right);
351 }
352 
354 {
355  return ((G4VParticleChange *)this != (G4VParticleChange *) &right);
356 }
357 
358 //----------------------------------------------------------------
359 // methods for updating G4Step
360 //
361 
363 {
364  // A physics process always calculates the final state of the particle
365 
366  // Take note that the return type of GetMomentumChange is a
367  // pointer to G4ParticleMometum. Also it is a normalized
368  // momentum vector.
369 
370  // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
371  G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
372  G4Track* aTrack = pStep->GetTrack();
373  // G4double mass = aTrack->GetDynamicParticle()->GetMass();
374 
375  // update kinetic energy and momentum direction
376  pPostStepPoint->SetMomentumDirection(theMomentumChange);
377  pPostStepPoint->SetKineticEnergy( theEnergyChange );
378 
379  // update polarization
380  pPostStepPoint->SetPolarization( thePolarizationChange );
381 
382  // update position and time
383  pPostStepPoint->SetPosition( thePositionChange );
384  pPostStepPoint->SetGlobalTime( theTimeChange );
385  pPostStepPoint->AddLocalTime( theTimeChange
386  - aTrack->GetGlobalTime());
387  pPostStepPoint->SetProperTime( theProperTimeChange );
388 
389  // update weight
390  pPostStepPoint->SetWeight( theWeightChange );
391 
392  if (debugFlag) CheckIt(*aTrack);
393 
394 
395  // Update the G4Step specific attributes
396  return UpdateStepInfo(pStep);
397 }
398 
400 {
401  // A physics process always calculates the final state of the particle
402 
403  // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
404  G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
405  G4Track* aTrack = pStep->GetTrack();
406  // G4double mass = aTrack->GetDynamicParticle()->GetMass();
407 
408  // update kinetic energy and momentum direction
409  pPostStepPoint->SetMomentumDirection(theMomentumChange);
410  pPostStepPoint->SetKineticEnergy( theEnergyChange );
411 
412  // update polarization
413  pPostStepPoint->SetPolarization( thePolarizationChange );
414 
415  // update position and time
416  pPostStepPoint->SetPosition( thePositionChange );
417  pPostStepPoint->SetGlobalTime( theTimeChange );
418  pPostStepPoint->AddLocalTime( theTimeChange
419  - aTrack->GetGlobalTime());
420  pPostStepPoint->SetProperTime( theProperTimeChange );
421 
422  // update weight
423  pPostStepPoint->SetWeight( theWeightChange );
424 
425  if (debugFlag) CheckIt(*aTrack);
426 
427  // Update the G4Step specific attributes
428  return UpdateStepInfo(pStep);
429 }
430 
431 //----------------------------------------------------------------
432 // methods for printing messages
433 //
434 
436 {
437 // use base-class DumpInfo
439 
440  G4cout.precision(3);
441  G4cout << " Position - x (mm) : "
442  << std::setw(20) << thePositionChange.x()/mm
443  << G4endl;
444  G4cout << " Position - y (mm) : "
445  << std::setw(20) << thePositionChange.y()/mm
446  << G4endl;
447  G4cout << " Position - z (mm) : "
448  << std::setw(20) << thePositionChange.z()/mm
449  << G4endl;
450  G4cout << " Time (ns) : "
451  << std::setw(20) << theTimeChange/ns
452  << G4endl;
453  G4cout << " Proper Time (ns) : "
454  << std::setw(20) << theProperTimeChange/ns
455  << G4endl;
456  G4cout << " Momentum Direct - x : "
457  << std::setw(20) << theMomentumChange.x()
458  << G4endl;
459  G4cout << " Momentum Direct - y : "
460  << std::setw(20) << theMomentumChange.y()
461  << G4endl;
462  G4cout << " Momentum Direct - z : "
463  << std::setw(20) << theMomentumChange.z()
464  << G4endl;
465  G4cout << " Kinetic Energy (MeV): "
466  << std::setw(20) << theEnergyChange/MeV
467  << G4endl;
468  G4cout << " Polarization - x : "
469  << std::setw(20) << thePolarizationChange.x()
470  << G4endl;
471  G4cout << " Polarization - y : "
472  << std::setw(20) << thePolarizationChange.y()
473  << G4endl;
474  G4cout << " Polarization - z : "
475  << std::setw(20) << thePolarizationChange.z()
476  << G4endl;
477 }
478 
480 {
481  //
482  // In the G4FastStep::CheckIt
483  // We only check a bit
484  //
485  // If the user violates the energy,
486  // We don't care, we agree.
487  //
488  // But for theMomentumDirectionChange,
489  // We do pay attention.
490  // And if too large is its range,
491  // We issue an Exception.
492  //
493  //
494  // It means, the G4FastStep::CheckIt issues an exception only for the
495  // theMomentumDirectionChange which should be an unit vector
496  // and it corrects it because it could cause problems for the ulterior
497  // tracking.For the rest, only warning are issued.
498 
499  G4bool itsOK = true;
500  G4bool exitWithError = false;
501  G4double accuracy;
502 
503  // Energy should not be larger than the initial value
504  accuracy = ( theEnergyChange - aTrack.GetKineticEnergy())/MeV;
505  if (accuracy > GetAccuracyForWarning())
506  {
508  ed << "The energy becomes larger than the initial value, difference = " << accuracy << " MeV" << G4endl;
509  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
510  "FastSim006",
511  JustWarning, ed);
512  itsOK = false;
513  if (accuracy > GetAccuracyForException()) {exitWithError = true;}
514  }
515 
516  G4bool itsOKforMomentum = true;
517  if ( theEnergyChange >0.)
518  {
519  accuracy = std::abs(theMomentumChange.mag2()-1.0);
520  if (accuracy > GetAccuracyForWarning())
521  {
523  ed << "The Momentum Change is not a unit vector, difference = " << accuracy << G4endl;
524  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
525  "FastSim007",
526  JustWarning, ed);
527  itsOK = itsOKforMomentum = false;
528  if (accuracy > GetAccuracyForException()) {exitWithError = true;}
529  }
530  }
531 
532  accuracy = (aTrack.GetGlobalTime()- theTimeChange)/ns;
533  if (accuracy > GetAccuracyForWarning())
534  {
536  ed << "The global time is getting backward, difference = " << accuracy << " ns" << G4endl;
537  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
538  "FastSim008",
539  JustWarning, ed);
540  itsOK = false;
541  }
542 
543  accuracy = (aTrack.GetProperTime() - theProperTimeChange )/ns;
544  if (accuracy > GetAccuracyForWarning())
545  {
547  ed << "The proper time is getting backward, difference = " << accuracy << " ns" << G4endl;
548  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
549  "FastSim009",
550  JustWarning, ed);
551  itsOK = false;
552  }
553 
554  if (!itsOK)
555  {
556  G4cout << "ERROR - G4FastStep::CheckIt() " << G4endl;
557  G4cout << " Pointer : " << this << G4endl ;
558  DumpInfo();
559  }
560 
561  // Exit with error
562  if (exitWithError)
563  {
565  ed << "An inaccuracy in G4FastStep is beyond tolerance." << G4endl;
566  G4Exception("G4FastStep::CheckIt(const G4Track& aTrack)",
567  "FastSim010",
568  FatalException, ed);
569  }
570 
571  //correction for Momentum only.
572  if (!itsOKforMomentum) {
573  G4double vmag = theMomentumChange.mag();
574  theMomentumChange = (1./vmag)*theMomentumChange;
575  }
576 
577  itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
578  return itsOK;
579 }