Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Types | Public Member Functions | List of all members
G4NucleiModel Class Reference

#include <G4NucleiModel.hh>

Public Types

typedef std::pair< std::vector< G4CascadParticle >, std::vector< G4InuclElementaryParticle > > modelLists
 

Public Member Functions

 G4NucleiModel ()
 
 G4NucleiModel (G4int a, G4int z)
 
 G4NucleiModel (G4InuclNuclei *nuclei)
 
 ~G4NucleiModel ()
 
void setVerboseLevel (G4int verbose)
 
void generateModel (G4InuclNuclei *nuclei)
 
void generateModel (G4int a, G4int z)
 
void reset (G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)
 
void printModel () const
 
G4double getDensity (G4int ip, G4int izone) const
 
G4double getFermiMomentum (G4int ip, G4int izone) const
 
G4double getFermiKinetic (G4int ip, G4int izone) const
 
G4double getPotential (G4int ip, G4int izone) const
 
G4double getRadiusUnits () const
 
G4double getRadius () const
 
G4double getRadius (G4int izone) const
 
G4double getVolume (G4int izone) const
 
G4int getNumberOfZones () const
 
G4int getZone (G4double r) const
 
G4int getNumberOfNeutrons () const
 
G4int getNumberOfProtons () const
 
G4bool empty () const
 
G4bool stillInside (const G4CascadParticle &cparticle)
 
G4CascadParticle initializeCascad (G4InuclElementaryParticle *particle)
 
void initializeCascad (G4InuclNuclei *bullet, G4InuclNuclei *target, modelLists &output)
 
std::pair< G4int, G4int > getTypesOfNucleonsInvolved () const
 
void generateParticleFate (G4CascadParticle &cparticle, G4ElementaryParticleCollider *theEPCollider, std::vector< G4CascadParticle > &cascade)
 
G4bool isProjectile (const G4CascadParticle &cparticle) const
 
G4bool worthToPropagate (const G4CascadParticle &cparticle) const
 
G4InuclElementaryParticle generateNucleon (G4int type, G4int zone) const
 
G4LorentzVector generateNucleonMomentum (G4int type, G4int zone) const
 
G4double absorptionCrossSection (G4double e, G4int type) const
 
G4double totalCrossSection (G4double ke, G4int rtype) const
 

Detailed Description

Definition at line 82 of file G4NucleiModel.hh.

Member Typedef Documentation

◆ modelLists

typedef std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > G4NucleiModel::modelLists

Definition at line 152 of file G4NucleiModel.hh.

Constructor & Destructor Documentation

◆ G4NucleiModel() [1/3]

G4NucleiModel::G4NucleiModel ( )

Definition at line 194 of file G4NucleiModel.cc.

195 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
196 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
197 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
198 current_nucl2(0) {}

◆ G4NucleiModel() [2/3]

G4NucleiModel::G4NucleiModel ( G4int  a,
G4int  z 
)

Definition at line 200 of file G4NucleiModel.cc.

201 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
202 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
203 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
204 current_nucl2(0) {
205 generateModel(a,z);
206}
G4NucleiModel::generateModel
void generateModel(G4InuclNuclei *nuclei)
Definition: G4NucleiModel.cc:237

◆ G4NucleiModel() [3/3]

G4NucleiModel::G4NucleiModel ( G4InuclNuclei *  nuclei)
explicit

Definition at line 208 of file G4NucleiModel.cc.

209 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
210 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
211 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
212 current_nucl2(0) {
213 generateModel(nuclei);
214}

◆ ~G4NucleiModel()

G4NucleiModel::~G4NucleiModel ( )

Definition at line 216 of file G4NucleiModel.cc.

216 {
217 delete theNucleus;
218 theNucleus = 0;
219}

Member Function Documentation

◆ absorptionCrossSection()

G4double G4NucleiModel::absorptionCrossSection ( G4double  e,
G4int  type 
) const

Definition at line 1646 of file G4NucleiModel.cc.

1646 {
1647 if (type != pionPlus && type != pionMinus && type != pionZero &&
1648 type != photon) {
1649 G4cerr << " absorptionCrossSection only valid for incident pions" << G4endl;
1650 return 0.;
1651 }
1652
1653 G4double csec = 0.;
1654
1655 // Pion absorption is parametrized for low vs. medium energy
1656 if (type == pionPlus || type == pionMinus || type == pionZero) {
1657 if (ke < 0.3) csec = (0.1106 / std::sqrt(ke) - 0.8
1658 + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056) );
1659 else if (ke < 1.0) csec = 3.6735 * (1.0-ke)*(1.0-ke);
1660 }
1661
1662 // Photon cross-section is binned from parametrization by Mi. Kossov
1663 if (type == photon) {
1664 static const G4double gammaQDscale = G4CascadeParameters::gammaQDScale();
1665 static const G4double kebins[] =
1666 { 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
1667 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
1668 2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
1669 static const G4double gammaD[] =
1670 { 2.8, 1.3, 0.89, 0.56, 0.38, 0.27, 0.19, 0.14, 0.098,
1671 0.071, 0.054, 0.0003, 0.0007, 0.0027, 0.0014, 0.001, 0.0012, 0.0005,
1672 0.0003, 0.0002,0.0002, 0.0002, 0.0002, 0.0002, 0.0001, 0.0001, 0.0001,
1673 0.0001, 0.0001, 0.0001 };
1674 static G4CascadeInterpolator<30> interp(kebins);
1675 csec = interp.interpolate(ke, gammaD) * gammaQDscale;
1676 }
1677
1678 if (csec < 0.0) csec = 0.0;
1679
1680 if (verboseLevel > 2) {
1681 G4cout << " ekin " << ke << " abs. csec " << csec << " mb" << G4endl;
1682 }
1683
1684 return crossSectionUnits * csec;
1685}
photon
@ photon
Definition: G4DataQuestionaire.hh:35
G4double
double G4double
Definition: G4Types.hh:64
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4cout
G4DLLIMPORT std::ostream G4cout
G4CascadeParameters::gammaQDScale
static G4double gammaQDScale()
Definition: G4CascadeParameters.hh:59

◆ empty()

G4bool G4NucleiModel::empty ( ) const
inline

Definition at line 141 of file G4NucleiModel.hh.

141 {
142 return neutronNumberCurrent < 1 && protonNumberCurrent < 1;
143 }

Referenced by G4IntraNucleiCascader::generateCascade().

◆ generateModel() [1/2]

void G4NucleiModel::generateModel ( G4int  a,
G4int  z 
)

Definition at line 241 of file G4NucleiModel.cc.

241 {
242 if (verboseLevel) {
243 G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
244 << G4endl;
245 }
246
247 // If model already built, just return; otherwise intialize everything
248 if (a == A && z == Z) {
249 if (verboseLevel > 1) G4cout << " model already generated" << z << G4endl;
250 reset();
251 return;
252 }
253
254 A = a;
255 Z = z;
256 delete theNucleus;
257 theNucleus = new G4InuclNuclei(A,Z); // For conservation checking
258
259 neutronNumber = A - Z;
260 protonNumber = Z;
261 reset();
262
263 if (verboseLevel > 3) {
264 G4cout << " crossSectionUnits = " << crossSectionUnits << G4endl
265 << " radiusUnits = " << radiusUnits << G4endl
266 << " skinDepth = " << skinDepth << G4endl
267 << " radiusScale = " << radiusScale << G4endl
268 << " radiusScale2 = " << radiusScale2 << G4endl
269 << " radiusForSmall = " << radiusForSmall << G4endl
270 << " radScaleAlpha = " << radScaleAlpha << G4endl
271 << " fermiMomentum = " << fermiMomentum << G4endl
272 << " piTimes4thirds = " << piTimes4thirds << G4endl;
273 }
274
275 G4double nuclearRadius; // Nuclear radius computed from A
276 if (A>4) nuclearRadius = radiusScale*G4cbrt(A) + radiusScale2/G4cbrt(A);
277 else nuclearRadius = radiusForSmall * (A==4 ? radScaleAlpha : 1.);
278
279 // This will be used to pre-allocate lots of arrays below
280 number_of_zones = (A < 5) ? 1 : (A < 100) ? 3 : 6;
281
282 // Clear all parameters arrays for reloading
283 binding_energies.clear();
284 nucleon_densities.clear();
285 zone_potentials.clear();
286 fermi_momenta.clear();
287 zone_radii.clear();
288 zone_volumes.clear();
289
290 fillBindingEnergies();
291 fillZoneRadii(nuclearRadius);
292
293 G4double tot_vol = fillZoneVolumes(nuclearRadius); // Woods-Saxon integral
294
295 fillPotentials(proton, tot_vol); // Protons
296 fillPotentials(neutron, tot_vol); // Neutrons
297
298 // Additional flat zone potentials for other hadrons
299 const std::vector<G4double> vp(number_of_zones, (A>4)?pion_vp:pion_vp_small);
300 const std::vector<G4double> kp(number_of_zones, kaon_vp);
301 const std::vector<G4double> hp(number_of_zones, hyperon_vp);
302
303 zone_potentials.push_back(vp);
304 zone_potentials.push_back(kp);
305 zone_potentials.push_back(hp);
306
307 nuclei_radius = zone_radii.back();
308 nuclei_volume = std::accumulate(zone_volumes.begin(),zone_volumes.end(),0.);
309
310 if (verboseLevel > 3) printModel();
311}
neutron
@ neutron
Definition: G4DataQuestionaire.hh:35
G4NucleiModel::printModel
void printModel() const
Definition: G4NucleiModel.cc:540
G4NucleiModel::reset
void reset(G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)
Definition: G4NucleiModel.cc:224

◆ generateModel() [2/2]

void G4NucleiModel::generateModel ( G4InuclNuclei *  nuclei)

Definition at line 237 of file G4NucleiModel.cc.

237 {
238 generateModel(nuclei->getA(), nuclei->getZ());
239}

Referenced by G4NucleiModel(), generateModel(), and G4IntraNucleiCascader::initialize().

◆ generateNucleon()

G4InuclElementaryParticle G4NucleiModel::generateNucleon ( G4int  type,
G4int  zone 
) const

Definition at line 585 of file G4NucleiModel.cc.

585 {
586 if (verboseLevel > 1) {
587 G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
588 }
589
590 G4LorentzVector mom = generateNucleonMomentum(type, zone);
591 return G4InuclElementaryParticle(mom, type);
592}
G4NucleiModel::generateNucleonMomentum
G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const
Definition: G4NucleiModel.cc:576

◆ generateNucleonMomentum()

G4LorentzVector G4NucleiModel::generateNucleonMomentum ( G4int  type,
G4int  zone 
) const

Definition at line 576 of file G4NucleiModel.cc.

576 {
577 G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
578 G4double mass = G4InuclElementaryParticle::getParticleMass(type);
579
580 return generateWithRandomAngles(pmod, mass);
581}
G4InuclElementaryParticle::getParticleMass
static G4double getParticleMass(G4int type)
Definition: G4InuclElementaryParticle.cc:217
G4NucleiModel::getFermiMomentum
G4double getFermiMomentum(G4int ip, G4int izone) const
Definition: G4NucleiModel.hh:105
G4InuclSpecialFunctions::generateWithRandomAngles
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
Definition: G4InuclSpecialFunctions.cc:124

Referenced by generateNucleon().

◆ generateParticleFate()

void G4NucleiModel::generateParticleFate ( G4CascadParticle &  cparticle,
G4ElementaryParticleCollider *  theEPCollider,
std::vector< G4CascadParticle > &  cascade 
)

Definition at line 789 of file G4NucleiModel.cc.

792 {
793 if (verboseLevel > 1)
794 G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
795
796 if (verboseLevel > 2) G4cout << " cparticle: " << cparticle << G4endl;
797
798 // Create four-vector checking
799#ifdef G4CASCADE_CHECK_ECONS
800 G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel"); // Second arg is in GeV
801 balance.setVerboseLevel(verboseLevel);
802#endif
803
804 outgoing_cparticles.clear(); // Clear return buffer for this event
805 generateInteractionPartners(cparticle); // Fills "thePartners" data
806
807 if (thePartners.empty()) { // smth. is wrong -> needs special treatment
808 if (verboseLevel)
809 G4cerr << " generateParticleFate-> got empty interaction-partners list "
810 << G4endl;
811 return;
812 }
813
814 G4int npart = thePartners.size(); // Last item is a total-path placeholder
815
816 if (npart == 1) { // cparticle is on the next zone entry
817 cparticle.propagateAlongThePath(thePartners[0].second);
818 cparticle.incrementCurrentPath(thePartners[0].second);
819 boundaryTransition(cparticle);
820 outgoing_cparticles.push_back(cparticle);
821
822 if (verboseLevel > 2) G4cout << " next zone \n" << cparticle << G4endl;
823 } else { // there are possible interactions
824 if (verboseLevel > 1)
825 G4cout << " processing " << npart-1 << " possible interactions" << G4endl;
826
827 G4ThreeVector old_position = cparticle.getPosition();
828 G4InuclElementaryParticle& bullet = cparticle.getParticle();
829 G4bool no_interaction = true;
830 G4int zone = cparticle.getCurrentZone();
831
832 for (G4int i=0; i<npart-1; i++) { // Last item is a total-path placeholder
833 if (i > 0) cparticle.updatePosition(old_position);
834
835 G4InuclElementaryParticle& target = thePartners[i].first;
836
837 if (verboseLevel > 3) {
838 if (target.quasi_deutron()) G4cout << " try absorption: ";
839 G4cout << " target " << target.type() << " bullet " << bullet.type()
840 << G4endl;
841 }
842
843 EPCoutput.reset();
844 theEPCollider->collide(&bullet, &target, EPCoutput);
845
846 // If collision failed, exit loop over partners
847 if (EPCoutput.numberOfOutgoingParticles() == 0) break;
848
849 if (verboseLevel > 2) {
850 EPCoutput.printCollisionOutput();
851#ifdef G4CASCADE_CHECK_ECONS
852 balance.collide(&bullet, &target, EPCoutput);
853 balance.okay(); // Do checks, but ignore result
854#endif
855 }
856
857 // Get list of outgoing particles for evaluation
858 std::vector<G4InuclElementaryParticle>& outgoing_particles =
859 EPCoutput.getOutgoingParticles();
860
861 if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails
862
863 // Trailing effect: reject interaction at previously hit nucleon
864 cparticle.propagateAlongThePath(thePartners[i].second);
865 G4ThreeVector new_position = cparticle.getPosition();
866
867 if (!passTrailing(new_position)) continue;
868 collisionPts.push_back(new_position);
869
870 // Sort particles according to beta (fastest first)
871 std::sort(outgoing_particles.begin(), outgoing_particles.end(),
872 G4ParticleLargerBeta() );
873
874 if (verboseLevel > 2)
875 G4cout << " adding " << outgoing_particles.size()
876 << " output particles" << G4endl;
877
878 // NOTE: Embedded temporary is optimized away (no copying gets done)
879 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) {
880 outgoing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip],
881 new_position, zone,
882 0.0, 0));
883 }
884
885 no_interaction = false;
886 current_nucl1 = 0;
887 current_nucl2 = 0;
888
889#ifdef G4CASCADE_DEBUG_CHARGE
890 {
891 G4double out_charge = 0.0;
892
893 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++)
894 out_charge += outgoing_particles[ip].getCharge();
895
896 G4cout << " multiplicity " << outgoing_particles.size()
897 << " bul type " << bullet.type()
898 << " targ type " << target.type()
899 << "\n initial charge "
900 << bullet.getCharge() + target.getCharge()
901 << " out charge " << out_charge << G4endl;
902 }
903#endif
904
905 if (verboseLevel > 2)
906 G4cout << " partner type " << target.type() << G4endl;
907
908 if (target.nucleon()) {
909 current_nucl1 = target.type();
910 } else {
911 if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
912 current_nucl1 = (target.type() - 100) / 10;
913 current_nucl2 = target.type() - 100 - 10 * current_nucl1;
914 }
915
916 if (current_nucl1 == 1) {
917 if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
918 protonNumberCurrent--;
919 } else {
920 if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
921 neutronNumberCurrent--;
922 }
923
924 if (current_nucl2 == 1) {
925 if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
926 protonNumberCurrent--;
927 } else if (current_nucl2 == 2) {
928 if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
929 neutronNumberCurrent--;
930 }
931
932 break;
933 } // loop over partners
934
935 if (no_interaction) { // still no interactions
936 if (verboseLevel > 1) G4cout << " no interaction " << G4endl;
937
938 // For conservation checking (below), get particle before updating
939 static G4InuclElementaryParticle prescatCP; // Avoid memory churn
940 prescatCP = cparticle.getParticle();
941
942 // Last "partner" is just a total-path placeholder
943 cparticle.updatePosition(old_position);
944 cparticle.propagateAlongThePath(thePartners[npart-1].second);
945 cparticle.incrementCurrentPath(thePartners[npart-1].second);
946 boundaryTransition(cparticle);
947 outgoing_cparticles.push_back(cparticle);
948
949 // Check conservation for simple scattering (ignore target nucleus!)
950#ifdef G4CASCADE_CHECK_ECONS
951 if (verboseLevel > 2) {
952 balance.collide(&prescatCP, 0, outgoing_cparticles);
953 balance.okay(); // Report violations, but don't act on them
954 }
955#endif
956 } // if (no_interaction)
957 } // if (npart == 1) [else]
958
959 return;
960}
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4CascadParticle::incrementCurrentPath
void incrementCurrentPath(G4double npath)
Definition: G4CascadParticle.hh:92
G4CascadParticle::getParticle
const G4InuclElementaryParticle & getParticle() const
Definition: G4CascadParticle.hh:65
G4CascadParticle::updatePosition
void updatePosition(const G4ThreeVector &pos)
Definition: G4CascadParticle.hh:80
G4CascadParticle::propagateAlongThePath
void propagateAlongThePath(G4double path)
Definition: G4CascadParticle.cc:129
G4CascadParticle::getCurrentZone
G4int getCurrentZone() const
Definition: G4CascadParticle.hh:96
G4CascadParticle::getPosition
const G4ThreeVector & getPosition() const
Definition: G4CascadParticle.hh:79
G4CollisionOutput::numberOfOutgoingParticles
G4int numberOfOutgoingParticles() const
Definition: G4CollisionOutput.hh:123
G4CollisionOutput::printCollisionOutput
void printCollisionOutput(std::ostream &os=G4cout) const
Definition: G4CollisionOutput.cc:276
G4CollisionOutput::getOutgoingParticles
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
Definition: G4CollisionOutput.hh:125
G4ElementaryParticleCollider::collide
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
Definition: G4ElementaryParticleCollider.cc:115
G4InuclParticle::getCharge
G4double getCharge() const
Definition: G4InuclParticle.hh:110

Referenced by G4IntraNucleiCascader::generateCascade().

◆ getDensity()

G4double G4NucleiModel::getDensity ( G4int  ip,
G4int  izone 
) const
inline

Definition at line 101 of file G4NucleiModel.hh.

101 {
102 return nucleon_densities[ip - 1][izone];
103 }

Referenced by printModel().

◆ getFermiKinetic()

G4double G4NucleiModel::getFermiKinetic ( G4int  ip,
G4int  izone 
) const

Definition at line 562 of file G4NucleiModel.cc.

562 {
563 G4double ekin = 0.0;
564
565 if (ip < 3 && izone < number_of_zones) { // ip for proton/neutron only
566 G4double pfermi = fermi_momenta[ip - 1][izone];
567 G4double mass = G4InuclElementaryParticle::getParticleMass(ip);
568
569 ekin = std::sqrt(pfermi*pfermi + mass*mass) - mass;
570 }
571 return ekin;
572}

Referenced by worthToPropagate().

◆ getFermiMomentum()

G4double G4NucleiModel::getFermiMomentum ( G4int  ip,
G4int  izone 
) const
inline

Definition at line 105 of file G4NucleiModel.hh.

105 {
106 return fermi_momenta[ip - 1][izone];
107 }

Referenced by generateNucleonMomentum(), and printModel().

◆ getNumberOfNeutrons()

G4int G4NucleiModel::getNumberOfNeutrons ( ) const
inline

Definition at line 138 of file G4NucleiModel.hh.

138{ return neutronNumberCurrent; }

Referenced by G4IntraNucleiCascader::generateCascade().

◆ getNumberOfProtons()

G4int G4NucleiModel::getNumberOfProtons ( ) const
inline

Definition at line 139 of file G4NucleiModel.hh.

139{ return protonNumberCurrent; }

Referenced by G4IntraNucleiCascader::generateCascade().

◆ getNumberOfZones()

G4int G4NucleiModel::getNumberOfZones ( ) const
inline

Definition at line 132 of file G4NucleiModel.hh.

132{ return number_of_zones; }

◆ getPotential()

G4double G4NucleiModel::getPotential ( G4int  ip,
G4int  izone 
) const
inline

Definition at line 111 of file G4NucleiModel.hh.

111 {
112 if (ip == 9) return 0.0; // Special case for photons
113 G4int ip0 = ip < 3 ? ip - 1 : 2;
114 if (ip > 10 && ip < 18) ip0 = 3;
115 if (ip > 20) ip0 = 4;
116 return izone < number_of_zones ? zone_potentials[ip0][izone] : 0.0;
117 }

Referenced by printModel().

◆ getRadius() [1/2]

G4double G4NucleiModel::getRadius ( ) const
inline

Definition at line 122 of file G4NucleiModel.hh.

122{ return nuclei_radius; }

◆ getRadius() [2/2]

G4double G4NucleiModel::getRadius ( G4int  izone) const
inline

Definition at line 123 of file G4NucleiModel.hh.

123 {
124 return ( (izone<0) ? 0.
125 : (izone<number_of_zones) ? zone_radii[izone] : nuclei_radius);
126 }

◆ getRadiusUnits()

G4double G4NucleiModel::getRadiusUnits ( ) const
inline

Definition at line 120 of file G4NucleiModel.hh.

120{ return radiusUnits*CLHEP::fermi; }

Referenced by G4IntraNucleiCascader::processSecondary().

◆ getTypesOfNucleonsInvolved()

std::pair< G4int, G4int > G4NucleiModel::getTypesOfNucleonsInvolved ( ) const
inline

Definition at line 158 of file G4NucleiModel.hh.

158 {
159 return std::pair<G4int, G4int>(current_nucl1, current_nucl2);
160 }

Referenced by G4IntraNucleiCascader::generateCascade().

◆ getVolume()

G4double G4NucleiModel::getVolume ( G4int  izone) const
inline

Definition at line 127 of file G4NucleiModel.hh.

127 {
128 return ( (izone<0) ? 0.
129 : (izone<number_of_zones) ? zone_volumes[izone] : nuclei_volume);
130 }

◆ getZone()

G4int G4NucleiModel::getZone ( G4double  r) const
inline

Definition at line 133 of file G4NucleiModel.hh.

133 {
134 for (G4int iz=0; iz<number_of_zones; iz++) if (r<zone_radii[iz]) return iz;
135 return number_of_zones;
136 }

Referenced by G4IntraNucleiCascader::processSecondary().

◆ initializeCascad() [1/2]

G4CascadParticle G4NucleiModel::initializeCascad ( G4InuclElementaryParticle *  particle)

Definition at line 1164 of file G4NucleiModel.cc.

1164 {
1165 if (verboseLevel > 1) {
1166 G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl;
1167 }
1168
1169 const G4double large = 1000.0;
1170
1171 // FIXME: Previous version generated random sin(theta), then used -cos(theta)
1172 // Using generateWithRandomAngles changes result!
1173 // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect();
1174 G4double costh = std::sqrt(1.0 - inuclRndm());
1175 G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius);
1176
1177 G4CascadParticle cpart(*particle, pos, number_of_zones, large, 0);
1178
1179 if (verboseLevel > 2) G4cout << cpart << G4endl;
1180
1181 return cpart;
1182}
G4InuclSpecialFunctions::generateWithFixedTheta
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
Definition: G4InuclSpecialFunctions.cc:109

Referenced by G4IntraNucleiCascader::setupCascade().

◆ initializeCascad() [2/2]

void G4NucleiModel::initializeCascad ( G4InuclNuclei *  bullet,
G4InuclNuclei *  target,
modelLists &  output 
)

Definition at line 1184 of file G4NucleiModel.cc.

1186 {
1187 if (verboseLevel) {
1188 G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)"
1189 << G4endl;
1190 }
1191
1192 const G4double large = 1000.0;
1193 const G4double max_a_for_cascad = 5.0;
1194 const G4double ekin_cut = 2.0;
1195 const G4double small_ekin = 1.0e-6;
1196 const G4double r_large2for3 = 62.0;
1197 const G4double r0forAeq3 = 3.92;
1198 const G4double s3max = 6.5;
1199 const G4double r_large2for4 = 69.14;
1200 const G4double r0forAeq4 = 4.16;
1201 const G4double s4max = 7.0;
1202 const G4int itry_max = 100;
1203
1204 // Convenient references to input buffer contents
1205 std::vector<G4CascadParticle>& casparticles = output.first;
1206 std::vector<G4InuclElementaryParticle>& particles = output.second;
1207
1208 casparticles.clear(); // Reset input buffer to fill this event
1209 particles.clear();
1210
1211 // first decide whether it will be cascad or compound final nuclei
1212 G4int ab = bullet->getA();
1213 G4int zb = bullet->getZ();
1214 G4int at = target->getA();
1215 G4int zt = target->getZ();
1216
1217 G4double massb = bullet->getMass(); // For creating LorentzVectors below
1218
1219 if (ab < max_a_for_cascad) {
1220
1221 G4double benb = bindingEnergy(ab,zb)/GeV / G4double(ab);
1222 G4double bent = bindingEnergy(at,zt)/GeV / G4double(at);
1223 G4double ben = benb < bent ? bent : benb;
1224
1225 if (bullet->getKineticEnergy()/ab > ekin_cut*ben) {
1226 G4int itryg = 0;
1227
1228 while (casparticles.size() == 0 && itryg < itry_max) {
1229 itryg++;
1230 particles.clear();
1231
1232 // nucleons coordinates and momenta in nuclei rest frame
1233 coordinates.clear();
1234 momentums.clear();
1235
1236 if (ab < 3) { // deuteron, simplest case
1237 G4double r = 2.214 - 3.4208 * std::log(1.0 - 0.981 * inuclRndm());
1238 G4ThreeVector coord1 = generateWithRandomAngles(r).vect();
1239 coordinates.push_back(coord1);
1240 coordinates.push_back(-coord1);
1241
1242 G4double p = 0.0;
1243 G4bool bad = true;
1244 G4int itry = 0;
1245
1246 while (bad && itry < itry_max) {
1247 itry++;
1248 p = 456.0 * inuclRndm();
1249
1250 if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
1251 p * r > 312.0) bad = false;
1252 }
1253
1254 if (itry == itry_max)
1255 if (verboseLevel > 2){
1256 G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl;
1257 }
1258
1259 p = 0.0005 * p;
1260
1261 if (verboseLevel > 2){
1262 G4cout << " p nuc " << p << G4endl;
1263 }
1264
1265 G4LorentzVector mom = generateWithRandomAngles(p, massb);
1266
1267 momentums.push_back(mom);
1268 mom.setVect(-mom.vect());
1269 momentums.push_back(-mom);
1270 } else {
1271 G4ThreeVector coord1;
1272
1273 G4bool badco = true;
1274
1275 G4int itry = 0;
1276
1277 if (ab == 3) {
1278 while (badco && itry < itry_max) {
1279 if (itry > 0) coordinates.clear();
1280 itry++;
1281 G4int i(0);
1282
1283 for (i = 0; i < 2; i++) {
1284 G4int itry1 = 0;
1285 G4double ss, u, rho;
1286 G4double fmax = std::exp(-0.5) / std::sqrt(0.5);
1287
1288 while (itry1 < itry_max) {
1289 itry1++;
1290 ss = -std::log(inuclRndm());
1291 u = fmax * inuclRndm();
1292 rho = std::sqrt(ss) * std::exp(-ss);
1293
1294 if (rho > u && ss < s3max) {
1295 ss = r0forAeq3 * std::sqrt(ss);
1296 coord1 = generateWithRandomAngles(ss).vect();
1297 coordinates.push_back(coord1);
1298
1299 if (verboseLevel > 2){
1300 G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1301 }
1302 break;
1303 }
1304 }
1305
1306 if (itry1 == itry_max) { // bad case
1307 coord1.set(10000.,10000.,10000.);
1308 coordinates.push_back(coord1);
1309 break;
1310 }
1311 }
1312
1313 coord1 = -coordinates[0] - coordinates[1];
1314 if (verboseLevel > 2) {
1315 G4cout << " 3 r " << coord1.mag() << G4endl;
1316 }
1317
1318 coordinates.push_back(coord1);
1319
1320 G4bool large_dist = false;
1321
1322 for (i = 0; i < 2; i++) {
1323 for (G4int j = i+1; j < 3; j++) {
1324 G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1325
1326 if (verboseLevel > 2) {
1327 G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1328 }
1329
1330 if (r2 > r_large2for3) {
1331 large_dist = true;
1332
1333 break;
1334 }
1335 }
1336
1337 if (large_dist) break;
1338 }
1339
1340 if(!large_dist) badco = false;
1341
1342 }
1343
1344 } else { // a >= 4
1345 G4double b = 3./(ab - 2.0);
1346 G4double b1 = 1.0 - b / 2.0;
1347 G4double u = b1 + std::sqrt(b1 * b1 + b);
1348 G4double fmax = (1.0 + u/b) * u * std::exp(-u);
1349
1350 while (badco && itry < itry_max) {
1351
1352 if (itry > 0) coordinates.clear();
1353 itry++;
1354 G4int i(0);
1355
1356 for (i = 0; i < ab-1; i++) {
1357 G4int itry1 = 0;
1358 G4double ss;
1359
1360 while (itry1 < itry_max) {
1361 itry1++;
1362 ss = -std::log(inuclRndm());
1363 u = fmax * inuclRndm();
1364
1365 if (std::sqrt(ss) * std::exp(-ss) * (1.0 + ss/b) > u
1366 && ss < s4max) {
1367 ss = r0forAeq4 * std::sqrt(ss);
1368 coord1 = generateWithRandomAngles(ss).vect();
1369 coordinates.push_back(coord1);
1370
1371 if (verboseLevel > 2) {
1372 G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1373 }
1374
1375 break;
1376 }
1377 }
1378
1379 if (itry1 == itry_max) { // bad case
1380 coord1.set(10000.,10000.,10000.);
1381 coordinates.push_back(coord1);
1382 break;
1383 }
1384 }
1385
1386 coord1 *= 0.0; // Cheap way to reset
1387 for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j];
1388
1389 coordinates.push_back(coord1);
1390
1391 if (verboseLevel > 2){
1392 G4cout << " last r " << coord1.mag() << G4endl;
1393 }
1394
1395 G4bool large_dist = false;
1396
1397 for (i = 0; i < ab-1; i++) {
1398 for (G4int j = i+1; j < ab; j++) {
1399
1400 G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1401
1402 if (verboseLevel > 2){
1403 G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1404 }
1405
1406 if (r2 > r_large2for4) {
1407 large_dist = true;
1408
1409 break;
1410 }
1411 }
1412
1413 if (large_dist) break;
1414 }
1415
1416 if (!large_dist) badco = false;
1417 }
1418 }
1419
1420 if(badco) {
1421 G4cout << " can not generate the nucleons coordinates for a "
1422 << ab << G4endl;
1423
1424 casparticles.clear(); // Return empty buffer on error
1425 particles.clear();
1426 return;
1427
1428 } else { // momentums
1429 G4double p, u, x;
1430 G4LorentzVector mom;
1431 //G4bool badp = True;
1432
1433 for (G4int i = 0; i < ab - 1; i++) {
1434 G4int itry2 = 0;
1435
1436 while(itry2 < itry_max) {
1437 itry2++;
1438 u = -std::log(0.879853 - 0.8798502 * inuclRndm());
1439 x = u * std::exp(-u);
1440
1441 if(x > inuclRndm()) {
1442 p = std::sqrt(0.01953 * u);
1443 mom = generateWithRandomAngles(p, massb);
1444 momentums.push_back(mom);
1445
1446 break;
1447 }
1448 } // while (itry2 < itry_max)
1449
1450 if(itry2 == itry_max) {
1451 G4cout << " can not generate proper momentum for a "
1452 << ab << G4endl;
1453
1454 casparticles.clear(); // Return empty buffer on error
1455 particles.clear();
1456 return;
1457 }
1458 } // for (i=0 ...
1459 // last momentum
1460
1461 mom *= 0.; // Cheap way to reset
1462 mom.setE(bullet->getEnergy()+target->getEnergy());
1463
1464 for(G4int j=0; j< ab-1; j++) mom -= momentums[j];
1465
1466 momentums.push_back(mom);
1467 }
1468 }
1469
1470 // Coordinates and momenta at rest are generated, now back to the lab
1471 G4double rb = 0.0;
1472 G4int i(0);
1473
1474 for(i = 0; i < G4int(coordinates.size()); i++) {
1475 G4double rp = coordinates[i].mag2();
1476
1477 if(rp > rb) rb = rp;
1478 }
1479
1480 // nuclei i.p. as a whole
1481 G4double s1 = std::sqrt(inuclRndm());
1482 G4double phi = randomPHI();
1483 G4double rz = (nuclei_radius + rb) * s1;
1484 G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi),
1485 -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1));
1486
1487 for (i = 0; i < G4int(coordinates.size()); i++) {
1488 coordinates[i] += global_pos;
1489 }
1490
1491 // all nucleons at rest
1492 raw_particles.clear();
1493
1494 for (G4int ipa = 0; ipa < ab; ipa++) {
1495 G4int knd = ipa < zb ? 1 : 2;
1496 raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
1497 }
1498
1499 G4InuclElementaryParticle dummy(small_ekin, 1);
1500 G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet);
1501 toTheBulletRestFrame.toTheTargetRestFrame();
1502
1503 particleIterator ipart;
1504
1505 for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
1506 ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum()));
1507 }
1508
1509 // fill cascad particles and outgoing particles
1510
1511 for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
1512 G4LorentzVector mom = raw_particles[ip].getMomentum();
1513 G4double pmod = mom.rho();
1514 G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod;
1515 G4double det = t0 * t0 + nuclei_radius * nuclei_radius
1516 - coordinates[ip].mag2();
1517 G4double tr = -1.0;
1518
1519 if(det > 0.0) {
1520 G4double t1 = t0 + std::sqrt(det);
1521 G4double t2 = t0 - std::sqrt(det);
1522
1523 if(std::fabs(t1) <= std::fabs(t2)) {
1524 if(t1 > 0.0) {
1525 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1526 }
1527
1528 if(tr < 0.0 && t2 > 0.0) {
1529
1530 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1531 }
1532
1533 } else {
1534 if(t2 > 0.0) {
1535
1536 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1537 }
1538
1539 if(tr < 0.0 && t1 > 0.0) {
1540 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1541 }
1542 }
1543
1544 }
1545
1546 if(tr >= 0.0) { // cascad particle
1547 coordinates[ip] += mom.vect()*tr / pmod;
1548 casparticles.push_back(G4CascadParticle(raw_particles[ip],
1549 coordinates[ip],
1550 number_of_zones, large, 0));
1551
1552 } else {
1553 particles.push_back(raw_particles[ip]);
1554 }
1555 }
1556 }
1557
1558 if(casparticles.size() == 0) {
1559 particles.clear();
1560
1561 G4cout << " can not generate proper distribution for " << itry_max
1562 << " steps " << G4endl;
1563 }
1564 }
1565 }
1566
1567 if(verboseLevel > 2){
1568 G4int ip(0);
1569
1570 G4cout << " cascad particles: " << casparticles.size() << G4endl;
1571 for(ip = 0; ip < G4int(casparticles.size()); ip++)
1572 G4cout << casparticles[ip] << G4endl;
1573
1574 G4cout << " outgoing particles: " << particles.size() << G4endl;
1575 for(ip = 0; ip < G4int(particles.size()); ip++)
1576 G4cout << particles[ip] << G4endl;
1577 }
1578
1579 return; // Buffer has been filled
1580}
particleIterator
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:60
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
CLHEP::Hep3Vector::mag
double mag() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
CLHEP::HepLorentzVector::setVect
void setVect(const Hep3Vector &)
G4InuclNuclei::getZ
G4int getZ() const
Definition: G4InuclNuclei.hh:138
G4InuclNuclei::getA
G4int getA() const
Definition: G4InuclNuclei.hh:137
G4InuclParticle::getKineticEnergy
G4double getKineticEnergy() const
Definition: G4InuclParticle.hh:114
G4InuclParticle::getMass
G4double getMass() const
Definition: G4InuclParticle.hh:106
G4InuclParticle::getEnergy
G4double getEnergy() const
Definition: G4InuclParticle.hh:118
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37

◆ isProjectile()

G4bool G4NucleiModel::isProjectile ( const G4CascadParticle &  cparticle) const

Definition at line 1071 of file G4NucleiModel.cc.

1071 {
1072 if (verboseLevel > 2) {
1073 G4cout << " isProjectile(cparticle): zone "
1074 << cparticle.getCurrentZone() << " of " << number_of_zones
1075 << " path " << cparticle.getCurrentPath()
1076 << " movingInside " << cparticle.movingInsideNuclei()
1077 << " nReflections " << cparticle.getNumberOfReflections()
1078 << G4endl;
1079 }
1080
1081 return (cparticle.getCurrentZone() == number_of_zones-1 && // At surface
1082 cparticle.getCurrentPath() == 1000. && // Input primary
1083 cparticle.getNumberOfReflections() <= 0 && // first pass
1084 (cparticle.movingInsideNuclei()||number_of_zones==1) // inbound
1085 );
1086}
G4CascadParticle::getCurrentPath
G4double getCurrentPath() const
Definition: G4CascadParticle.hh:93
G4CascadParticle::movingInsideNuclei
G4bool movingInsideNuclei() const
Definition: G4CascadParticle.hh:99
G4CascadParticle::getNumberOfReflections
G4int getNumberOfReflections() const
Definition: G4CascadParticle.hh:86

◆ printModel()

void G4NucleiModel::printModel ( ) const

Definition at line 540 of file G4NucleiModel.cc.

540 {
541 if (verboseLevel > 1) {
542 G4cout << " >>> G4NucleiModel::printModel" << G4endl;
543 }
544
545 G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
546 << " proton binding energy " << binding_energies[0]
547 << " neutron binding energy " << binding_energies[1] << G4endl
548 << " Nuclei radius " << nuclei_radius << " volume " << nuclei_volume
549 << " number of zones " << number_of_zones << G4endl;
550
551 for (G4int i = 0; i < number_of_zones; i++)
552 G4cout << " zone " << i+1 << " radius " << zone_radii[i]
553 << " volume " << zone_volumes[i] << G4endl
554 << " protons: density " << getDensity(1,i) << " PF " <<
555 getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
556 << " neutrons: density " << getDensity(2,i) << " PF " <<
557 getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
558 << " pions: VP " << getPotential(3,i) << G4endl;
559}
G4NucleiModel::getDensity
G4double getDensity(G4int ip, G4int izone) const
Definition: G4NucleiModel.hh:101
G4NucleiModel::getPotential
G4double getPotential(G4int ip, G4int izone) const
Definition: G4NucleiModel.hh:111

Referenced by generateModel().

◆ reset()

void G4NucleiModel::reset ( G4int  nHitNeutrons = 0,
G4int  nHitProtons = 0,
const std::vector< G4ThreeVector > *  hitPoints = 0 
)

Definition at line 224 of file G4NucleiModel.cc.

225 {
226 neutronNumberCurrent = neutronNumber - nHitNeutrons;
227 protonNumberCurrent = protonNumber - nHitProtons;
228
229 // zero or copy collision point array for trailing effect
230 if (!hitPoints || !hitPoints->empty()) collisionPts.clear();
231 else collisionPts = *hitPoints;
232}

Referenced by G4IntraNucleiCascader::copyWoundedNucleus(), generateModel(), and G4IntraNucleiCascader::newCascade().

◆ setVerboseLevel()

void G4NucleiModel::setVerboseLevel ( G4int  verbose)
inline

Definition at line 90 of file G4NucleiModel.hh.

90{ verboseLevel = verbose; }

Referenced by G4IntraNucleiCascader::setVerboseLevel().

◆ stillInside()

G4bool G4NucleiModel::stillInside ( const G4CascadParticle &  cparticle)
inline

Definition at line 145 of file G4NucleiModel.hh.

145 {
146 return cparticle.getCurrentZone() < number_of_zones;
147 }

Referenced by G4IntraNucleiCascader::generateCascade().

◆ totalCrossSection()

G4double G4NucleiModel::totalCrossSection ( G4double  ke,
G4int  rtype 
) const

Definition at line 1688 of file G4NucleiModel.cc.

1689{
1690 // All scattering cross-sections are available from tables
1691 const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(rtype);
1692 if (!xsecTable) {
1693 G4cerr << " unknown collison type = " << rtype << G4endl;
1694 return 0.;
1695 }
1696
1697 return (crossSectionUnits * xsecTable->getCrossSection(ke));
1698}
G4CascadeChannelTables::GetTable
static const G4CascadeChannel * GetTable(G4int initialState)
Definition: G4CascadeChannelTables.cc:55
G4CascadeChannel::getCrossSection
virtual G4double getCrossSection(double ke) const =0

◆ worthToPropagate()

G4bool G4NucleiModel::worthToPropagate ( const G4CascadParticle &  cparticle) const

Definition at line 1088 of file G4NucleiModel.cc.

1088 {
1089 if (verboseLevel > 1) {
1090 G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
1091 }
1092
1093 const G4double ekin_scale = 2.0;
1094
1095 G4bool worth = true;
1096
1097 if (cparticle.reflectedNow()) { // Just reflected -- keep going?
1098 G4int zone = cparticle.getCurrentZone();
1099 G4int ip = cparticle.getParticle().type();
1100
1101 // NOTE: Temporarily backing out use of potential for non-nucleons
1102 G4double ekin_cut = (cparticle.getParticle().nucleon()) ?
1103 getFermiKinetic(ip, zone) : 0.; //*** getPotential(ip, zone);
1104
1105 worth = cparticle.getParticle().getKineticEnergy()/ekin_scale > ekin_cut;
1106
1107 if (verboseLevel > 3) {
1108 G4cout << " type=" << ip
1109 << " ekin=" << cparticle.getParticle().getKineticEnergy()
1110 << " potential=" << ekin_cut
1111 << " : worth? " << worth << G4endl;
1112 }
1113 }
1114
1115 return worth;
1116}
G4CascadParticle::reflectedNow
G4bool reflectedNow() const
Definition: G4CascadParticle.hh:89
G4NucleiModel::getFermiKinetic
G4double getFermiKinetic(G4int ip, G4int izone) const
Definition: G4NucleiModel.cc:562

Referenced by G4IntraNucleiCascader::generateCascade().

The documentation for this class was generated from the following files: