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