157 outFile <<
"G4WilsonAbrasionModel is a macroscopic treatment of\n"
158 <<
"nucleus-nucleus collisions using simple geometric arguments.\n"
159 <<
"The smaller projectile nucleus gouges out a part of the larger\n"
160 <<
"target nucleus, leaving a residual nucleus and a fireball\n"
161 <<
"region where the projectile and target intersect. The fireball"
162 <<
"is then treated as a highly excited nuclear fragment. This\n"
163 <<
"model is based on the NUCFRG2 model and is valid for all\n"
164 <<
"projectile energies between 70 MeV/n and 10.1 GeV/n. \n";
252 G4cout <<
"########################################"
253 <<
"########################################"
257 G4cout <<
"Initial projectile A=" <<AP
259 <<
", radius = " <<rP/
fermi <<
" fm"
261 G4cout <<
"Initial target A=" <<AT
263 <<
", radius = " <<rT/
fermi <<
" fm"
265 G4cout <<
"Projectile momentum and Energy/nuc = " <<pP <<
" ," <<E <<
G4endl;
290 G4bool skipInteraction =
false;
291 const G4int maxNumberOfLoops = 1000;
292 G4int loopCounter = -1;
293 while (Dabr == 0 && ++loopCounter < maxNumberOfLoops)
311 skipInteraction =
true;
320 while (r > rPT && ++evtcnt < 1000)
323 r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
329 if (evtcnt >= 1000) {
330 skipInteraction =
true;
340 G4double x = (rPsq + rsq - rTsq) / 2.0 / r;
341 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
342 else CT = 2.0 * std::sqrt(rTsq - rsq);
346 G4double x = (rTsq + rsq - rPsq) / 2.0 / r;
347 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
358 delete theAbrasionGeometry;
360 F = theAbrasionGeometry->
F();
364 for (
G4int i = 0; i<10; ++i)
369 if (n>AP) Dabr = (
G4int) AP;
370 else Dabr = (
G4int) n;
376 if ( loopCounter >= maxNumberOfLoops || skipInteraction ) {
382 G4cout <<
"Particle energy too low to overcome repulsion." <<
G4endl;
383 G4cout <<
"Event rejected and original track maintained" <<
G4endl;
384 G4cout <<
"########################################"
385 <<
"########################################"
388 delete theAbrasionGeometry;
419 for (i=0; i<nSecP; ++i)
422 GetParticle()->GetTotalEnergy();
430 if (DspcP <= 0) DspcP = 0;
431 else if (DspcP > AP-Dabr) DspcP = ((
G4int) AP) - Dabr;
439 G4bool excitationAbsorbedByProjectile =
false;
440 if (fragmentP !=
nullptr)
446 if (excitationAbsorbedByProjectile)
449 if (xP >
B*(AP-Dabr)) xP =
B*(AP-Dabr);
451 lorentzVector.
setE(lorentzVector.e()+xP);
453 TotalEPost += lorentzVector.e();
467 for (i=nSecP; i<nSec; ++i)
470 GetParticle()->GetTotalEnergy();
478 if (DspcT <= 0) DspcT = 0;
479 else if (DspcT > AP-Dabr) DspcT = ((
G4int) AT) - Dabr;
487 if (fragmentT !=
nullptr)
491 if (!excitationAbsorbedByProjectile)
494 if (xT >
B*(AT-Dabr)) xT =
B*(AT-Dabr);
496 lorentzVector.
setE(lorentzVector.e()+xT);
498 TotalEPost += lorentzVector.e();
506 G4double deltaE = TotalEPre - TotalEPost;
510 boost = boost / boost.
mag() * beta;
517 for (i=0; i<nSecP; ++i)
522 lorentzVector.
boost(-boost);
524 pBalance -= lorentzVector.
vect();
536 if (fragmentP !=
nullptr)
546 fragmentP->
SetMomentum(lorentzVector.
boost(-boost * fragmentGroundStateM/fragmentM));
556 G4cout <<
"-----------------------------------" <<
G4endl;
557 G4cout <<
"Secondary nucleons from projectile:" <<
G4endl;
558 G4cout <<
"-----------------------------------" <<
G4endl;
560 for (i=0; i<nSecP; ++i)
572 if (fragmentP !=
nullptr)
581 for (i=nSecP; i<nSec; ++i)
593 if (fragmentT !=
nullptr)
603 if (fragmentP !=
nullptr)
613 G4ReactionProductVector::iterator iter;
614 for (iter = products->begin(); iter != products->end(); ++iter)
618 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
620 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
622 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
627 G4cout <<
" fragmentP = " <<particleName
640 if (fragmentT !=
nullptr)
650 G4ReactionProductVector::iterator iter;
651 for (iter = products->begin(); iter != products->end(); ++iter)
655 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
657 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
659 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
664 G4cout <<
" fragmentT = " <<particleName
673 G4cout <<
"########################################"
674 <<
"########################################"
677 delete theAbrasionGeometry;
717 G4bool isForLoopExitAnticipated =
false;
718 for (
G4int i=0; i<Dabr; ++i)
727 const G4int maxNumberOfLoops = 100000;
728 G4int loopCounter = -1;
729 while (!found && ++loopCounter < maxNumberOfLoops)
734 C2*
G4Exp(-psq/p2sq/2.0) + C3*
G4Exp(-psq/p3sq/2.0) + p/gamma/(0.5*(
G4Exp(p/gamma)-
G4Exp(-p/gamma)));
736 if ( loopCounter >= maxNumberOfLoops )
738 isForLoopExitAnticipated =
true;
763 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
765 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
767 G4double E = std::sqrt(p*p + nucleonMass*nucleonMass)-nucleonMass;
780 if ( ! isForLoopExitAnticipated && Z-Zabr>=1.0 )
810 if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r*rT - rsq - rTsq);
820 if (rT > rP && rsq < rTsq - rPsq) Ct = 2.0 * rP;
821 else if (rP > rT && rsq < rPsq - rTsq) Ct = 2.0 * rT;
823 G4double bP = (rPsq+rsq-rTsq)/2.0/r;
826 G4cerr <<
"########################################"
827 <<
"########################################"
829 G4cerr <<
"ERROR IN G4WilsonAbrasionModel::GetNucleonInducedExcitation"
831 G4cerr <<
"rPsq - bP*bP < 0.0 and cannot be square-rooted" <<
G4endl;
833 G4cerr <<
"########################################"
834 <<
"########################################"
837 Ct = 2.0*std::sqrt(x);
842 Ex += 13.0 * Cl /
fermi /3.0 * (Ct/
fermi - 1.5);
873 G4cout <<
" *****************************************************************"
875 G4cout <<
" Nuclear abrasion model for nuclear-nuclear interactions activated"
877 G4cout <<
" (Written by QinetiQ Ltd for the European Space Agency)"
879 G4cout <<
" *****************************************************************"