source: trunk/source/processes/parameterisation/src/G4FastStep.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

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