#include <G4RKPropagation.hh>
Inheritance diagram for G4RKPropagation:
Public Member Functions | |
G4RKPropagation () | |
virtual | ~G4RKPropagation () |
virtual void | Init (G4V3DNucleus *nucleus) |
virtual void | Transport (G4KineticTrackVector &theActive, const G4KineticTrackVector &theSpectators, G4double theTimeStep) |
G4bool | GetSphereIntersectionTimes (const G4KineticTrack *track, G4double &t1, G4double &t2) |
G4ThreeVector | GetMomentumTransfer () const |
G4double | GetBarrier (G4int encoding) |
G4double | GetField (G4int encoding, G4ThreeVector pos) |
Definition at line 38 of file G4RKPropagation.hh.
G4RKPropagation::G4RKPropagation | ( | ) |
Definition at line 80 of file G4RKPropagation.cc.
00080 : 00081 theOuterRadius(0), theNucleus(0), 00082 theFieldMap(0), theEquationMap(0), 00083 theField(0) 00084 { }
G4RKPropagation::~G4RKPropagation | ( | ) | [virtual] |
Definition at line 87 of file G4RKPropagation.cc.
00088 { 00089 // free theFieldMap memory 00090 if(theFieldMap) delete_FieldsAndMap(theFieldMap); 00091 00092 // free theEquationMap memory 00093 if(theEquationMap) delete_EquationsAndMap(theEquationMap); 00094 00095 if (theField) delete theField; 00096 }
Definition at line 84 of file G4RKPropagation.hh.
00085 { 00086 std::map <G4int, G4VNuclearField *, std::less<G4int> >::iterator iter; 00087 iter = theFieldMap->find(encoding); 00088 if(iter == theFieldMap->end()) return 0; 00089 return (*theFieldMap)[encoding]->GetBarrier(); 00090 }
G4double G4RKPropagation::GetField | ( | G4int | encoding, | |
G4ThreeVector | pos | |||
) | [inline] |
Definition at line 92 of file G4RKPropagation.hh.
00093 { 00094 std::map <G4int, G4VNuclearField *, std::less<G4int> >::iterator iter; 00095 iter = theFieldMap->find(encoding); 00096 if(iter == theFieldMap->end()) return 0; 00097 return (*theFieldMap)[encoding]->GetField(pos); 00098 }
G4ThreeVector G4RKPropagation::GetMomentumTransfer | ( | ) | const [virtual] |
G4bool G4RKPropagation::GetSphereIntersectionTimes | ( | const G4KineticTrack * | track, | |
G4double & | t1, | |||
G4double & | t2 | |||
) |
Definition at line 607 of file G4RKPropagation.cc.
References G4KineticTrack::GetPosition(), and G4KineticTrack::GetTrackingMomentum().
Referenced by Transport().
00609 { 00610 G4double radius = theOuterRadius + 3*fermi; // "safety" of 3 fermi 00611 G4ThreeVector speed = kt->GetTrackingMomentum().vect()/kt->GetTrackingMomentum().e(); // bost vector 00612 G4double scalarProd = kt->GetPosition().dot(speed); 00613 G4double speedMag2 = speed.mag2(); 00614 G4double sqrtArg = scalarProd*scalarProd - 00615 speedMag2*(kt->GetPosition().mag2()-radius*radius); 00616 if(sqrtArg <= 0.) // particle will not intersect the sphere 00617 { 00618 return false; 00619 } 00620 t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light; 00621 t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light; 00622 return true; 00623 }
void G4RKPropagation::Init | ( | G4V3DNucleus * | nucleus | ) | [virtual] |
Implements G4VFieldPropagation.
Definition at line 99 of file G4RKPropagation.cc.
References G4AntiProton::AntiProton(), G4V3DNucleus::GetOuterRadius(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZero::KaonZero(), G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4Proton::Proton(), G4KM_OpticalEqRhs::SetFactor(), G4KM_NucleonEqRhs::SetMass(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), and G4SigmaZero::SigmaZero().
00101 { 00102 00103 // free theFieldMap memory 00104 if(theFieldMap) delete_FieldsAndMap(theFieldMap); 00105 00106 // free theEquationMap memory 00107 if(theEquationMap) delete_EquationsAndMap(theEquationMap); 00108 00109 if (theField) delete theField; 00110 00111 // Initialize the nuclear field map. 00112 theNucleus = nucleus; 00113 theOuterRadius = theNucleus->GetOuterRadius(); 00114 00115 theFieldMap = new std::map <G4int, G4VNuclearField*, std::less<G4int> >; 00116 00117 (*theFieldMap)[G4Proton::Proton()->GetPDGEncoding()] = new G4ProtonField(theNucleus); 00118 (*theFieldMap)[G4Neutron::Neutron()->GetPDGEncoding()] = new G4NeutronField(theNucleus); 00119 (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = new G4AntiProtonField(theNucleus); 00120 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = new G4KaonPlusField(theNucleus); 00121 (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = new G4KaonMinusField(theNucleus); 00122 (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = new G4KaonZeroField(theNucleus); 00123 (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = new G4PionPlusField(theNucleus); 00124 (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = new G4PionMinusField(theNucleus); 00125 (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()] = new G4PionZeroField(theNucleus); 00126 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = new G4SigmaPlusField(theNucleus); 00127 (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = new G4SigmaMinusField(theNucleus); 00128 (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = new G4SigmaZeroField(theNucleus); 00129 00130 theEquationMap = new std::map <G4int, G4Mag_EqRhs*, std::less<G4int> >; 00131 00132 // theField needed by the design of G4Mag_eqRhs 00133 theField = new G4KM_DummyField; //Field not needed for integration 00134 G4KM_OpticalEqRhs * opticalEq; 00135 G4KM_NucleonEqRhs * nucleonEq; 00136 G4double mass; 00137 G4double opticalCoeff; 00138 00139 nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus); 00140 mass = G4Proton::Proton()->GetPDGMass(); 00141 nucleonEq->SetMass(mass); 00142 (*theEquationMap)[G4Proton::Proton()->GetPDGEncoding()] = nucleonEq; 00143 00144 nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus); 00145 mass = G4Neutron::Neutron()->GetPDGMass(); 00146 nucleonEq->SetMass(mass); 00147 (*theEquationMap)[G4Neutron::Neutron()->GetPDGEncoding()] = nucleonEq; 00148 00149 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00150 mass = G4AntiProton::AntiProton()->GetPDGMass(); 00151 opticalCoeff = 00152 (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()]->GetCoeff(); 00153 opticalEq->SetFactor(mass,opticalCoeff); 00154 (*theEquationMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = opticalEq; 00155 00156 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00157 mass = G4KaonPlus::KaonPlus()->GetPDGMass(); 00158 opticalCoeff = 00159 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()]->GetCoeff(); 00160 opticalEq->SetFactor(mass,opticalCoeff); 00161 (*theEquationMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = opticalEq; 00162 00163 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00164 mass = G4KaonMinus::KaonMinus()->GetPDGMass(); 00165 opticalCoeff = 00166 (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()]->GetCoeff(); 00167 opticalEq->SetFactor(mass,opticalCoeff); 00168 (*theEquationMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = opticalEq; 00169 00170 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00171 mass = G4KaonZero::KaonZero()->GetPDGMass(); 00172 opticalCoeff = 00173 (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()]->GetCoeff(); 00174 opticalEq->SetFactor(mass,opticalCoeff); 00175 (*theEquationMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = opticalEq; 00176 00177 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00178 mass = G4PionPlus::PionPlus()->GetPDGMass(); 00179 opticalCoeff = 00180 (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()]->GetCoeff(); 00181 opticalEq->SetFactor(mass,opticalCoeff); 00182 (*theEquationMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = opticalEq; 00183 00184 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00185 mass = G4PionMinus::PionMinus()->GetPDGMass(); 00186 opticalCoeff = 00187 (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()]->GetCoeff(); 00188 opticalEq->SetFactor(mass,opticalCoeff); 00189 (*theEquationMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = opticalEq; 00190 00191 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00192 mass = G4PionZero::PionZero()->GetPDGMass(); 00193 opticalCoeff = 00194 (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()]->GetCoeff(); 00195 opticalEq->SetFactor(mass,opticalCoeff); 00196 (*theEquationMap)[G4PionZero::PionZero()->GetPDGEncoding()] = opticalEq; 00197 00198 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00199 mass = G4SigmaPlus::SigmaPlus()->GetPDGMass(); 00200 opticalCoeff = 00201 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()]->GetCoeff(); 00202 opticalEq->SetFactor(mass,opticalCoeff); 00203 (*theEquationMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = opticalEq; 00204 00205 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00206 mass = G4SigmaMinus::SigmaMinus()->GetPDGMass(); 00207 opticalCoeff = 00208 (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()]->GetCoeff(); 00209 opticalEq->SetFactor(mass,opticalCoeff); 00210 (*theEquationMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = opticalEq; 00211 00212 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 00213 mass = G4SigmaZero::SigmaZero()->GetPDGMass(); 00214 opticalCoeff = 00215 (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()]->GetCoeff(); 00216 opticalEq->SetFactor(mass,opticalCoeff); 00217 (*theEquationMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = opticalEq; 00218 }
void G4RKPropagation::Transport | ( | G4KineticTrackVector & | theActive, | |
const G4KineticTrackVector & | theSpectators, | |||
G4double | theTimeStep | |||
) | [virtual] |
Implements G4VFieldPropagation.
Definition at line 223 of file G4RKPropagation.cc.
References G4KineticTrack::captured, DBL_MAX, encoding, G4cerr, G4cout, G4endl, G4KineticTrack::GetActualMass(), G4KineticTrack::GetDefinition(), G4VNuclearField::GetField(), G4V3DNucleus::GetMass(), G4ParticleDefinition::GetPDGEncoding(), G4KineticTrack::GetPosition(), G4KineticTrack::GetProjectilePotential(), GetSphereIntersectionTimes(), G4KineticTrack::GetState(), G4KineticTrack::GetTrackingMomentum(), G4KineticTrack::gone_out, G4KineticTrack::inside, G4KineticTrack::miss_nucleus, G4Neutron::Neutron(), G4KineticTrack::outside, G4Proton::Proton(), G4KineticTrack::SetState(), G4KineticTrack::SetTrackingMomentum(), and sqr().
00227 { 00228 // reset momentum transfer to field 00229 theMomentumTranfer=G4ThreeVector(0,0,0); 00230 00231 // Loop over tracks 00232 00233 std::vector<G4KineticTrack *>::iterator i; 00234 for(i = active.begin(); i != active.end(); ++i) 00235 { 00236 G4double currTimeStep = timeStep; 00237 G4KineticTrack * kt = *i; 00238 G4int encoding = kt->GetDefinition()->GetPDGEncoding(); 00239 00240 std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter= theFieldMap->find(encoding); 00241 00242 G4VNuclearField* currentField=0; 00243 if ( fieldIter != theFieldMap->end() ) currentField=fieldIter->second; 00244 00245 // debug 00246 // if ( timeStep > 1e30 ) { 00247 // G4cout << " Name :" << kt->GetDefinition()->GetParticleName() << G4endl; 00248 // } 00249 00250 // Get the time of intersections with the nucleus surface. 00251 G4double t_enter, t_leave; 00252 // if the particle does not intersecate with the nucleus go to next particle 00253 if(!GetSphereIntersectionTimes(kt, t_enter, t_leave)) 00254 { 00255 kt->SetState(G4KineticTrack::miss_nucleus); 00256 continue; 00257 } 00258 00259 00260 #ifdef debug_1_RKPropagation 00261 G4cout <<" kt,timeStep, Intersection times tenter, tleave " 00262 <<kt<< " / state= " <<kt->GetState() <<" / " <<" "<< currTimeStep << " / " << t_enter << " / " << t_leave <<G4endl; 00263 #endif 00264 00265 // if the particle is already outside nucleus go to next @@GF should never happen? check! 00266 // does happen for particles added as late.... 00267 // if(t_leave < 0 ) 00268 // { 00269 // throw G4HadronicException(__FILE__, __LINE__, "G4RKPropagation:: Attempt to track particle past a nucleus"); 00270 // continue; 00271 // } 00272 00273 // Apply a straight line propagation for particle types 00274 // not included in the model 00275 if( ! currentField ) 00276 { 00277 if(currTimeStep == DBL_MAX)currTimeStep = t_leave*1.05; 00278 FreeTransport(kt, currTimeStep); 00279 if ( currTimeStep >= t_leave ) 00280 { 00281 if ( kt->GetState() == G4KineticTrack::inside ) 00282 { kt->SetState(G4KineticTrack::gone_out); } 00283 else 00284 { kt->SetState(G4KineticTrack::miss_nucleus);} 00285 } else if (kt->GetState() == G4KineticTrack::outside && currTimeStep >= t_enter ){ 00286 kt->SetState(G4KineticTrack::inside); 00287 } 00288 00289 continue; 00290 } 00291 00292 if(t_enter > 0) // the particle is out. Transport free to the surface 00293 { 00294 if(t_enter > currTimeStep) // the particle won't enter the nucleus 00295 { 00296 FreeTransport(kt, currTimeStep); 00297 continue; 00298 } 00299 else 00300 { 00301 FreeTransport(kt, t_enter); // go to surface 00302 currTimeStep -= t_enter; 00303 t_leave -= t_enter; // time left to leave nucleus 00304 // on the surface the particle loose the barrier energy 00305 // G4double newE = mom.e()-(*theFieldMap)[encoding]->GetBarrier(); 00306 // GetField = Barrier + FermiPotential 00307 G4double newE = kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition()); 00308 00309 if(newE <= kt->GetActualMass()) // the particle cannot enter the nucleus 00310 { 00311 // FixMe: should be "pushed back?" 00312 // for the moment take it past the nucleus, so we'll not worry next time.. 00313 FreeTransport(kt, 1.1*t_leave); // take past nucleus 00314 kt->SetState(G4KineticTrack::miss_nucleus); 00315 // G4cout << "G4RKPropagation: Warning particle cannot enter Nucleus :" << G4endl; 00316 // G4cout << " enter nucleus, E out/in: " << kt->GetTrackingMomentum().e() << " / " << newE <<G4endl; 00317 // G4cout << " the Field "<< currentField->GetField(kt->GetPosition()) << " "<< kt->GetPosition()<<G4endl; 00318 // G4cout << " the particle "<<kt->GetDefinition()->GetParticleName()<<G4endl; 00319 continue; 00320 } 00321 // 00322 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass())); 00323 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE); 00324 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect()); 00325 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass())); 00326 new4Mom*=G4LorentzRotation(boost); 00327 kt->SetTrackingMomentum(new4Mom); 00328 kt->SetState(G4KineticTrack::inside); 00329 00330 /* 00331 G4cout <<" Enter Nucleus - E/Field/Sum: " <<kt->GetTrackingMomentum().e() << " / " 00332 << (*theFieldMap)[encoding]->GetField(kt->GetPosition()) << " / " 00333 << kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition()) 00334 << G4endl 00335 << " Barrier / field just inside nucleus (0.9999*kt->GetPosition())" 00336 << (*theFieldMap)[encoding]->GetBarrier() << " / " 00337 << (*theFieldMap)[encoding]->GetField(0.9999*kt->GetPosition()) 00338 << G4endl; 00339 */ 00340 } 00341 } 00342 00343 // FixMe: should I add a control on theCutOnP here? 00344 // Transport the particle into the nucleus 00345 // G4cerr << "RKPropagation t_leave, curTimeStep " <<t_leave << " " <<currTimeStep<<G4endl; 00346 G4bool is_exiting=false; 00347 if(currTimeStep > t_leave) // particle will exit from the nucleus 00348 { 00349 currTimeStep = t_leave; 00350 is_exiting=true; 00351 } 00352 00353 #ifdef debug_1_RKPropagation 00354 G4cerr << "RKPropagation is_exiting?, t_leave, curTimeStep " <<is_exiting<<" "<<t_leave << " " <<currTimeStep<<G4endl; 00355 G4cout << "RKPropagation Ekin, field, projectile potential, p " 00356 << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " " 00357 << kt->GetPosition()<<" " 00358 << G4endl << currentField->GetField(kt->GetPosition()) << " " 00359 << kt->GetProjectilePotential()<< G4endl 00360 << kt->GetTrackingMomentum() 00361 << G4endl; 00362 #endif 00363 00364 G4LorentzVector momold=kt->GetTrackingMomentum(); 00365 G4ThreeVector posold=kt->GetPosition(); 00366 00367 // if (currentField->GetField(kt->GetPosition()) > kt->GetProjectilePotential() || 00368 if (currTimeStep > 0 && 00369 ! FieldTransport(kt, currTimeStep)) { 00370 FreeTransport(kt,currTimeStep); 00371 } 00372 00373 #ifdef debug_1_RKPropagation 00374 G4cout << "RKPropagation Ekin, field, p " 00375 << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " " 00376 << G4endl << currentField->GetField(kt->GetPosition())<< G4endl 00377 << kt->GetTrackingMomentum() 00378 << G4endl 00379 << "delta p " << momold-kt->GetTrackingMomentum() << G4endl 00380 << "del pos " << posold-kt->GetPosition() 00381 << G4endl; 00382 #endif 00383 00384 // complete the transport 00385 // FixMe: in some cases there could be a significant 00386 // part to do still in the nucleus, or we stepped to far... depending on 00387 // slope of potential 00388 G4double t_in=-1, t_out=0; // set onto boundary. 00389 00390 // should go out, or are already out by a too long step.. 00391 if(is_exiting || 00392 (GetSphereIntersectionTimes(kt, t_in, t_out) &&t_in<0 && t_out<=0 )) // particle is exiting 00393 { 00394 if(t_in < 0 && t_out >= 0) //still inside, transport safely out. 00395 { 00396 // transport free to a position that is surely out of the nucleus, to avoid 00397 // a new transportation and a new adding the barrier next loop. 00398 G4ThreeVector savePos = kt->GetPosition(); 00399 FreeTransport(kt, t_out); 00400 // and evaluate the right the energy 00401 G4double newE=kt->GetTrackingMomentum().e(); 00402 00403 // G4cout << " V pos/savePos << " 00404 // << (*theFieldMap)[encoding]->GetField(kt->GetPosition())<< " / " 00405 // << (*theFieldMap)[encoding]->GetField(savePos) 00406 // << G4endl; 00407 00408 if ( std::abs(currentField->GetField(savePos)) > 0. && 00409 std::abs(currentField->GetField(kt->GetPosition())) > 0.) 00410 { // FixMe GF: savePos/pos may be out of nucleus, where GetField(..)=0 00411 // This wrongly adds or subtracts the Barrier here while 00412 // this is done later. 00413 newE += currentField->GetField(savePos) 00414 - currentField->GetField(kt->GetPosition()); 00415 } 00416 00417 // G4cout << " go border nucleus, E in/border: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl; 00418 00419 if(newE < kt->GetActualMass()) 00420 { 00421 #ifdef debug_1_RKPropagation 00422 G4cout << "RKPropagation-Transport: problem with particle exiting - ignored" << G4endl; 00423 G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl; 00424 #endif 00425 if (kt->GetDefinition() == G4Proton::Proton() || 00426 kt->GetDefinition() == G4Neutron::Neutron() ) { 00427 kt->SetState(G4KineticTrack::captured); 00428 } else { 00429 kt->SetState(G4KineticTrack::gone_out); //@@GF tofix 00430 } 00431 continue; // the particle cannot exit the nucleus 00432 } 00433 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass())); 00434 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE); 00435 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect()); 00436 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass())); 00437 new4Mom*=G4LorentzRotation(boost); 00438 kt->SetTrackingMomentum(new4Mom); 00439 } 00440 // add the potential barrier 00441 // FixMe the Coulomb field is not parallel to mom, this is simple approximation 00442 G4double newE = kt->GetTrackingMomentum().e()+currentField->GetField(kt->GetPosition()); 00443 if(newE < kt->GetActualMass()) 00444 { // the particle cannot exit the nucleus @@@ GF check. 00445 #ifdef debug_1_RKPropagation 00446 G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl; 00447 #endif 00448 if (kt->GetDefinition() == G4Proton::Proton() || 00449 kt->GetDefinition() == G4Neutron::Neutron() ) { 00450 kt->SetState(G4KineticTrack::captured); 00451 } else { 00452 kt->SetState(G4KineticTrack::gone_out); //@@GF tofix 00453 } 00454 continue; 00455 } 00456 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass())); 00457 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE); 00458 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect()); 00459 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass())); 00460 new4Mom*=G4LorentzRotation(boost); 00461 kt->SetTrackingMomentum(new4Mom); 00462 kt->SetState(G4KineticTrack::gone_out); 00463 } 00464 00465 } 00466 00467 }