#include <G4NucleiModel.hh>
Definition at line 82 of file G4NucleiModel.hh.
| typedef std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > G4NucleiModel::modelLists |
Definition at line 152 of file G4NucleiModel.hh.
| G4NucleiModel::G4NucleiModel | ( | ) |
Definition at line 194 of file G4NucleiModel.cc.
00195 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0), 00196 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0), 00197 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0), 00198 current_nucl2(0) {}
Definition at line 200 of file G4NucleiModel.cc.
References generateModel().
00201 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0), 00202 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0), 00203 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0), 00204 current_nucl2(0) { 00205 generateModel(a,z); 00206 }
| G4NucleiModel::G4NucleiModel | ( | G4InuclNuclei * | nuclei | ) | [explicit] |
Definition at line 208 of file G4NucleiModel.cc.
References generateModel(), and G4InuclParticleNames::nuclei.
00209 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0), 00210 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0), 00211 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0), 00212 current_nucl2(0) { 00213 generateModel(nuclei); 00214 }
| G4NucleiModel::~G4NucleiModel | ( | ) |
Definition at line 1646 of file G4NucleiModel.cc.
References G4cerr, G4cout, G4CascadeParameters::gammaQDScale(), G4CascadeInterpolator< NBINS >::interpolate(), photon, G4InuclParticleNames::pionMinus, G4InuclParticleNames::pionPlus, and G4InuclParticleNames::pionZero.
01646 { 01647 if (type != pionPlus && type != pionMinus && type != pionZero && 01648 type != photon) { 01649 G4cerr << " absorptionCrossSection only valid for incident pions" << G4endl; 01650 return 0.; 01651 } 01652 01653 G4double csec = 0.; 01654 01655 // Pion absorption is parametrized for low vs. medium energy 01656 if (type == pionPlus || type == pionMinus || type == pionZero) { 01657 if (ke < 0.3) csec = (0.1106 / std::sqrt(ke) - 0.8 01658 + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056) ); 01659 else if (ke < 1.0) csec = 3.6735 * (1.0-ke)*(1.0-ke); 01660 } 01661 01662 // Photon cross-section is binned from parametrization by Mi. Kossov 01663 if (type == photon) { 01664 static const G4double gammaQDscale = G4CascadeParameters::gammaQDScale(); 01665 static const G4double kebins[] = 01666 { 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1, 01667 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8, 01668 2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 }; 01669 static const G4double gammaD[] = 01670 { 2.8, 1.3, 0.89, 0.56, 0.38, 0.27, 0.19, 0.14, 0.098, 01671 0.071, 0.054, 0.0003, 0.0007, 0.0027, 0.0014, 0.001, 0.0012, 0.0005, 01672 0.0003, 0.0002,0.0002, 0.0002, 0.0002, 0.0002, 0.0001, 0.0001, 0.0001, 01673 0.0001, 0.0001, 0.0001 }; 01674 static G4CascadeInterpolator<30> interp(kebins); 01675 csec = interp.interpolate(ke, gammaD) * gammaQDscale; 01676 } 01677 01678 if (csec < 0.0) csec = 0.0; 01679 01680 if (verboseLevel > 2) { 01681 G4cout << " ekin " << ke << " abs. csec " << csec << " mb" << G4endl; 01682 } 01683 01684 return crossSectionUnits * csec; 01685 }
| G4bool G4NucleiModel::empty | ( | ) | const [inline] |
Definition at line 141 of file G4NucleiModel.hh.
Referenced by G4IntraNucleiCascader::generateCascade().
Definition at line 241 of file G4NucleiModel.cc.
References G4InuclSpecialFunctions::G4cbrt(), G4cout, G4endl, neutron, printModel(), G4InuclParticleNames::proton, and reset().
00241 { 00242 if (verboseLevel) { 00243 G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z 00244 << G4endl; 00245 } 00246 00247 // If model already built, just return; otherwise intialize everything 00248 if (a == A && z == Z) { 00249 if (verboseLevel > 1) G4cout << " model already generated" << z << G4endl; 00250 reset(); 00251 return; 00252 } 00253 00254 A = a; 00255 Z = z; 00256 delete theNucleus; 00257 theNucleus = new G4InuclNuclei(A,Z); // For conservation checking 00258 00259 neutronNumber = A - Z; 00260 protonNumber = Z; 00261 reset(); 00262 00263 if (verboseLevel > 3) { 00264 G4cout << " crossSectionUnits = " << crossSectionUnits << G4endl 00265 << " radiusUnits = " << radiusUnits << G4endl 00266 << " skinDepth = " << skinDepth << G4endl 00267 << " radiusScale = " << radiusScale << G4endl 00268 << " radiusScale2 = " << radiusScale2 << G4endl 00269 << " radiusForSmall = " << radiusForSmall << G4endl 00270 << " radScaleAlpha = " << radScaleAlpha << G4endl 00271 << " fermiMomentum = " << fermiMomentum << G4endl 00272 << " piTimes4thirds = " << piTimes4thirds << G4endl; 00273 } 00274 00275 G4double nuclearRadius; // Nuclear radius computed from A 00276 if (A>4) nuclearRadius = radiusScale*G4cbrt(A) + radiusScale2/G4cbrt(A); 00277 else nuclearRadius = radiusForSmall * (A==4 ? radScaleAlpha : 1.); 00278 00279 // This will be used to pre-allocate lots of arrays below 00280 number_of_zones = (A < 5) ? 1 : (A < 100) ? 3 : 6; 00281 00282 // Clear all parameters arrays for reloading 00283 binding_energies.clear(); 00284 nucleon_densities.clear(); 00285 zone_potentials.clear(); 00286 fermi_momenta.clear(); 00287 zone_radii.clear(); 00288 zone_volumes.clear(); 00289 00290 fillBindingEnergies(); 00291 fillZoneRadii(nuclearRadius); 00292 00293 G4double tot_vol = fillZoneVolumes(nuclearRadius); // Woods-Saxon integral 00294 00295 fillPotentials(proton, tot_vol); // Protons 00296 fillPotentials(neutron, tot_vol); // Neutrons 00297 00298 // Additional flat zone potentials for other hadrons 00299 const std::vector<G4double> vp(number_of_zones, (A>4)?pion_vp:pion_vp_small); 00300 const std::vector<G4double> kp(number_of_zones, kaon_vp); 00301 const std::vector<G4double> hp(number_of_zones, hyperon_vp); 00302 00303 zone_potentials.push_back(vp); 00304 zone_potentials.push_back(kp); 00305 zone_potentials.push_back(hp); 00306 00307 nuclei_radius = zone_radii.back(); 00308 nuclei_volume = std::accumulate(zone_volumes.begin(),zone_volumes.end(),0.); 00309 00310 if (verboseLevel > 3) printModel(); 00311 }
| void G4NucleiModel::generateModel | ( | G4InuclNuclei * | nuclei | ) |
Definition at line 237 of file G4NucleiModel.cc.
References G4InuclParticleNames::nuclei.
Referenced by G4NucleiModel(), and G4IntraNucleiCascader::initialize().
00237 { 00238 generateModel(nuclei->getA(), nuclei->getZ()); 00239 }
| G4InuclElementaryParticle G4NucleiModel::generateNucleon | ( | G4int | type, | |
| G4int | zone | |||
| ) | const |
Definition at line 585 of file G4NucleiModel.cc.
References G4cout, and generateNucleonMomentum().
00585 { 00586 if (verboseLevel > 1) { 00587 G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl; 00588 } 00589 00590 G4LorentzVector mom = generateNucleonMomentum(type, zone); 00591 return G4InuclElementaryParticle(mom, type); 00592 }
| G4LorentzVector G4NucleiModel::generateNucleonMomentum | ( | G4int | type, | |
| G4int | zone | |||
| ) | const |
Definition at line 576 of file G4NucleiModel.cc.
References G4InuclSpecialFunctions::G4cbrt(), G4InuclSpecialFunctions::generateWithRandomAngles(), getFermiMomentum(), G4InuclElementaryParticle::getParticleMass(), and G4InuclSpecialFunctions::inuclRndm().
Referenced by generateNucleon().
00576 { 00577 G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm()); 00578 G4double mass = G4InuclElementaryParticle::getParticleMass(type); 00579 00580 return generateWithRandomAngles(pmod, mass); 00581 }
| void G4NucleiModel::generateParticleFate | ( | G4CascadParticle & | cparticle, | |
| G4ElementaryParticleCollider * | theEPCollider, | |||
| std::vector< G4CascadParticle > & | cascade | |||
| ) |
Definition at line 790 of file G4NucleiModel.cc.
References G4CascadeCheckBalance::collide(), G4ElementaryParticleCollider::collide(), G4cerr, G4cout, G4InuclParticle::getCharge(), G4CascadParticle::getCurrentZone(), G4CollisionOutput::getOutgoingParticles(), G4CascadParticle::getParticle(), G4CascadParticle::getPosition(), G4CascadParticle::incrementCurrentPath(), G4InuclElementaryParticle::nucleon(), G4CollisionOutput::numberOfOutgoingParticles(), G4CascadeCheckBalance::okay(), G4CollisionOutput::printCollisionOutput(), G4CascadParticle::propagateAlongThePath(), G4InuclElementaryParticle::quasi_deutron(), G4CollisionOutput::reset(), G4VCascadeCollider::setVerboseLevel(), G4InuclElementaryParticle::type(), and G4CascadParticle::updatePosition().
Referenced by G4IntraNucleiCascader::generateCascade().
00792 { 00793 if (verboseLevel > 1) 00794 G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl; 00795 00796 if (verboseLevel > 2) G4cout << " cparticle: " << cparticle << G4endl; 00797 00798 // Create four-vector checking 00799 #ifdef G4CASCADE_CHECK_ECONS 00800 G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel"); // Second arg is in GeV 00801 balance.setVerboseLevel(verboseLevel); 00802 #endif 00803 00804 outgoing_cparticles.clear(); // Clear return buffer for this event 00805 generateInteractionPartners(cparticle); // Fills "thePartners" data 00806 00807 if (thePartners.empty()) { // smth. is wrong -> needs special treatment 00808 if (verboseLevel) 00809 G4cerr << " generateParticleFate-> got empty interaction-partners list " 00810 << G4endl; 00811 return; 00812 } 00813 00814 G4int npart = thePartners.size(); // Last item is a total-path placeholder 00815 00816 if (npart == 1) { // cparticle is on the next zone entry 00817 cparticle.propagateAlongThePath(thePartners[0].second); 00818 cparticle.incrementCurrentPath(thePartners[0].second); 00819 boundaryTransition(cparticle); 00820 outgoing_cparticles.push_back(cparticle); 00821 00822 if (verboseLevel > 2) G4cout << " next zone \n" << cparticle << G4endl; 00823 } else { // there are possible interactions 00824 if (verboseLevel > 1) 00825 G4cout << " processing " << npart-1 << " possible interactions" << G4endl; 00826 00827 G4ThreeVector old_position = cparticle.getPosition(); 00828 G4InuclElementaryParticle& bullet = cparticle.getParticle(); 00829 G4bool no_interaction = true; 00830 G4int zone = cparticle.getCurrentZone(); 00831 00832 for (G4int i=0; i<npart-1; i++) { // Last item is a total-path placeholder 00833 if (i > 0) cparticle.updatePosition(old_position); 00834 00835 G4InuclElementaryParticle& target = thePartners[i].first; 00836 00837 if (verboseLevel > 3) { 00838 if (target.quasi_deutron()) G4cout << " try absorption: "; 00839 G4cout << " target " << target.type() << " bullet " << bullet.type() 00840 << G4endl; 00841 } 00842 00843 EPCoutput.reset(); 00844 theEPCollider->collide(&bullet, &target, EPCoutput); 00845 00846 // If collision failed, exit loop over partners 00847 if (EPCoutput.numberOfOutgoingParticles() == 0) break; 00848 00849 if (verboseLevel > 2) { 00850 EPCoutput.printCollisionOutput(); 00851 #ifdef G4CASCADE_CHECK_ECONS 00852 balance.collide(&bullet, &target, EPCoutput); 00853 balance.okay(); // Do checks, but ignore result 00854 #endif 00855 } 00856 00857 // Get list of outgoing particles for evaluation 00858 std::vector<G4InuclElementaryParticle>& outgoing_particles = 00859 EPCoutput.getOutgoingParticles(); 00860 00861 if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails 00862 00863 // Trailing effect: reject interaction at previously hit nucleon 00864 cparticle.propagateAlongThePath(thePartners[i].second); 00865 G4ThreeVector new_position = cparticle.getPosition(); 00866 00867 if (!passTrailing(new_position)) continue; 00868 collisionPts.push_back(new_position); 00869 00870 // Sort particles according to beta (fastest first) 00871 std::sort(outgoing_particles.begin(), outgoing_particles.end(), 00872 G4ParticleLargerBeta() ); 00873 00874 if (verboseLevel > 2) 00875 G4cout << " adding " << outgoing_particles.size() 00876 << " output particles" << G4endl; 00877 00878 // NOTE: Embedded temporary is optimized away (no copying gets done) 00879 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) { 00880 outgoing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip], 00881 new_position, zone, 00882 0.0, 0)); 00883 } 00884 00885 no_interaction = false; 00886 current_nucl1 = 0; 00887 current_nucl2 = 0; 00888 00889 #ifdef G4CASCADE_DEBUG_CHARGE 00890 { 00891 G4double out_charge = 0.0; 00892 00893 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) 00894 out_charge += outgoing_particles[ip].getCharge(); 00895 00896 G4cout << " multiplicity " << outgoing_particles.size() 00897 << " bul type " << bullet.type() 00898 << " targ type " << target.type() 00899 << "\n initial charge " 00900 << bullet.getCharge() + target.getCharge() 00901 << " out charge " << out_charge << G4endl; 00902 } 00903 #endif 00904 00905 if (verboseLevel > 2) 00906 G4cout << " partner type " << target.type() << G4endl; 00907 00908 if (target.nucleon()) { 00909 current_nucl1 = target.type(); 00910 } else { 00911 if (verboseLevel > 2) G4cout << " good absorption " << G4endl; 00912 current_nucl1 = (target.type() - 100) / 10; 00913 current_nucl2 = target.type() - 100 - 10 * current_nucl1; 00914 } 00915 00916 if (current_nucl1 == 1) { 00917 if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl; 00918 protonNumberCurrent--; 00919 } else { 00920 if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl; 00921 neutronNumberCurrent--; 00922 } 00923 00924 if (current_nucl2 == 1) { 00925 if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl; 00926 protonNumberCurrent--; 00927 } else if (current_nucl2 == 2) { 00928 if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl; 00929 neutronNumberCurrent--; 00930 } 00931 00932 break; 00933 } // loop over partners 00934 00935 if (no_interaction) { // still no interactions 00936 if (verboseLevel > 1) G4cout << " no interaction " << G4endl; 00937 00938 // For conservation checking (below), get particle before updating 00939 static G4InuclElementaryParticle prescatCP; // Avoid memory churn 00940 prescatCP = cparticle.getParticle(); 00941 00942 // Last "partner" is just a total-path placeholder 00943 cparticle.updatePosition(old_position); 00944 cparticle.propagateAlongThePath(thePartners[npart-1].second); 00945 cparticle.incrementCurrentPath(thePartners[npart-1].second); 00946 boundaryTransition(cparticle); 00947 outgoing_cparticles.push_back(cparticle); 00948 00949 // Check conservation for simple scattering (ignore target nucleus!) 00950 #ifdef G4CASCADE_CHECK_ECONS 00951 if (verboseLevel > 2) { 00952 balance.collide(&prescatCP, 0, outgoing_cparticles); 00953 balance.okay(); // Report violations, but don't act on them 00954 } 00955 #endif 00956 } // if (no_interaction) 00957 } // if (npart == 1) [else] 00958 00959 return; 00960 }
Definition at line 562 of file G4NucleiModel.cc.
References G4InuclElementaryParticle::getParticleMass().
Referenced by worthToPropagate().
00562 { 00563 G4double ekin = 0.0; 00564 00565 if (ip < 3 && izone < number_of_zones) { // ip for proton/neutron only 00566 G4double pfermi = fermi_momenta[ip - 1][izone]; 00567 G4double mass = G4InuclElementaryParticle::getParticleMass(ip); 00568 00569 ekin = std::sqrt(pfermi*pfermi + mass*mass) - mass; 00570 } 00571 return ekin; 00572 }
Definition at line 105 of file G4NucleiModel.hh.
Referenced by generateNucleonMomentum(), and printModel().
| G4int G4NucleiModel::getNumberOfNeutrons | ( | ) | const [inline] |
Definition at line 138 of file G4NucleiModel.hh.
Referenced by G4IntraNucleiCascader::generateCascade().
| G4int G4NucleiModel::getNumberOfProtons | ( | ) | const [inline] |
Definition at line 139 of file G4NucleiModel.hh.
Referenced by G4IntraNucleiCascader::generateCascade().
| G4int G4NucleiModel::getNumberOfZones | ( | ) | const [inline] |
Definition at line 111 of file G4NucleiModel.hh.
Referenced by printModel().
00111 { 00112 if (ip == 9) return 0.0; // Special case for photons 00113 G4int ip0 = ip < 3 ? ip - 1 : 2; 00114 if (ip > 10 && ip < 18) ip0 = 3; 00115 if (ip > 20) ip0 = 4; 00116 return izone < number_of_zones ? zone_potentials[ip0][izone] : 0.0; 00117 }
Definition at line 123 of file G4NucleiModel.hh.
00123 { 00124 return ( (izone<0) ? 0. 00125 : (izone<number_of_zones) ? zone_radii[izone] : nuclei_radius); 00126 }
| G4double G4NucleiModel::getRadius | ( | ) | const [inline] |
| G4double G4NucleiModel::getRadiusUnits | ( | ) | const [inline] |
Definition at line 120 of file G4NucleiModel.hh.
Referenced by G4IntraNucleiCascader::processSecondary().
Definition at line 158 of file G4NucleiModel.hh.
Referenced by G4IntraNucleiCascader::generateCascade().
Definition at line 127 of file G4NucleiModel.hh.
00127 { 00128 return ( (izone<0) ? 0. 00129 : (izone<number_of_zones) ? zone_volumes[izone] : nuclei_volume); 00130 }
Definition at line 133 of file G4NucleiModel.hh.
Referenced by G4IntraNucleiCascader::processSecondary().
00133 { 00134 for (G4int iz=0; iz<number_of_zones; iz++) if (r<zone_radii[iz]) return iz; 00135 return number_of_zones; 00136 }
| void G4NucleiModel::initializeCascad | ( | G4InuclNuclei * | bullet, | |
| G4InuclNuclei * | target, | |||
| modelLists & | output | |||
| ) |
Definition at line 1184 of file G4NucleiModel.cc.
References G4LorentzConvertor::backToTheLab(), G4InuclSpecialFunctions::bindingEnergy(), G4cout, G4InuclSpecialFunctions::generateWithRandomAngles(), G4InuclNuclei::getA(), G4InuclParticle::getEnergy(), G4InuclParticle::getKineticEnergy(), G4InuclParticle::getMass(), G4InuclNuclei::getZ(), G4InuclSpecialFunctions::inuclRndm(), G4InuclSpecialFunctions::randomPHI(), and G4LorentzConvertor::toTheTargetRestFrame().
01186 { 01187 if (verboseLevel) { 01188 G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)" 01189 << G4endl; 01190 } 01191 01192 const G4double large = 1000.0; 01193 const G4double max_a_for_cascad = 5.0; 01194 const G4double ekin_cut = 2.0; 01195 const G4double small_ekin = 1.0e-6; 01196 const G4double r_large2for3 = 62.0; 01197 const G4double r0forAeq3 = 3.92; 01198 const G4double s3max = 6.5; 01199 const G4double r_large2for4 = 69.14; 01200 const G4double r0forAeq4 = 4.16; 01201 const G4double s4max = 7.0; 01202 const G4int itry_max = 100; 01203 01204 // Convenient references to input buffer contents 01205 std::vector<G4CascadParticle>& casparticles = output.first; 01206 std::vector<G4InuclElementaryParticle>& particles = output.second; 01207 01208 casparticles.clear(); // Reset input buffer to fill this event 01209 particles.clear(); 01210 01211 // first decide whether it will be cascad or compound final nuclei 01212 G4int ab = bullet->getA(); 01213 G4int zb = bullet->getZ(); 01214 G4int at = target->getA(); 01215 G4int zt = target->getZ(); 01216 01217 G4double massb = bullet->getMass(); // For creating LorentzVectors below 01218 01219 if (ab < max_a_for_cascad) { 01220 01221 G4double benb = bindingEnergy(ab,zb)/GeV / G4double(ab); 01222 G4double bent = bindingEnergy(at,zt)/GeV / G4double(at); 01223 G4double ben = benb < bent ? bent : benb; 01224 01225 if (bullet->getKineticEnergy()/ab > ekin_cut*ben) { 01226 G4int itryg = 0; 01227 01228 while (casparticles.size() == 0 && itryg < itry_max) { 01229 itryg++; 01230 particles.clear(); 01231 01232 // nucleons coordinates and momenta in nuclei rest frame 01233 coordinates.clear(); 01234 momentums.clear(); 01235 01236 if (ab < 3) { // deuteron, simplest case 01237 G4double r = 2.214 - 3.4208 * std::log(1.0 - 0.981 * inuclRndm()); 01238 G4ThreeVector coord1 = generateWithRandomAngles(r).vect(); 01239 coordinates.push_back(coord1); 01240 coordinates.push_back(-coord1); 01241 01242 G4double p = 0.0; 01243 G4bool bad = true; 01244 G4int itry = 0; 01245 01246 while (bad && itry < itry_max) { 01247 itry++; 01248 p = 456.0 * inuclRndm(); 01249 01250 if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() && 01251 p * r > 312.0) bad = false; 01252 } 01253 01254 if (itry == itry_max) 01255 if (verboseLevel > 2){ 01256 G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl; 01257 } 01258 01259 p = 0.0005 * p; 01260 01261 if (verboseLevel > 2){ 01262 G4cout << " p nuc " << p << G4endl; 01263 } 01264 01265 G4LorentzVector mom = generateWithRandomAngles(p, massb); 01266 01267 momentums.push_back(mom); 01268 mom.setVect(-mom.vect()); 01269 momentums.push_back(-mom); 01270 } else { 01271 G4ThreeVector coord1; 01272 01273 G4bool badco = true; 01274 01275 G4int itry = 0; 01276 01277 if (ab == 3) { 01278 while (badco && itry < itry_max) { 01279 if (itry > 0) coordinates.clear(); 01280 itry++; 01281 G4int i(0); 01282 01283 for (i = 0; i < 2; i++) { 01284 G4int itry1 = 0; 01285 G4double ss, u, rho; 01286 G4double fmax = std::exp(-0.5) / std::sqrt(0.5); 01287 01288 while (itry1 < itry_max) { 01289 itry1++; 01290 ss = -std::log(inuclRndm()); 01291 u = fmax * inuclRndm(); 01292 rho = std::sqrt(ss) * std::exp(-ss); 01293 01294 if (rho > u && ss < s3max) { 01295 ss = r0forAeq3 * std::sqrt(ss); 01296 coord1 = generateWithRandomAngles(ss).vect(); 01297 coordinates.push_back(coord1); 01298 01299 if (verboseLevel > 2){ 01300 G4cout << " i " << i << " r " << coord1.mag() << G4endl; 01301 } 01302 break; 01303 } 01304 } 01305 01306 if (itry1 == itry_max) { // bad case 01307 coord1.set(10000.,10000.,10000.); 01308 coordinates.push_back(coord1); 01309 break; 01310 } 01311 } 01312 01313 coord1 = -coordinates[0] - coordinates[1]; 01314 if (verboseLevel > 2) { 01315 G4cout << " 3 r " << coord1.mag() << G4endl; 01316 } 01317 01318 coordinates.push_back(coord1); 01319 01320 G4bool large_dist = false; 01321 01322 for (i = 0; i < 2; i++) { 01323 for (G4int j = i+1; j < 3; j++) { 01324 G4double r2 = (coordinates[i]-coordinates[j]).mag2(); 01325 01326 if (verboseLevel > 2) { 01327 G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl; 01328 } 01329 01330 if (r2 > r_large2for3) { 01331 large_dist = true; 01332 01333 break; 01334 } 01335 } 01336 01337 if (large_dist) break; 01338 } 01339 01340 if(!large_dist) badco = false; 01341 01342 } 01343 01344 } else { // a >= 4 01345 G4double b = 3./(ab - 2.0); 01346 G4double b1 = 1.0 - b / 2.0; 01347 G4double u = b1 + std::sqrt(b1 * b1 + b); 01348 G4double fmax = (1.0 + u/b) * u * std::exp(-u); 01349 01350 while (badco && itry < itry_max) { 01351 01352 if (itry > 0) coordinates.clear(); 01353 itry++; 01354 G4int i(0); 01355 01356 for (i = 0; i < ab-1; i++) { 01357 G4int itry1 = 0; 01358 G4double ss; 01359 01360 while (itry1 < itry_max) { 01361 itry1++; 01362 ss = -std::log(inuclRndm()); 01363 u = fmax * inuclRndm(); 01364 01365 if (std::sqrt(ss) * std::exp(-ss) * (1.0 + ss/b) > u 01366 && ss < s4max) { 01367 ss = r0forAeq4 * std::sqrt(ss); 01368 coord1 = generateWithRandomAngles(ss).vect(); 01369 coordinates.push_back(coord1); 01370 01371 if (verboseLevel > 2) { 01372 G4cout << " i " << i << " r " << coord1.mag() << G4endl; 01373 } 01374 01375 break; 01376 } 01377 } 01378 01379 if (itry1 == itry_max) { // bad case 01380 coord1.set(10000.,10000.,10000.); 01381 coordinates.push_back(coord1); 01382 break; 01383 } 01384 } 01385 01386 coord1 *= 0.0; // Cheap way to reset 01387 for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j]; 01388 01389 coordinates.push_back(coord1); 01390 01391 if (verboseLevel > 2){ 01392 G4cout << " last r " << coord1.mag() << G4endl; 01393 } 01394 01395 G4bool large_dist = false; 01396 01397 for (i = 0; i < ab-1; i++) { 01398 for (G4int j = i+1; j < ab; j++) { 01399 01400 G4double r2 = (coordinates[i]-coordinates[j]).mag2(); 01401 01402 if (verboseLevel > 2){ 01403 G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl; 01404 } 01405 01406 if (r2 > r_large2for4) { 01407 large_dist = true; 01408 01409 break; 01410 } 01411 } 01412 01413 if (large_dist) break; 01414 } 01415 01416 if (!large_dist) badco = false; 01417 } 01418 } 01419 01420 if(badco) { 01421 G4cout << " can not generate the nucleons coordinates for a " 01422 << ab << G4endl; 01423 01424 casparticles.clear(); // Return empty buffer on error 01425 particles.clear(); 01426 return; 01427 01428 } else { // momentums 01429 G4double p, u, x; 01430 G4LorentzVector mom; 01431 //G4bool badp = True; 01432 01433 for (G4int i = 0; i < ab - 1; i++) { 01434 G4int itry2 = 0; 01435 01436 while(itry2 < itry_max) { 01437 itry2++; 01438 u = -std::log(0.879853 - 0.8798502 * inuclRndm()); 01439 x = u * std::exp(-u); 01440 01441 if(x > inuclRndm()) { 01442 p = std::sqrt(0.01953 * u); 01443 mom = generateWithRandomAngles(p, massb); 01444 momentums.push_back(mom); 01445 01446 break; 01447 } 01448 } // while (itry2 < itry_max) 01449 01450 if(itry2 == itry_max) { 01451 G4cout << " can not generate proper momentum for a " 01452 << ab << G4endl; 01453 01454 casparticles.clear(); // Return empty buffer on error 01455 particles.clear(); 01456 return; 01457 } 01458 } // for (i=0 ... 01459 // last momentum 01460 01461 mom *= 0.; // Cheap way to reset 01462 mom.setE(bullet->getEnergy()+target->getEnergy()); 01463 01464 for(G4int j=0; j< ab-1; j++) mom -= momentums[j]; 01465 01466 momentums.push_back(mom); 01467 } 01468 } 01469 01470 // Coordinates and momenta at rest are generated, now back to the lab 01471 G4double rb = 0.0; 01472 G4int i(0); 01473 01474 for(i = 0; i < G4int(coordinates.size()); i++) { 01475 G4double rp = coordinates[i].mag2(); 01476 01477 if(rp > rb) rb = rp; 01478 } 01479 01480 // nuclei i.p. as a whole 01481 G4double s1 = std::sqrt(inuclRndm()); 01482 G4double phi = randomPHI(); 01483 G4double rz = (nuclei_radius + rb) * s1; 01484 G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi), 01485 -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1)); 01486 01487 for (i = 0; i < G4int(coordinates.size()); i++) { 01488 coordinates[i] += global_pos; 01489 } 01490 01491 // all nucleons at rest 01492 raw_particles.clear(); 01493 01494 for (G4int ipa = 0; ipa < ab; ipa++) { 01495 G4int knd = ipa < zb ? 1 : 2; 01496 raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd)); 01497 } 01498 01499 G4InuclElementaryParticle dummy(small_ekin, 1); 01500 G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet); 01501 toTheBulletRestFrame.toTheTargetRestFrame(); 01502 01503 particleIterator ipart; 01504 01505 for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) { 01506 ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum())); 01507 } 01508 01509 // fill cascad particles and outgoing particles 01510 01511 for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) { 01512 G4LorentzVector mom = raw_particles[ip].getMomentum(); 01513 G4double pmod = mom.rho(); 01514 G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod; 01515 G4double det = t0 * t0 + nuclei_radius * nuclei_radius 01516 - coordinates[ip].mag2(); 01517 G4double tr = -1.0; 01518 01519 if(det > 0.0) { 01520 G4double t1 = t0 + std::sqrt(det); 01521 G4double t2 = t0 - std::sqrt(det); 01522 01523 if(std::fabs(t1) <= std::fabs(t2)) { 01524 if(t1 > 0.0) { 01525 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1; 01526 } 01527 01528 if(tr < 0.0 && t2 > 0.0) { 01529 01530 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2; 01531 } 01532 01533 } else { 01534 if(t2 > 0.0) { 01535 01536 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2; 01537 } 01538 01539 if(tr < 0.0 && t1 > 0.0) { 01540 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1; 01541 } 01542 } 01543 01544 } 01545 01546 if(tr >= 0.0) { // cascad particle 01547 coordinates[ip] += mom.vect()*tr / pmod; 01548 casparticles.push_back(G4CascadParticle(raw_particles[ip], 01549 coordinates[ip], 01550 number_of_zones, large, 0)); 01551 01552 } else { 01553 particles.push_back(raw_particles[ip]); 01554 } 01555 } 01556 } 01557 01558 if(casparticles.size() == 0) { 01559 particles.clear(); 01560 01561 G4cout << " can not generate proper distribution for " << itry_max 01562 << " steps " << G4endl; 01563 } 01564 } 01565 } 01566 01567 if(verboseLevel > 2){ 01568 G4int ip(0); 01569 01570 G4cout << " cascad particles: " << casparticles.size() << G4endl; 01571 for(ip = 0; ip < G4int(casparticles.size()); ip++) 01572 G4cout << casparticles[ip] << G4endl; 01573 01574 G4cout << " outgoing particles: " << particles.size() << G4endl; 01575 for(ip = 0; ip < G4int(particles.size()); ip++) 01576 G4cout << particles[ip] << G4endl; 01577 } 01578 01579 return; // Buffer has been filled 01580 }
| G4CascadParticle G4NucleiModel::initializeCascad | ( | G4InuclElementaryParticle * | particle | ) |
Definition at line 1164 of file G4NucleiModel.cc.
References G4cout, G4InuclSpecialFunctions::generateWithFixedTheta(), and G4InuclSpecialFunctions::inuclRndm().
Referenced by G4IntraNucleiCascader::setupCascade().
01164 { 01165 if (verboseLevel > 1) { 01166 G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl; 01167 } 01168 01169 const G4double large = 1000.0; 01170 01171 // FIXME: Previous version generated random sin(theta), then used -cos(theta) 01172 // Using generateWithRandomAngles changes result! 01173 // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect(); 01174 G4double costh = std::sqrt(1.0 - inuclRndm()); 01175 G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius); 01176 01177 G4CascadParticle cpart(*particle, pos, number_of_zones, large, 0); 01178 01179 if (verboseLevel > 2) G4cout << cpart << G4endl; 01180 01181 return cpart; 01182 }
| G4bool G4NucleiModel::isProjectile | ( | const G4CascadParticle & | cparticle | ) | const |
Definition at line 1071 of file G4NucleiModel.cc.
References G4cout, G4CascadParticle::getCurrentPath(), G4CascadParticle::getCurrentZone(), G4CascadParticle::getNumberOfReflections(), and G4CascadParticle::movingInsideNuclei().
01071 { 01072 if (verboseLevel > 2) { 01073 G4cout << " isProjectile(cparticle): zone " 01074 << cparticle.getCurrentZone() << " of " << number_of_zones 01075 << " path " << cparticle.getCurrentPath() 01076 << " movingInside " << cparticle.movingInsideNuclei() 01077 << " nReflections " << cparticle.getNumberOfReflections() 01078 << G4endl; 01079 } 01080 01081 return (cparticle.getCurrentZone() == number_of_zones-1 && // At surface 01082 cparticle.getCurrentPath() == 1000. && // Input primary 01083 cparticle.getNumberOfReflections() <= 0 && // first pass 01084 (cparticle.movingInsideNuclei()||number_of_zones==1) // inbound 01085 ); 01086 }
| void G4NucleiModel::printModel | ( | ) | const |
Definition at line 540 of file G4NucleiModel.cc.
References G4cout, getDensity(), getFermiMomentum(), and getPotential().
Referenced by generateModel().
00540 { 00541 if (verboseLevel > 1) { 00542 G4cout << " >>> G4NucleiModel::printModel" << G4endl; 00543 } 00544 00545 G4cout << " nuclei model for A " << A << " Z " << Z << G4endl 00546 << " proton binding energy " << binding_energies[0] 00547 << " neutron binding energy " << binding_energies[1] << G4endl 00548 << " Nuclei radius " << nuclei_radius << " volume " << nuclei_volume 00549 << " number of zones " << number_of_zones << G4endl; 00550 00551 for (G4int i = 0; i < number_of_zones; i++) 00552 G4cout << " zone " << i+1 << " radius " << zone_radii[i] 00553 << " volume " << zone_volumes[i] << G4endl 00554 << " protons: density " << getDensity(1,i) << " PF " << 00555 getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl 00556 << " neutrons: density " << getDensity(2,i) << " PF " << 00557 getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl 00558 << " pions: VP " << getPotential(3,i) << G4endl; 00559 }
| void G4NucleiModel::reset | ( | G4int | nHitNeutrons = 0, |
|
| G4int | nHitProtons = 0, |
|||
| const std::vector< G4ThreeVector > * | hitPoints = 0 | |||
| ) |
Definition at line 224 of file G4NucleiModel.cc.
Referenced by G4IntraNucleiCascader::copyWoundedNucleus(), generateModel(), and G4IntraNucleiCascader::newCascade().
00225 { 00226 neutronNumberCurrent = neutronNumber - nHitNeutrons; 00227 protonNumberCurrent = protonNumber - nHitProtons; 00228 00229 // zero or copy collision point array for trailing effect 00230 if (!hitPoints || !hitPoints->empty()) collisionPts.clear(); 00231 else collisionPts = *hitPoints; 00232 }
| void G4NucleiModel::setVerboseLevel | ( | G4int | verbose | ) | [inline] |
Definition at line 90 of file G4NucleiModel.hh.
Referenced by G4IntraNucleiCascader::setVerboseLevel().
| G4bool G4NucleiModel::stillInside | ( | const G4CascadParticle & | cparticle | ) | [inline] |
Definition at line 145 of file G4NucleiModel.hh.
References G4CascadParticle::getCurrentZone().
Referenced by G4IntraNucleiCascader::generateCascade().
00145 { 00146 return cparticle.getCurrentZone() < number_of_zones; 00147 }
Definition at line 1688 of file G4NucleiModel.cc.
References G4cerr, G4CascadeChannel::getCrossSection(), and G4CascadeChannelTables::GetTable().
01689 { 01690 // All scattering cross-sections are available from tables 01691 const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(rtype); 01692 if (!xsecTable) { 01693 G4cerr << " unknown collison type = " << rtype << G4endl; 01694 return 0.; 01695 } 01696 01697 return (crossSectionUnits * xsecTable->getCrossSection(ke)); 01698 }
| G4bool G4NucleiModel::worthToPropagate | ( | const G4CascadParticle & | cparticle | ) | const |
Definition at line 1088 of file G4NucleiModel.cc.
References G4cout, G4CascadParticle::getCurrentZone(), getFermiKinetic(), G4InuclParticle::getKineticEnergy(), G4CascadParticle::getParticle(), G4InuclElementaryParticle::nucleon(), G4CascadParticle::reflectedNow(), and G4InuclElementaryParticle::type().
Referenced by G4IntraNucleiCascader::generateCascade().
01088 { 01089 if (verboseLevel > 1) { 01090 G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl; 01091 } 01092 01093 const G4double ekin_scale = 2.0; 01094 01095 G4bool worth = true; 01096 01097 if (cparticle.reflectedNow()) { // Just reflected -- keep going? 01098 G4int zone = cparticle.getCurrentZone(); 01099 G4int ip = cparticle.getParticle().type(); 01100 01101 // NOTE: Temporarily backing out use of potential for non-nucleons 01102 G4double ekin_cut = (cparticle.getParticle().nucleon()) ? 01103 getFermiKinetic(ip, zone) : 0.; //*** getPotential(ip, zone); 01104 01105 worth = cparticle.getParticle().getKineticEnergy()/ekin_scale > ekin_cut; 01106 01107 if (verboseLevel > 3) { 01108 G4cout << " type=" << ip 01109 << " ekin=" << cparticle.getParticle().getKineticEnergy() 01110 << " potential=" << ekin_cut 01111 << " : worth? " << worth << G4endl; 01112 } 01113 } 01114 01115 return worth; 01116 }
1.4.7