287 G4cout <<
"########################################"
288 <<
"########################################"
292 G4cout <<
"Initial projectile A=" <<AP
294 <<
", radius = " <<rP/
fermi <<
" fm"
296 G4cout <<
"Initial target A=" <<AT
298 <<
", radius = " <<rT/
fermi <<
" fm"
300 G4cout <<
"Projectile momentum and Energy/nuc = " <<pP <<
" ," <<E <<
G4endl;
325 G4bool skipInteraction =
false;
326 const G4int maxNumberOfLoops = 1000;
327 G4int loopCounter = -1;
328 while (Dabr == 0 && ++loopCounter < maxNumberOfLoops)
331 if (theAbrasionGeometry)
333 delete theAbrasionGeometry;
334 theAbrasionGeometry = NULL;
351 if (rm >= fradius * rPT) {
352 skipInteraction =
true;
361 while (r > rPT && ++evtcnt < 1000)
364 r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
370 if (evtcnt >= 1000) {
371 skipInteraction =
true;
381 G4double x = (rPsq + rsq - rTsq) / 2.0 / r;
382 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
383 else CT = 2.0 * std::sqrt(rTsq - rsq);
387 G4double x = (rTsq + rsq - rPsq) / 2.0 / r;
388 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
400 F = theAbrasionGeometry->
F();
404 for (
G4int i = 0; i<10; i++)
409 if (n>AP) Dabr = (
G4int) AP;
410 else Dabr = (
G4int) n;
416 if ( loopCounter >= maxNumberOfLoops || skipInteraction ) {
422 G4cout <<
"Particle energy too low to overcome repulsion." <<
G4endl;
423 G4cout <<
"Event rejected and original track maintained" <<
G4endl;
424 G4cout <<
"########################################"
425 <<
"########################################"
455 G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP);
458 for (i=0; i<nSecP; i++)
461 GetParticle()->GetTotalEnergy();
469 if (DspcP <= 0) DspcP = 0;
470 else if (DspcP > AP-Dabr) DspcP = ((
G4int) AP) - Dabr;
478 G4bool excitationAbsorbedByProjectile =
false;
479 if (fragmentP != NULL)
485 if (excitationAbsorbedByProjectile)
486 ExP = GetNucleonInducedExcitation(rP, rT, r);
488 if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr);
490 lorentzVector.
setE(lorentzVector.e()+xP);
492 TotalEPost += lorentzVector.e();
504 G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT);
506 for (i=nSecP; i<nSec; i++)
509 GetParticle()->GetTotalEnergy();
517 if (DspcT <= 0) DspcT = 0;
518 else if (DspcT > AP-Dabr) DspcT = ((
G4int) AT) - Dabr;
526 if (fragmentT != NULL)
530 if (!excitationAbsorbedByProjectile)
531 ExT = GetNucleonInducedExcitation(rT, rP, r);
533 if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr);
535 lorentzVector.
setE(lorentzVector.e()+xT);
537 TotalEPost += lorentzVector.e();
545 G4double deltaE = TotalEPre - TotalEPost;
546 if (deltaE > 0.0 && conserveEnergy)
549 boost = boost / boost.
mag() * beta;
556 for (i=0; i<nSecP; i++)
561 lorentzVector.
boost(-boost);
563 pBalance -= lorentzVector.
vect();
575 if (fragmentP != NULL)
579 if (conserveMomentum)
585 fragmentP->
SetMomentum(lorentzVector.
boost(-boost * fragmentGroundStateM/fragmentM));
595 G4cout <<
"-----------------------------------" <<
G4endl;
596 G4cout <<
"Secondary nucleons from projectile:" <<
G4endl;
597 G4cout <<
"-----------------------------------" <<
G4endl;
599 for (i=0; i<nSecP; i++)
611 if (fragmentP != NULL)
620 for (i=nSecP; i<nSec; i++)
632 if (fragmentT != NULL)
642 if (fragmentP !=NULL)
646 products = theExcitationHandler->
BreakItUp(*fragmentP);
648 products = theExcitationHandlerx->
BreakItUp(*fragmentP);
652 G4ReactionProductVector::iterator iter;
653 for (iter = products->begin(); iter != products->end(); ++iter)
657 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
659 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
661 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
666 G4cout <<
" fragmentP = " <<particleName
679 if (fragmentT != NULL)
683 products = theExcitationHandler->
BreakItUp(*fragmentT);
685 products = theExcitationHandlerx->
BreakItUp(*fragmentT);
689 G4ReactionProductVector::iterator iter;
690 for (iter = products->begin(); iter != products->end(); ++iter)
694 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
696 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
698 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
703 G4cout <<
" fragmentT = " <<particleName
712 G4cout <<
"########################################"
713 <<
"########################################"
716 delete theAbrasionGeometry;
G4double GetExcitationEnergyOfTarget()
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
G4double AtomicMass(const G4double A, const G4double Z) const
G4long G4Poisson(G4double mean)
G4HadSecondary * GetSecondary(size_t i)
G4double GetKineticEnergy() const
void DumpInfo(G4int mode=0) const
G4ParticleDefinition * GetDefinition() const
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
const G4String & GetParticleName() const
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
const G4LorentzVector & GetMomentum() const
HepLorentzVector & boost(double, double, double)
void SetMomentum(const G4LorentzVector &value)
G4double GetKineticEnergy() const
G4double GetEnergyDeposit()
static constexpr double eV
G4double GetGroundStateMass() const
const G4LorentzVector & Get4Momentum() const
G4LorentzVector Get4Momentum() const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
void Set4Momentum(const G4LorentzVector &momentum)
void SetEnergyChange(G4double anEnergy)
Hep3Vector findBoostToCM() const
G4DynamicParticle * GetParticle()
G4double GetExcitationEnergyOfProjectile()
G4HadFinalState theParticleChange
static constexpr double elm_coupling
static constexpr double MeV
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
static constexpr double fermi
G4double GetPDGCharge() const
void SetMomentumChange(const G4ThreeVector &aV)
G4int GetNumberOfSecondaries() const
G4double GetWilsonRadius(G4double A)
G4int GetBaryonNumber() const
G4double GetTotalEnergy() const
CLHEP::HepLorentzVector G4LorentzVector