81 theOuterRadius(0), theNucleus(0),
82 theFieldMap(0), theEquationMap(0),
115 theFieldMap =
new std::map <G4int, G4VNuclearField*, std::less<G4int> >;
130 theEquationMap =
new std::map <G4int, G4Mag_EqRhs*, std::less<G4int> >;
233 std::vector<G4KineticTrack *>::iterator i;
234 for(i = active.begin(); i != active.end(); ++i)
240 std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter=
theFieldMap->find(encoding);
243 if ( fieldIter !=
theFieldMap->end() ) currentField=fieldIter->second;
260 #ifdef debug_1_RKPropagation
261 G4cout <<
" kt,timeStep, Intersection times tenter, tleave "
262 <<kt<<
" / state= " <<kt->
GetState() <<
" / " <<
" "<< currTimeStep <<
" / " << t_enter <<
" / " << t_leave <<
G4endl;
277 if(currTimeStep ==
DBL_MAX)currTimeStep = t_leave*1.05;
279 if ( currTimeStep >= t_leave )
294 if(t_enter > currTimeStep)
302 currTimeStep -= t_enter;
309 if(newE <= kt->GetActualMass())
347 if(currTimeStep > t_leave)
349 currTimeStep = t_leave;
353 #ifdef debug_1_RKPropagation
354 G4cerr <<
"RKPropagation is_exiting?, t_leave, curTimeStep " <<is_exiting<<
" "<<t_leave <<
" " <<currTimeStep<<
G4endl;
355 G4cout <<
"RKPropagation Ekin, field, projectile potential, p "
368 if (currTimeStep > 0 &&
373 #ifdef debug_1_RKPropagation
374 G4cout <<
"RKPropagation Ekin, field, p "
394 if(t_in < 0 && t_out >= 0)
408 if ( std::abs(currentField->
GetField(savePos)) > 0. &&
413 newE += currentField->
GetField(savePos)
419 if(newE < kt->GetActualMass())
421 #ifdef debug_1_RKPropagation
422 G4cout <<
"RKPropagation-Transport: problem with particle exiting - ignored" <<
G4endl;
443 if(newE < kt->GetActualMass())
445 #ifdef debug_1_RKPropagation
506 #ifdef debug_1_RKPropagation
507 std::cerr <<
"G4RKPropagation::FieldTransport() warning: integration error."
509 <<
G4endl <<
" timestep " <<timeStep
591 G4double scalarProd = currentPos.dot(speed);
593 G4double sqrtArg = scalarProd*scalarProd -
594 speedMag2*(currentPos.mag2()-radius*radius);
600 t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
601 t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
613 G4double sqrtArg = scalarProd*scalarProd -
614 speedMag2*(kt->
GetPosition().mag2()-radius*radius);
619 t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
620 t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
633 std::map <G4int, G4VNuclearField *, std::less<G4int> >::iterator cur;
634 for(cur = aMap->begin(); cur != aMap->end(); ++cur)
635 delete (*cur).second;
650 std::map <G4int, G4Mag_EqRhs *, std::less<G4int> >::iterator cur;
651 for(cur = aMap->begin(); cur != aMap->end(); ++cur)
652 delete (*cur).second;
static G4KaonZero * KaonZero()
G4bool FieldTransport(G4KineticTrack *track, const G4double timestep)
CascadeState GetState() const
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepLorentzRotation G4LorentzRotation
const G4ThreeVector & GetPosition() const
G4ThreeVector GetMomentumTransfer() const
virtual G4double GetField(const G4ThreeVector &aPosition)=0
G4int GetPDGEncoding() const
G4double GetActualMass() const
static const G4double eps
G4KM_DummyField * theField
void delete_FieldsAndMap(std::map< G4int, G4VNuclearField *, std::less< G4int > > *aMap)
static G4SigmaZero * SigmaZero()
virtual ~G4RKPropagation()
static G4KaonMinus * KaonMinus()
std::map< G4int, G4Mag_EqRhs *, std::less< G4int > > * theEquationMap
virtual void Init(G4V3DNucleus *nucleus)
virtual G4double GetOuterRadius()=0
G4GLOB_DLL std::ostream G4cout
CascadeState SetState(const CascadeState new_state)
G4V3DNucleus * theNucleus
static G4AntiProton * AntiProton()
void delete_EquationsAndMap(std::map< G4int, G4Mag_EqRhs *, std::less< G4int > > *aMap)
void SetPosition(const G4ThreeVector aPosition)
G4ThreeVector theMomentumTranfer
static G4Proton * Proton()
static G4PionPlus * PionPlus()
static const double second
void SetMass(G4double aMass)
static G4Neutron * Neutron()
static G4PionZero * PionZero()
static G4SigmaMinus * SigmaMinus()
G4double GetPDGMass() const
const G4LorentzVector & GetTrackingMomentum() const
static G4PionMinus * PionMinus()
void SetFactor(G4double mass, G4double opticalParameter)
void SetTrackingMomentum(const G4LorentzVector &a4Momentum)
virtual G4double GetMass()=0
static G4SigmaPlus * SigmaPlus()
static G4KaonPlus * KaonPlus()
G4bool GetSphereIntersectionTimes(const G4KineticTrack *track, G4double &t1, G4double &t2)
virtual void Transport(G4KineticTrackVector &theActive, const G4KineticTrackVector &theSpectators, G4double theTimeStep)
std::map< G4int, G4VNuclearField *, std::less< G4int > > * theFieldMap
const G4ParticleDefinition * GetDefinition() const
G4bool FreeTransport(G4KineticTrack *track, const G4double timestep)
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
static const double fermi
G4GLOB_DLL std::ostream G4cerr
G4double GetProjectilePotential() const
CLHEP::HepLorentzVector G4LorentzVector