Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/track/src/G4ParticleChange.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /track/src/G4ParticleChange.cc (Version 11.3.0) and /track/src/G4ParticleChange.cc (Version 10.3.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4ParticleChange class implementation       << 
 27 //                                                 26 //
 28 // Author: Hisaya Kurashige, 23 March 1998     <<  27 // $Id: G4ParticleChange.cc 92776 2015-09-16 06:57:55Z gcosmo $
 29 // ------------------------------------------- <<  28 //
                                                   >>  29 // 
                                                   >>  30 // --------------------------------------------------------------
                                                   >>  31 //  GEANT 4 class implementation file 
                                                   >>  32 //
                                                   >>  33 //  
                                                   >>  34 //  
                                                   >>  35 // ------------------------------------------------------------
                                                   >>  36 //   Implemented for the new scheme                 23 Mar. 1998  H.Kurahige
                                                   >>  37 //   Change default debug flag to false             10 May. 1998  H.Kurahige
                                                   >>  38 //   Add Track weight                               12 Nov. 1998  H.Kurashige
                                                   >>  39 //   Activate CheckIt method for VERBOSE mode       14 Dec. 1998 H.Kurashige
                                                   >>  40 //   Modified CheckIt method for time                9 Feb. 1999 H.Kurashige
                                                   >>  41 //   Rename SetXXX methods to ProposeXXX   DynamicCharge  Oct. 2005 H.Kurashige
                                                   >>  42 //   Add get/ProposeMagneticMoment                  Mar 2007 H.Kurashige
                                                   >>  43 // --------------------------------------------------------------
 30                                                    44 
 31 #include "G4ParticleChange.hh"                     45 #include "G4ParticleChange.hh"
 32 #include "G4PhysicalConstants.hh"                  46 #include "G4PhysicalConstants.hh"
 33 #include "G4SystemOfUnits.hh"                      47 #include "G4SystemOfUnits.hh"
 34 #include "G4Track.hh"                              48 #include "G4Track.hh"
 35 #include "G4Step.hh"                               49 #include "G4Step.hh"
                                                   >>  50 #include "G4TrackFastVector.hh"
 36 #include "G4DynamicParticle.hh"                    51 #include "G4DynamicParticle.hh"
 37 #include "G4ExceptionSeverity.hh"                  52 #include "G4ExceptionSeverity.hh"
 38                                                    53 
 39 // ------------------------------------------- << 
 40 G4ParticleChange::G4ParticleChange()               54 G4ParticleChange::G4ParticleChange()
 41 {}                                             <<  55   : G4VParticleChange(),  
                                                   >>  56     theMomentumDirectionChange(),
                                                   >>  57     thePolarizationChange(),
                                                   >>  58     theEnergyChange(0.), 
                                                   >>  59     theVelocityChange(0.), isVelocityChanged(false),
                                                   >>  60     thePositionChange(),
                                                   >>  61     theGlobalTime0(0.),    theLocalTime0(0.),
                                                   >>  62     theTimeChange(0.),     theProperTimeChange(0.), 
                                                   >>  63     theMassChange(0.), theChargeChange(0.),
                                                   >>  64     theMagneticMomentChange(0.), theCurrentTrack(0)
                                                   >>  65 {
                                                   >>  66 }
                                                   >>  67 
                                                   >>  68 G4ParticleChange::~G4ParticleChange() 
                                                   >>  69 {
                                                   >>  70 #ifdef G4VERBOSE
                                                   >>  71   if (verboseLevel>2) {
                                                   >>  72     G4cout << "G4ParticleChange::~G4ParticleChange() " << G4endl;
                                                   >>  73   }
                                                   >>  74 #endif
                                                   >>  75 }
                                                   >>  76 
                                                   >>  77 // copy constructor
                                                   >>  78 G4ParticleChange::G4ParticleChange(const G4ParticleChange &right)
                                                   >>  79   : G4VParticleChange(right)
                                                   >>  80 {
                                                   >>  81    if (verboseLevel>1) {
                                                   >>  82     G4cout << "G4ParticleChange::  copy constructor is called " << G4endl;
                                                   >>  83    }
                                                   >>  84    theCurrentTrack = right.theCurrentTrack;
                                                   >>  85 
                                                   >>  86    theMomentumDirectionChange = right.theMomentumDirectionChange;
                                                   >>  87    thePolarizationChange = right.thePolarizationChange;
                                                   >>  88    thePositionChange   = right.thePositionChange;
                                                   >>  89    theGlobalTime0      = right.theGlobalTime0;
                                                   >>  90    theLocalTime0       = right.theLocalTime0;
                                                   >>  91    theTimeChange       = right.theTimeChange;
                                                   >>  92    theProperTimeChange = right.theProperTimeChange;
                                                   >>  93    theEnergyChange     = right.theEnergyChange;
                                                   >>  94    theVelocityChange   = right.theVelocityChange;
                                                   >>  95    isVelocityChanged   = true;
                                                   >>  96    theMassChange       = right.theMassChange;
                                                   >>  97    theChargeChange     = right.theChargeChange;
                                                   >>  98    theMagneticMomentChange = right.theMagneticMomentChange;
                                                   >>  99 }
                                                   >> 100 
                                                   >> 101 // assignemnt operator
                                                   >> 102 G4ParticleChange & G4ParticleChange::operator=(const G4ParticleChange &right)
                                                   >> 103 {
                                                   >> 104 #ifdef G4VERBOSE
                                                   >> 105    if (verboseLevel>1) {
                                                   >> 106     G4cout << "G4ParticleChange:: assignment operator is called " << G4endl;
                                                   >> 107    }
                                                   >> 108 #endif
                                                   >> 109    if (this != &right){
                                                   >> 110      if (theNumberOfSecondaries>0) {
                                                   >> 111 #ifdef G4VERBOSE
                                                   >> 112        if (verboseLevel>0) {
                                                   >> 113    G4cout << "G4ParticleChange: assignment operator Warning  ";
                                                   >> 114    G4cout << "theListOfSecondaries is not empty ";
                                                   >> 115        }
                                                   >> 116 #endif
                                                   >> 117        for (G4int index= 0; index<theNumberOfSecondaries; index++){
                                                   >> 118    if ( (*theListOfSecondaries)[index] ) delete (*theListOfSecondaries)[index] ;
                                                   >> 119        }
                                                   >> 120      }
                                                   >> 121      delete theListOfSecondaries; 
                                                   >> 122      
                                                   >> 123     theListOfSecondaries =  new G4TrackFastVector();
                                                   >> 124     theNumberOfSecondaries = right.theNumberOfSecondaries;
                                                   >> 125     for (G4int index = 0; index<theNumberOfSecondaries; index++){
                                                   >> 126       G4Track* newTrack =  new G4Track(*((*right.theListOfSecondaries)[index] ));
                                                   >> 127       theListOfSecondaries->SetElement(index, newTrack);          }
                                                   >> 128 
                                                   >> 129      theStatusChange = right.theStatusChange;
                                                   >> 130      theCurrentTrack = right.theCurrentTrack;
                                                   >> 131      
                                                   >> 132      theMomentumDirectionChange = right.theMomentumDirectionChange;
                                                   >> 133      thePolarizationChange = right.thePolarizationChange;
                                                   >> 134      thePositionChange = right.thePositionChange;
                                                   >> 135      theGlobalTime0      = right.theGlobalTime0;
                                                   >> 136      theLocalTime0       = right.theLocalTime0;
                                                   >> 137      theTimeChange       = right.theTimeChange;
                                                   >> 138      theProperTimeChange = right.theProperTimeChange;
                                                   >> 139      theEnergyChange     = right.theEnergyChange;
                                                   >> 140      theVelocityChange   = right.theVelocityChange;
                                                   >> 141      isVelocityChanged   = true;
                                                   >> 142      theMassChange       = right.theMassChange;
                                                   >> 143      theChargeChange     = right.theChargeChange;
                                                   >> 144      theMagneticMomentChange = right.theMagneticMomentChange;
                                                   >> 145      
                                                   >> 146      theTrueStepLength = right.theTrueStepLength;
                                                   >> 147      theLocalEnergyDeposit = right.theLocalEnergyDeposit;
                                                   >> 148      theSteppingControlFlag = right.theSteppingControlFlag;
                                                   >> 149    }
                                                   >> 150    return *this;
                                                   >> 151 }
                                                   >> 152 
                                                   >> 153 G4bool G4ParticleChange::operator==(const G4ParticleChange &right) const
                                                   >> 154 {
                                                   >> 155    return ((G4VParticleChange *)this == (G4VParticleChange *) &right);
                                                   >> 156 }
                                                   >> 157 
                                                   >> 158 G4bool G4ParticleChange::operator!=(const G4ParticleChange &right) const
                                                   >> 159 {
                                                   >> 160    return ((G4VParticleChange *)this != (G4VParticleChange *) &right);
                                                   >> 161 }
                                                   >> 162 
                                                   >> 163 
                                                   >> 164 //----------------------------------------------------------------
                                                   >> 165 // methods for handling secondaries 
                                                   >> 166 //
 42                                                   167 
 43 // ------------------------------------------- << 168 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle, 
 44 void G4ParticleChange::AddSecondary(G4DynamicP << 169             G4bool   IsGoodForTracking    )
 45                                     G4bool IsG << 
 46 {                                                 170 {
 47   // create track                              << 171   //  create track
 48   G4Track* aTrack = new G4Track(aParticle, Get    172   G4Track* aTrack = new G4Track(aParticle, GetGlobalTime(), thePositionChange);
 49                                                   173 
 50   // set IsGoodGorTrackingFlag                    174   // set IsGoodGorTrackingFlag
 51   if(IsGoodForTracking)                        << 175   if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
 52     aTrack->SetGoodForTrackingFlag();          << 
 53                                                   176 
 54   // touchable handle is copied to keep the po << 177   //   Touchable handle is copied to keep the pointer
 55   aTrack->SetTouchableHandle(theCurrentTrack->    178   aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
 56                                                << 179  
 57   // add a secondary                           << 180  //  add a secondary
 58   G4VParticleChange::AddSecondary(aTrack);        181   G4VParticleChange::AddSecondary(aTrack);
 59 }                                                 182 }
 60                                                   183 
 61 // ------------------------------------------- << 184 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle, 
 62 void G4ParticleChange::AddSecondary(G4DynamicP << 185             G4ThreeVector      newPosition,
 63                                     G4ThreeVec << 186             G4bool   IsGoodForTracking    )
 64                                     G4bool IsG << 
 65 {                                                 187 {
 66   // create track                              << 188   //  create track
 67   G4Track* aTrack = new G4Track(aParticle, Get << 189   G4Track*  aTrack = new G4Track(aParticle, GetGlobalTime(), newPosition);
 68                                                   190 
 69   // set IsGoodGorTrackingFlag                    191   // set IsGoodGorTrackingFlag
 70   if(IsGoodForTracking)                        << 192   if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
 71     aTrack->SetGoodForTrackingFlag();          << 
 72                                                   193 
 73   // touchable is a temporary object, so you c << 194   //   Touchable is a temporary object, so you cannot keep the pointer
 74   aTrack->SetTouchableHandle(nullptr);         << 195   aTrack->SetTouchableHandle((G4VTouchable*)0);
 75                                                   196 
 76   // add a secondary                           << 197   //  add a secondary
 77   G4VParticleChange::AddSecondary(aTrack);        198   G4VParticleChange::AddSecondary(aTrack);
 78 }                                                 199 }
 79                                                   200 
 80 // ------------------------------------------- << 201 void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle, 
 81 void G4ParticleChange::AddSecondary(G4DynamicP << 202             G4double           newTime,
 82                                     G4double n << 203             G4bool   IsGoodForTracking    )
 83 {                                                 204 {
 84   // create track                              << 205   //  create track
 85   G4Track* aTrack = new G4Track(aParticle, new << 206   G4Track*  aTrack = new G4Track(aParticle, newTime, thePositionChange); 
 86                                                   207 
 87   // set IsGoodGorTrackingFlag                    208   // set IsGoodGorTrackingFlag
 88   if(IsGoodForTracking)                        << 209   if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
 89     aTrack->SetGoodForTrackingFlag();          << 210  
 90                                                << 211   //   Touchable handle is copied to keep the pointer
 91   // touchable handle is copied to keep the po << 
 92   aTrack->SetTouchableHandle(theCurrentTrack->    212   aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
 93                                                   213 
 94   // add a secondary                           << 214   //  add a secondary
 95   G4VParticleChange::AddSecondary(aTrack);        215   G4VParticleChange::AddSecondary(aTrack);
 96 }                                                 216 }
 97                                                   217 
 98 // ------------------------------------------- << 
 99                                                << 
100 void G4ParticleChange::AddSecondary(G4Track* a    218 void G4ParticleChange::AddSecondary(G4Track* aTrack)
101 {                                                 219 {
102   // add a secondary                           << 220   //  add a secondary
103   G4VParticleChange::AddSecondary(aTrack);        221   G4VParticleChange::AddSecondary(aTrack);
104 }                                                 222 }
105                                                   223 
106 // ------------------------------------------- << 224 //----------------------------------------------------------------
                                                   >> 225 // functions for Initialization
                                                   >> 226 //
                                                   >> 227 
107 void G4ParticleChange::Initialize(const G4Trac    228 void G4ParticleChange::Initialize(const G4Track& track)
108 {                                                 229 {
109   // use base class's method at first             230   // use base class's method at first
110   G4VParticleChange::Initialize(track);           231   G4VParticleChange::Initialize(track);
                                                   >> 232   theCurrentTrack= &track;
111                                                   233 
112   // set Energy/Momentum etc. equal to those o    234   // set Energy/Momentum etc. equal to those of the parent particle
113   const G4DynamicParticle* pParticle = track.G << 235   const G4DynamicParticle*  pParticle = track.GetDynamicParticle();
114   theEnergyChange                    = pPartic << 236   theEnergyChange            = pParticle->GetKineticEnergy();
115   theVelocityChange                  = track.G << 237   theVelocityChange          = track.GetVelocity();
116   isVelocityChanged                  = false;  << 238   isVelocityChanged          = false;
117   theMomentumDirectionChange         = pPartic << 239   theMomentumDirectionChange = pParticle->GetMomentumDirection();
118   thePolarizationChange              = pPartic << 240   thePolarizationChange      = pParticle->GetPolarization();
119   theProperTimeChange                = pPartic << 241   theProperTimeChange        = pParticle->GetProperTime();
120                                                << 242 
121   // set mass/charge/MagneticMoment of Dynamic << 243   // Set mass/charge/MagneticMoment  of DynamicParticle
122   theMassChange           = pParticle->GetMass << 244   theMassChange = pParticle->GetMass();
123   theChargeChange         = pParticle->GetChar << 245   theChargeChange = pParticle->GetCharge();
124   theMagneticMomentChange = pParticle->GetMagn    246   theMagneticMomentChange = pParticle->GetMagneticMoment();
125                                                   247 
126   // set Position equal to those of the parent << 248   // set Position  equal to those of the parent track
127   thePositionChange = track.GetPosition();     << 249   thePositionChange      = track.GetPosition();
128                                                   250 
129   // set TimeChange equal to local time of the    251   // set TimeChange equal to local time of the parent track
130   theTimeChange = theLocalTime0 = track.GetLoc << 252   theTimeChange                = track.GetLocalTime();
                                                   >> 253 
                                                   >> 254   // set initial Local/Global time of the parent track
                                                   >> 255   theLocalTime0           = track.GetLocalTime();
                                                   >> 256   theGlobalTime0          = track.GetGlobalTime();
131                                                   257 
132   // set initial global time of the parent tra << 
133   theGlobalTime0 = track.GetGlobalTime();      << 
134 }                                                 258 }
135                                                   259 
136 // ------------------------------------------- << 260 //----------------------------------------------------------------
                                                   >> 261 // methods for updating G4Step 
                                                   >> 262 //
                                                   >> 263 
137 G4Step* G4ParticleChange::UpdateStepForAlongSt    264 G4Step* G4ParticleChange::UpdateStepForAlongStep(G4Step* pStep)
138 {                                                 265 {
139   // A physics process always calculates the f    266   // A physics process always calculates the final state of the
140   // particle relative to the initial state at    267   // particle relative to the initial state at the beginning
141   // of the Step, i.e., based on information o    268   // of the Step, i.e., based on information of G4Track (or
142   // equivalently the PreStepPoint).           << 269   // equivalently the PreStepPoint). 
143   // So, the differences (delta) between these    270   // So, the differences (delta) between these two states have to be
144   // calculated and be accumulated in PostStep << 271   // calculated and be accumulated in PostStepPoint. 
145                                                << 272   
146   // Take note that the return type of GetMome << 273   // Take note that the return type of GetMomentumDirectionChange is a
147   // pointer to G4ParticleMometum. Also it is  << 274   // pointer to G4ParticleMometum. Also it is a normalized 
148   // momentum vector                           << 275   // momentum vector.
                                                   >> 276 
                                                   >> 277   G4StepPoint* pPreStepPoint  = pStep->GetPreStepPoint(); 
                                                   >> 278   G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint(); 
                                                   >> 279   G4Track* pTrack = pStep->GetTrack();
                                                   >> 280   G4double     mass = theMassChange;
149                                                   281 
150   const G4StepPoint* pPreStepPoint = pStep->Ge << 282   // Set Mass/Charge/MagneticMoment 
151   G4StepPoint* pPostStepPoint = pStep->GetPost << 
152                                                << 
153   // set Mass/Charge/MagneticMoment            << 
154   pPostStepPoint->SetMass(theMassChange);         283   pPostStepPoint->SetMass(theMassChange);
155   pPostStepPoint->SetCharge(theChargeChange);  << 284   pPostStepPoint->SetCharge(theChargeChange);  
156   pPostStepPoint->SetMagneticMoment(theMagneti << 285   pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);  
157                                                << 286  
158   // calculate new kinetic energy                 287   // calculate new kinetic energy
159   G4double preEnergy = pPreStepPoint->GetKinet    288   G4double preEnergy = pPreStepPoint->GetKineticEnergy();
160   G4double energy =                            << 289   G4double energy = pPostStepPoint->GetKineticEnergy() 
161     pPostStepPoint->GetKineticEnergy() + (theE << 290                     + (theEnergyChange - preEnergy); 
162                                                   291 
163   // update kinetic energy and momentum direct    292   // update kinetic energy and momentum direction
164   if(energy > 0.0)                             << 293   if (energy > 0.0) {
165   {                                            << 
166     // calculate new momentum                     294     // calculate new momentum
167     G4ThreeVector pMomentum = pPostStepPoint-> << 295     G4ThreeVector pMomentum =  pPostStepPoint->GetMomentum() 
168       + (CalcMomentum(theEnergyChange, theMome << 296                 + ( CalcMomentum(theEnergyChange, theMomentumDirectionChange, mass)
169           theMassChange) - pPreStepPoint->GetM << 297               - pPreStepPoint->GetMomentum());
170     G4double tMomentum2 = pMomentum.mag2();    << 298     G4double      tMomentum = pMomentum.mag();
171     G4ThreeVector direction(1.0, 0.0, 0.0);    << 299     G4ThreeVector direction(1.0,0.0,0.0); 
172     if(tMomentum2 > 0.)                        << 300     if( tMomentum > 0. ){
173     {                                          << 301       G4double  inv_Momentum= 1.0 / tMomentum; 
174       direction = pMomentum / std::sqrt(tMomen << 302       direction= pMomentum * inv_Momentum;
175     }                                             303     }
176     pPostStepPoint->SetMomentumDirection(direc    304     pPostStepPoint->SetMomentumDirection(direction);
177     pPostStepPoint->SetKineticEnergy(energy);  << 305     pPostStepPoint->SetKineticEnergy( energy );
178                                                << 306   } else {
179     // if velocity is not set it should be rec << 
180     //                                         << 
181     if(!isVelocityChanged)                     << 
182     {                                          << 
183       if (theMassChange > 0.0)                 << 
184       {                                        << 
185   theVelocityChange = CLHEP::c_light *         << 
186     std::sqrt(energy*(energy + 2*theMassChange << 
187       }                                        << 
188       else                                     << 
189       {                                        << 
190   // zero mass particle                        << 
191   theVelocityChange = CLHEP::c_light;          << 
192   // optical photon case                       << 
193   if(theCurrentTrack->GetParticleDefinition()- << 
194   {                                            << 
195           G4Track* pTrack = pStep->GetTrack(); << 
196           G4double e = pTrack->GetKineticEnerg << 
197     pTrack->SetKineticEnergy(energy);          << 
198           theVelocityChange = pTrack->Calculat << 
199     pTrack->SetKineticEnergy(e);               << 
200   }                                            << 
201       }                                        << 
202     }                                          << 
203     pPostStepPoint->SetVelocity(theVelocityCha << 
204   }                                            << 
205   else                                         << 
206   {                                            << 
207     // stop case                                  307     // stop case
                                                   >> 308     //pPostStepPoint->SetMomentumDirection(G4ThreeVector(1., 0., 0.));
208     pPostStepPoint->SetKineticEnergy(0.0);        309     pPostStepPoint->SetKineticEnergy(0.0);
209     pPostStepPoint->SetVelocity(0.0);          << 
210   }                                               310   }
                                                   >> 311   // calculate velocity
                                                   >> 312   if (!isVelocityChanged) {
                                                   >> 313     if(energy > 0.0) {
                                                   >> 314       pTrack->SetKineticEnergy(energy);
                                                   >> 315       theVelocityChange = pTrack->CalculateVelocity();
                                                   >> 316       pTrack->SetKineticEnergy(preEnergy);
                                                   >> 317     } else if(theMassChange > 0.0) {
                                                   >> 318       theVelocityChange = 0.0;
                                                   >> 319     }
                                                   >> 320   }
                                                   >> 321   pPostStepPoint->SetVelocity(theVelocityChange);
211                                                   322 
212   // update polarization                          323   // update polarization
213   pPostStepPoint->AddPolarization(thePolarizat << 324   pPostStepPoint->AddPolarization( thePolarizationChange
214                                   pPreStepPoin << 325            - pPreStepPoint->GetPolarization());
215                                                << 326       
216   // update position and time                     327   // update position and time
217   pPostStepPoint->AddPosition(thePositionChang << 328   pPostStepPoint->AddPosition( thePositionChange 
                                                   >> 329              - pPreStepPoint->GetPosition() );
218   pPostStepPoint->AddGlobalTime(theTimeChange     330   pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
219   pPostStepPoint->AddLocalTime(theTimeChange - << 331   pPostStepPoint->AddLocalTime( theTimeChange - theLocalTime0 );         
220   pPostStepPoint->AddProperTime(theProperTimeC << 332   pPostStepPoint->AddProperTime( theProperTimeChange 
221                                 pPreStepPoint- << 333          - pPreStepPoint->GetProperTime());
222                                                   334 
223   if(isParentWeightProposed)                   << 335   if (isParentWeightProposed ){
224   {                                            << 336     pPostStepPoint->SetWeight( theParentWeight );
225     pPostStepPoint->SetWeight(theParentWeight) << 
226   }                                               337   }
227                                                   338 
228 #ifdef G4VERBOSE                                  339 #ifdef G4VERBOSE
229   if(debugFlag) { CheckIt( *theCurrentTrack ); << 340   G4Track*     aTrack  = pStep->GetTrack();
                                                   >> 341   if (debugFlag) CheckIt(*aTrack);
230 #endif                                            342 #endif
231                                                   343 
232   // update the G4Step specific attributes     << 344   //  Update the G4Step specific attributes 
233   return UpdateStepInfo(pStep);                   345   return UpdateStepInfo(pStep);
234 }                                                 346 }
235                                                   347 
236 // ------------------------------------------- << 
237 G4Step* G4ParticleChange::UpdateStepForPostSte    348 G4Step* G4ParticleChange::UpdateStepForPostStep(G4Step* pStep)
238 {                                              << 349 { 
239   // A physics process always calculates the f    350   // A physics process always calculates the final state of the particle
240                                                   351 
241   // Take note that the return type of GetMome << 352   // Take note that the return type of GetMomentumChange is a
242   // pointer to G4ParticleMometum. Also it is  << 353   // pointer to G4ParticleMometum. Also it is a normalized 
243   // momentum vector                           << 354   // momentum vector.
244                                                   355 
245   G4StepPoint* pPostStepPoint = pStep->GetPost << 356   G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint(); 
246   G4Track* pTrack             = pStep->GetTrac << 357   G4Track* pTrack = pStep->GetTrack();
247                                                   358 
248   // set Mass/Charge                           << 359   // Set Mass/Charge
249   pPostStepPoint->SetMass(theMassChange);         360   pPostStepPoint->SetMass(theMassChange);
250   pPostStepPoint->SetCharge(theChargeChange);  << 361   pPostStepPoint->SetCharge(theChargeChange);  
251   pPostStepPoint->SetMagneticMoment(theMagneti << 362   pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);  
252                                                << 363  
253   // update kinetic energy and momentum direct    364   // update kinetic energy and momentum direction
254   pPostStepPoint->SetMomentumDirection(theMome    365   pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
                                                   >> 366   pPostStepPoint->SetKineticEnergy( theEnergyChange );
255                                                   367 
256   // calculate velocity                           368   // calculate velocity
257   if(theEnergyChange > 0.0)                    << 369   pTrack->SetKineticEnergy( theEnergyChange );
258   {                                            << 370   if (!isVelocityChanged) {
259     pPostStepPoint->SetKineticEnergy(theEnergy << 371     if(theEnergyChange > 0.0) {
260     pTrack->SetKineticEnergy(theEnergyChange); << 
261     if(!isVelocityChanged)                     << 
262     {                                          << 
263       theVelocityChange = pTrack->CalculateVel    372       theVelocityChange = pTrack->CalculateVelocity();
                                                   >> 373     } else if(theMassChange > 0.0) {
                                                   >> 374       theVelocityChange = 0.0;
264     }                                             375     }
265     pPostStepPoint->SetVelocity(theVelocityCha << 
266   }                                               376   }
267   else                                         << 377   pPostStepPoint->SetVelocity(theVelocityChange);
268   {                                            << 378  
269     pPostStepPoint->SetKineticEnergy(0.0);     << 379    // update polarization
270     pPostStepPoint->SetVelocity(0.0);          << 380   pPostStepPoint->SetPolarization( thePolarizationChange );
271   }                                            << 381       
272                                                << 
273   // update polarization                       << 
274   pPostStepPoint->SetPolarization(thePolarizat << 
275                                                << 
276   // update position and time                     382   // update position and time
277   pPostStepPoint->SetPosition(thePositionChang << 383   pPostStepPoint->SetPosition( thePositionChange  );
278   pPostStepPoint->AddGlobalTime(theTimeChange     384   pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
279   pPostStepPoint->SetLocalTime(theTimeChange); << 385   pPostStepPoint->SetLocalTime( theTimeChange );         
280   pPostStepPoint->SetProperTime(theProperTimeC << 386   pPostStepPoint->SetProperTime( theProperTimeChange  );
281                                                   387 
282   if(isParentWeightProposed)                   << 388   if (isParentWeightProposed ){
283   {                                            << 389     pPostStepPoint->SetWeight( theParentWeight );
284     pPostStepPoint->SetWeight(theParentWeight) << 
285   }                                               390   }
286                                                   391 
287 #ifdef G4VERBOSE                                  392 #ifdef G4VERBOSE
288   if(debugFlag) { CheckIt( *theCurrentTrack ); << 393   G4Track*     aTrack  = pStep->GetTrack();
                                                   >> 394   if (debugFlag) CheckIt(*aTrack);
289 #endif                                            395 #endif
290                                                   396 
291   // update the G4Step specific attributes     << 397   //  Update the G4Step specific attributes 
292   return UpdateStepInfo(pStep);                   398   return UpdateStepInfo(pStep);
293 }                                                 399 }
294                                                   400 
295 // ------------------------------------------- << 401 
296 G4Step* G4ParticleChange::UpdateStepForAtRest(    402 G4Step* G4ParticleChange::UpdateStepForAtRest(G4Step* pStep)
297 {                                              << 403 { 
298   // A physics process always calculates the f    404   // A physics process always calculates the final state of the particle
299                                                   405 
300   G4StepPoint* pPostStepPoint = pStep->GetPost << 406   G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint(); 
301                                                   407 
302   // set Mass/Charge                           << 408   // Set Mass/Charge
303   pPostStepPoint->SetMass(theMassChange);         409   pPostStepPoint->SetMass(theMassChange);
304   pPostStepPoint->SetCharge(theChargeChange);  << 410   pPostStepPoint->SetCharge(theChargeChange);  
305   pPostStepPoint->SetMagneticMoment(theMagneti << 411   pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);  
306                                                << 412  
307   // update kinetic energy and momentum direct    413   // update kinetic energy and momentum direction
308   pPostStepPoint->SetMomentumDirection(theMome    414   pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
309   pPostStepPoint->SetKineticEnergy(theEnergyCh << 415   pPostStepPoint->SetKineticEnergy( theEnergyChange );
310   if(!isVelocityChanged)                       << 416   if (!isVelocityChanged) theVelocityChange = pStep->GetTrack()->CalculateVelocity();
311     theVelocityChange = theCurrentTrack->Calcu << 
312   pPostStepPoint->SetVelocity(theVelocityChang    417   pPostStepPoint->SetVelocity(theVelocityChange);
313                                                   418 
314   // update polarization                          419   // update polarization
315   pPostStepPoint->SetPolarization(thePolarizat << 420   pPostStepPoint->SetPolarization( thePolarizationChange );
316                                                << 421       
317   // update position and time                     422   // update position and time
318   pPostStepPoint->SetPosition(thePositionChang << 423   pPostStepPoint->SetPosition( thePositionChange  );
319   pPostStepPoint->AddGlobalTime(theTimeChange     424   pPostStepPoint->AddGlobalTime(theTimeChange - theLocalTime0);
320   pPostStepPoint->SetLocalTime(theTimeChange); << 425   pPostStepPoint->SetLocalTime( theTimeChange );         
321   pPostStepPoint->SetProperTime(theProperTimeC << 426   pPostStepPoint->SetProperTime( theProperTimeChange  );
322                                                   427 
323   if(isParentWeightProposed)                   << 428   if (isParentWeightProposed ){
324   {                                            << 429     pPostStepPoint->SetWeight( theParentWeight );
325     pPostStepPoint->SetWeight(theParentWeight) << 
326   }                                               430   }
327                                                   431 
328 #ifdef G4VERBOSE                                  432 #ifdef G4VERBOSE
329   if(debugFlag) { CheckIt( *theCurrentTrack ); << 433   G4Track*     aTrack  = pStep->GetTrack();
                                                   >> 434   if (debugFlag) CheckIt(*aTrack);
330 #endif                                            435 #endif
331                                                   436 
332   // update the G4Step specific attributes     << 437   //  Update the G4Step specific attributes 
333   return UpdateStepInfo(pStep);                   438   return UpdateStepInfo(pStep);
334 }                                                 439 }
335                                                   440 
336 // ------------------------------------------- << 441 //----------------------------------------------------------------
                                                   >> 442 // methods for printing messages  
                                                   >> 443 //
                                                   >> 444 
337 void G4ParticleChange::DumpInfo() const           445 void G4ParticleChange::DumpInfo() const
338 {                                                 446 {
339   // use base-class DumpInfo                   << 447 // use base-class DumpInfo
340   G4VParticleChange::DumpInfo();                  448   G4VParticleChange::DumpInfo();
341                                                   449 
342   G4long oldprc = G4cout.precision(8);         << 450   G4int oldprc = G4cout.precision(3);
343                                                   451 
344   G4cout << "        Mass (GeV)          : " < << 452   G4cout << "        Mass (GeV)   : " 
345          << G4endl;                            << 453        << std::setw(20) << theMassChange/GeV
346   G4cout << "        Charge (eplus)      : " < << 454        << G4endl; 
347          << theChargeChange / eplus << G4endl; << 455   G4cout << "        Charge (eplus)   : " 
348   G4cout << "        MagneticMoment      : " < << 456        << std::setw(20) << theChargeChange/eplus
349          << theMagneticMomentChange << G4endl; << 457        << G4endl; 
350   G4cout << "                         =  : " < << 458   G4cout << "        MagneticMoment   : " 
351          << theMagneticMomentChange            << 459          << std::setw(20) << theMagneticMomentChange << G4endl;
352             * 2. * theMassChange / c_squared / << 460   G4cout << "                :  = " << std::setw(20) 
353          << "*[e hbar]/[2 m]" << G4endl;       << 461          << theMagneticMomentChange*2.*theMassChange/c_squared/eplus/hbar_Planck 
354   G4cout << "        Position - x (mm)   : " < << 462          <<  "*[e hbar]/[2 m]" 
355          << thePositionChange.x() / mm << G4en << 463          << G4endl; 
356   G4cout << "        Position - y (mm)   : " < << 464   G4cout << "        Position - x (mm)   : " 
357          << thePositionChange.y() / mm << G4en << 465        << std::setw(20) << thePositionChange.x()/mm
358   G4cout << "        Position - z (mm)   : " < << 466        << G4endl; 
359          << thePositionChange.z() / mm << G4en << 467   G4cout << "        Position - y (mm)   : " 
360   G4cout << "        Time (ns)           : " < << 468        << std::setw(20) << thePositionChange.y()/mm
361          << theTimeChange / ns << G4endl;      << 469        << G4endl; 
362   G4cout << "        Proper Time (ns)    : " < << 470   G4cout << "        Position - z (mm)   : " 
363          << theProperTimeChange / ns << G4endl << 471        << std::setw(20) << thePositionChange.z()/mm
364   G4cout << "        Momentum Direct - x : " < << 472        << G4endl;
365          << theMomentumDirectionChange.x() <<  << 473   G4cout << "        Time (ns)           : " 
366   G4cout << "        Momentum Direct - y : " < << 474        << std::setw(20) << theTimeChange/ns
367          << theMomentumDirectionChange.y() <<  << 475        << G4endl;
368   G4cout << "        Momentum Direct - z : " < << 476   G4cout << "        Proper Time (ns)    : " 
369          << theMomentumDirectionChange.z() <<  << 477        << std::setw(20) << theProperTimeChange/ns
370   G4cout << "        Kinetic Energy (MeV): " < << 478        << G4endl;
371          << theEnergyChange / MeV << G4endl;   << 479   G4cout << "        Momentum Direct - x : " 
372   G4cout << "        Velocity  (/c)      : " < << 480        << std::setw(20) << theMomentumDirectionChange.x()
373          << theVelocityChange / c_light << G4e << 481        << G4endl;
374   G4cout << "        Polarization - x    : " < << 482   G4cout << "        Momentum Direct - y : " 
375          << thePolarizationChange.x() << G4end << 483        << std::setw(20) << theMomentumDirectionChange.y()
376   G4cout << "        Polarization - y    : " < << 484        << G4endl;
377          << thePolarizationChange.y() << G4end << 485   G4cout << "        Momentum Direct - z : " 
378   G4cout << "        Polarization - z    : " < << 486        << std::setw(20) << theMomentumDirectionChange.z()
379          << thePolarizationChange.z() << G4end << 487        << G4endl;
                                                   >> 488   G4cout << "        Kinetic Energy (MeV): " 
                                                   >> 489        << std::setw(20) << theEnergyChange/MeV
                                                   >> 490        << G4endl;
                                                   >> 491   G4cout << "        Velocity  (/c): " 
                                                   >> 492        << std::setw(20) << theVelocityChange/c_light
                                                   >> 493        << G4endl;
                                                   >> 494   G4cout << "        Polarization - x    : " 
                                                   >> 495        << std::setw(20) << thePolarizationChange.x()
                                                   >> 496        << G4endl;
                                                   >> 497   G4cout << "        Polarization - y    : " 
                                                   >> 498        << std::setw(20) << thePolarizationChange.y()
                                                   >> 499        << G4endl;
                                                   >> 500   G4cout << "        Polarization - z    : " 
                                                   >> 501        << std::setw(20) <<  thePolarizationChange.z()
                                                   >> 502        << G4endl;
380   G4cout.precision(oldprc);                       503   G4cout.precision(oldprc);
                                                   >> 504 }
                                                   >> 505 
                                                   >> 506 G4bool G4ParticleChange::CheckIt(const G4Track& aTrack)
                                                   >> 507 {
                                                   >> 508   G4bool    exitWithError = false;
                                                   >> 509   G4double  accuracy;
                                                   >> 510   static G4ThreadLocal G4int nError = 0;
                                                   >> 511 #ifdef G4VERBOSE
                                                   >> 512   const  G4int maxError = 30;
                                                   >> 513 #endif
                                                   >> 514 
                                                   >> 515   // No check in case of "fStopAndKill" 
                                                   >> 516   if (GetTrackStatus() ==   fStopAndKill )  return G4VParticleChange::CheckIt(aTrack);
                                                   >> 517 
                                                   >> 518   // MomentumDirection should be unit vector
                                                   >> 519   G4bool itsOKforMomentum = true;  
                                                   >> 520   if ( theEnergyChange >0.) {
                                                   >> 521     accuracy = std::fabs(theMomentumDirectionChange.mag2()-1.0);
                                                   >> 522     if (accuracy > accuracyForWarning) {
                                                   >> 523       itsOKforMomentum = false;
                                                   >> 524       nError += 1;
                                                   >> 525       exitWithError = exitWithError || (accuracy > accuracyForException);
                                                   >> 526 #ifdef G4VERBOSE
                                                   >> 527       if (nError < maxError) {
                                                   >> 528   G4cout << "  G4ParticleChange::CheckIt  : ";
                                                   >> 529   G4cout << "the Momentum Change is not unit vector !!" 
                                                   >> 530          << "  Difference:  " << accuracy << G4endl;
                                                   >> 531   G4cout << aTrack.GetDefinition()->GetParticleName()
                                                   >> 532          << " E=" << aTrack.GetKineticEnergy()/MeV
                                                   >> 533          << " pos=" << aTrack.GetPosition().x()/m
                                                   >> 534          << ", " << aTrack.GetPosition().y()/m
                                                   >> 535          << ", " << aTrack.GetPosition().z()/m
                                                   >> 536          <<G4endl;
                                                   >> 537       }
                                                   >> 538 #endif
                                                   >> 539     }
                                                   >> 540   }
                                                   >> 541 
                                                   >> 542   // Both global and proper time should not go back
                                                   >> 543   G4bool itsOKforGlobalTime = true;  
                                                   >> 544   accuracy = (aTrack.GetLocalTime()- theTimeChange)/ns;
                                                   >> 545   if (accuracy > accuracyForWarning) {
                                                   >> 546     itsOKforGlobalTime = false;
                                                   >> 547     nError += 1;
                                                   >> 548     exitWithError = exitWithError || (accuracy > accuracyForException);
                                                   >> 549 #ifdef G4VERBOSE
                                                   >> 550     if (nError < maxError) {
                                                   >> 551       G4cout << "  G4ParticleChange::CheckIt    : ";
                                                   >> 552       G4cout << "the local time goes back  !!" 
                                                   >> 553        << "  Difference:  " << accuracy  << "[ns] " <<G4endl;
                                                   >> 554       G4cout << aTrack.GetDefinition()->GetParticleName()
                                                   >> 555        << " E=" << aTrack.GetKineticEnergy()/MeV
                                                   >> 556        << " pos=" << aTrack.GetPosition().x()/m
                                                   >> 557        << ", " << aTrack.GetPosition().y()/m
                                                   >> 558        << ", " << aTrack.GetPosition().z()/m
                                                   >> 559        << " global time=" << aTrack.GetGlobalTime()/ns
                                                   >> 560        << " local time=" << aTrack.GetLocalTime()/ns
                                                   >> 561        << " proper time=" << aTrack.GetProperTime()/ns
                                                   >> 562        << G4endl;
                                                   >> 563     }
                                                   >> 564 #endif
                                                   >> 565   }
                                                   >> 566 
                                                   >> 567   G4bool itsOKforProperTime = true;
                                                   >> 568   accuracy = (aTrack.GetProperTime() - theProperTimeChange )/ns;
                                                   >> 569   if (accuracy > accuracyForWarning) {
                                                   >> 570     itsOKforProperTime = false;
                                                   >> 571     nError += 1;
                                                   >> 572     exitWithError = exitWithError ||  (accuracy > accuracyForException);
                                                   >> 573 #ifdef G4VERBOSE
                                                   >> 574     if (nError < maxError) {
                                                   >> 575       G4cout << "  G4ParticleChange::CheckIt    : ";
                                                   >> 576       G4cout << "the proper time goes back  !!" 
                                                   >> 577        << "  Difference:  " << accuracy  << "[ns] " <<G4endl;
                                                   >> 578       G4cout << aTrack.GetDefinition()->GetParticleName()
                                                   >> 579        << " E=" << aTrack.GetKineticEnergy()/MeV
                                                   >> 580        << " pos=" << aTrack.GetPosition().x()/m
                                                   >> 581        << ", " << aTrack.GetPosition().y()/m
                                                   >> 582        << ", " << aTrack.GetPosition().z()/m
                                                   >> 583        << " global time=" << aTrack.GetGlobalTime()/ns
                                                   >> 584        << " local time=" << aTrack.GetLocalTime()/ns
                                                   >> 585        << " proper time=" << aTrack.GetProperTime()/ns
                                                   >> 586        <<G4endl;
                                                   >> 587     }
                                                   >> 588 #endif
                                                   >> 589   }
                                                   >> 590 
                                                   >> 591   // Kinetic Energy should not be negative
                                                   >> 592   G4bool itsOKforEnergy = true;
                                                   >> 593   accuracy = -1.0*theEnergyChange/MeV;
                                                   >> 594   if (accuracy > accuracyForWarning) {
                                                   >> 595     itsOKforEnergy = false;
                                                   >> 596     nError += 1;
                                                   >> 597     exitWithError = exitWithError ||   (accuracy > accuracyForException);
                                                   >> 598 #ifdef G4VERBOSE
                                                   >> 599     if (nError < maxError) {
                                                   >> 600       G4cout << "  G4ParticleChange::CheckIt    : ";
                                                   >> 601       G4cout << "the kinetic energy is negative  !!" 
                                                   >> 602        << "  Difference:  " << accuracy  << "[MeV] " <<G4endl;
                                                   >> 603       G4cout << aTrack.GetDefinition()->GetParticleName()
                                                   >> 604        << " E=" << aTrack.GetKineticEnergy()/MeV
                                                   >> 605        << " pos=" << aTrack.GetPosition().x()/m
                                                   >> 606        << ", " << aTrack.GetPosition().y()/m
                                                   >> 607        << ", " << aTrack.GetPosition().z()/m
                                                   >> 608        <<G4endl;
                                                   >> 609     }
                                                   >> 610 #endif
                                                   >> 611   }
                                                   >> 612 
                                                   >> 613   // Velocity  should not be less than c_light
                                                   >> 614   G4bool itsOKforVelocity = true;
                                                   >> 615   if (theVelocityChange < 0.) {
                                                   >> 616     itsOKforVelocity = false;
                                                   >> 617     nError += 1;
                                                   >> 618     exitWithError = true;
                                                   >> 619 #ifdef G4VERBOSE
                                                   >> 620     if (nError < maxError) {
                                                   >> 621       G4cout << "  G4ParticleChange::CheckIt    : ";
                                                   >> 622       G4cout << "the velocity is negative  !!" 
                                                   >> 623        << "  Velocity:  " << theVelocityChange/c_light  <<G4endl;
                                                   >> 624       G4cout << aTrack.GetDefinition()->GetParticleName()
                                                   >> 625        << " E=" << aTrack.GetKineticEnergy()/MeV
                                                   >> 626        << " pos=" << aTrack.GetPosition().x()/m
                                                   >> 627        << ", " << aTrack.GetPosition().y()/m
                                                   >> 628        << ", " << aTrack.GetPosition().z()/m
                                                   >> 629        <<G4endl;
                                                   >> 630     }
                                                   >> 631 #endif
                                                   >> 632   }
                                                   >> 633 
                                                   >> 634   accuracy = theVelocityChange/c_light - 1.0;
                                                   >> 635   if (accuracy > accuracyForWarning) {
                                                   >> 636     itsOKforVelocity = false;
                                                   >> 637     nError += 1;
                                                   >> 638     exitWithError = exitWithError ||  (accuracy > accuracyForException);
                                                   >> 639 #ifdef G4VERBOSE
                                                   >> 640     if (nError < maxError) {
                                                   >> 641       G4cout << "  G4ParticleChange::CheckIt    : ";
                                                   >> 642       G4cout << "the velocity is greater than c_light  !!" << G4endl;
                                                   >> 643       G4cout << "  Velocity:  " << theVelocityChange/c_light  <<G4endl;
                                                   >> 644       G4cout << aTrack.GetDefinition()->GetParticleName()
                                                   >> 645        << " E=" << aTrack.GetKineticEnergy()/MeV
                                                   >> 646        << " pos=" << aTrack.GetPosition().x()/m
                                                   >> 647        << ", " << aTrack.GetPosition().y()/m
                                                   >> 648        << ", " << aTrack.GetPosition().z()/m
                                                   >> 649        <<G4endl;
                                                   >> 650     }
                                                   >> 651 #endif
                                                   >> 652   }
                                                   >> 653 
                                                   >> 654   G4bool itsOK = itsOKforMomentum && itsOKforEnergy && itsOKforVelocity && itsOKforProperTime && itsOKforGlobalTime;
                                                   >> 655   // dump out information of this particle change
                                                   >> 656 #ifdef G4VERBOSE
                                                   >> 657   if (!itsOK) { 
                                                   >> 658     DumpInfo();
                                                   >> 659   }
                                                   >> 660 #endif
                                                   >> 661 
                                                   >> 662   // Exit with error
                                                   >> 663   if (exitWithError) {
                                                   >> 664     G4Exception("G4ParticleChange::CheckIt",
                                                   >> 665     "TRACK003", EventMustBeAborted,
                                                   >> 666     "momentum, energy, and/or time was illegal");
                                                   >> 667   }
                                                   >> 668   //correction
                                                   >> 669   if (!itsOKforMomentum) {
                                                   >> 670     G4double vmag = theMomentumDirectionChange.mag();
                                                   >> 671     theMomentumDirectionChange = (1./vmag)*theMomentumDirectionChange;
                                                   >> 672   }
                                                   >> 673   if (!itsOKforGlobalTime) {
                                                   >> 674     theTimeChange = aTrack.GetLocalTime();
                                                   >> 675   }
                                                   >> 676   if (!itsOKforProperTime) {
                                                   >> 677     theProperTimeChange = aTrack.GetProperTime();
                                                   >> 678   }
                                                   >> 679   if (!itsOKforEnergy) {
                                                   >> 680     theEnergyChange = 0.0;
                                                   >> 681   }
                                                   >> 682   if (!itsOKforVelocity) {
                                                   >> 683     theVelocityChange = c_light;
                                                   >> 684   }
                                                   >> 685 
                                                   >> 686   itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
                                                   >> 687   return itsOK;
381 }                                                 688 }
382                                                   689