92 "G4ScreenedNuclearRecoil.cc,v 1.57 2008/05/07 11:51:26 marcus Exp GEANT4 tag ";
135 0, 1.007940, 4.002602, 6.941000, 9.012182, 10.811000, 12.010700,
136 14.006700, 15.999400, 18.998403, 20.179700, 22.989770, 24.305000, 26.981538,
138 30.973761, 32.065000, 35.453000, 39.948000, 39.098300, 40.078000, 44.955910,
140 50.941500, 51.996100, 54.938049, 55.845000, 58.933200, 58.693400, 63.546000,
142 69.723000, 72.640000, 74.921600, 78.960000, 79.904000, 83.798000, 85.467800,
144 88.905850, 91.224000, 92.906380, 95.940000, 98.000000, 101.070000, 102.905500,
146 107.868200, 112.411000, 114.818000, 118.710000, 121.760000, 127.600000,
147 126.904470, 131.293000,
148 132.905450, 137.327000, 138.905500, 140.116000, 140.907650, 144.240000,
149 145.000000, 150.360000,
150 151.964000, 157.250000, 158.925340, 162.500000, 164.930320, 167.259000,
151 168.934210, 173.040000,
152 174.967000, 178.490000, 180.947900, 183.840000, 186.207000, 190.230000,
153 192.217000, 195.078000,
154 196.966550, 200.590000, 204.383300, 207.200000, 208.980380, 209.000000,
155 210.000000, 222.000000,
156 223.000000, 226.000000, 227.000000, 232.038100, 231.035880, 238.028910,
157 237.000000, 244.000000,
158 243.000000, 247.000000, 247.000000, 251.000000, 252.000000, 257.000000,
159 258.000000, 259.000000,
160 262.000000, 261.000000, 262.000000, 266.000000, 264.000000, 277.000000,
161 268.000000, 281.000000,
162 272.000000, 285.000000, 282.500000, 289.000000, 287.500000, 292.000000};
177 if (nMatElements == 1)
179 element= (*elementVector)[0];
188 for (
G4int k=0 ;
k < nMatElements ;
k++ )
190 nsum+=atomDensities[
k];
191 element= (*elementVector)[
k];
192 if (nsum >= random)
break;
211 for(
G4int i=0; i<nIsotopes; i++) {
214 if(asum >= random)
break;
218 N=(
G4int)std::floor(element->
GetN()+0.5);
226 for(i=0; i<nIsotopes; i++) {
228 N=(*isoV)[i]->GetN();
229 if(asum >= random)
break;
236 ParticleCache::iterator
p=
targetMap.find(Z*1000+N);
248 const G4int nmfpvals=200;
250 std::vector<G4double> evals(nmfpvals), mfpvals(nmfpvals);
254 if (materialTable == 0) {
return; }
272 for (
G4int kel=0 ; kel < nMatElements ; kel++ )
274 element=elementVector[kel];
286 for (
G4int i=1; i<nmfpvals-1; i++) evals[i]=emin*std::exp(logint*i);
291 for (
G4int eidx=0; eidx < nmfpvals; eidx++) mfpvals[eidx] = 0.0;
295 for (
G4int kel=0 ; kel < nMatElements ; kel++ )
297 element=elementVector[kel];
300 G4double ndens = atomDensities[kel];
303 for (
G4int eidx=0; eidx < nmfpvals; eidx++) {
304 mfpvals[eidx] += ndens*sigma(evals[eidx]);
309 for (
G4int eidx=0; eidx < nmfpvals; eidx++) {
310 mfpvals[eidx] = 1.0/mfpvals[eidx];
314 mfpvals,
true,0,
true,0);
324 screeningKey(ScreeningKey),
325 generateRecoils(GenerateRecoils), avoidReactions(1),
326 recoilCutoff(RecoilCutoff), physicsCutoff(PhysicsCutoff),
327 hardeningFraction(0.0), hardeningFactor(1.0),
328 externalCrossSectionConstructor(0),
350 std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xt=
400 std::pow(a1,1.0/3.0)) + 1.4)*
fermi);
439 std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xh=
447 }
else xs=(*xh).second;
500 G4double lattice=0.5/std::pow(numberDensity,1.0/3.0);
515 if(sigopi < lattice*lattice) {
536 printf(
"ScreenedNuclear impact reject: length=%.3f P=%.4f limit=%.4f\n",
562 std::vector<G4ScreenedCollisionStage *>::iterator stage=
566 (*stage)->DoCollisionStep(
this,aTrack, aStep);
584 phifunc(c2.const_plugin_function()),
585 xovereps(c2.linear(0., 0., 0.)),
587 diff(c2.quadratic(0., 0., 0., 1.)-xovereps*phifunc)
603 G4double mlrho4=((((3.517e-4*y+1.401e-2)*y+2.393
e-1)*y+2.734)*y+2.220);
606 xx0=std::sqrt(bb2+std::sqrt(bb2*bb2+rho4));
609 xx0=ee+std::sqrt(ee*ee+beta*beta);
624 &root_error, &phip, &phip2)/au;
627 G4cout <<
"Screened Coulomb Root Finder Error" <<
G4endl;
628 G4cout <<
"au " << au <<
" A " << A <<
" a1 " << a1
629 <<
" xx1 " << xx1 <<
" eps " << eps
630 <<
" beta " << beta <<
G4endl;
637 <<
" target " << beta*beta*au*au ;
647 G4double lambda0=1.0/std::sqrt(0.5+beta*beta/(2.0*xx1*xx1)
648 -phiprime/(2.0*eps));
655 G4double xvals[]={0.98302349, 0.84652241, 0.53235309, 0.18347974};
656 G4double weights[]={0.03472124, 0.14769029, 0.23485003, 0.18602489};
660 ff=1.0/std::sqrt(1.0-
phifunc(x*au)/(x*eps)-beta*beta/(x*x));
661 alpha+=weights[
k]*
ff;
669 G4double sintheta=std::sin(thetac1);
670 G4double costheta=-std::cos(thetac1);
679 G4double zeta=std::atan2(sintheta, 1-costheta);
716 G4double Ec = incidentEnergy*(A/(A+a1));
736 if(incidentEnergy-eRecoil < master->GetRecoilCutoff()*a1) {
739 incidentEnergy-eRecoil);
771 recoilMomentumDirection=
772 recoilMomentumDirection.
rotateUz(incidentDirection);
774 recoilMomentumDirection*std::sqrt(2.0*eRecoil*kin.
a2*
amu_c2);
787 recoilMomentumDirection,eRecoil) ;
807 if(nam ==
"GenericIon" || nam ==
"proton"
808 || nam ==
"deuteron" || nam ==
"triton"
809 || nam ==
"alpha" || nam ==
"He3") {
837 static const size_t ncoef=4;
838 static G4double scales[ncoef]={-3.2, -0.9432, -0.4028, -0.2016};
839 static G4double coefs[ncoef]={0.1818,0.5099,0.2802,0.0281};
842 0.8854*
angstrom*0.529/(std::pow(z1, 0.23)+std::pow(z2,0.23));
843 std::vector<G4double>
r(npoints),
phi(npoints);
845 for(
size_t i=0; i<npoints; i++) {
846 G4double rr=(float)i/(
float)(npoints-1);
850 for(
size_t j=0; j<ncoef; j++)
851 sum+=coefs[j]*std::exp(scales[j]*
r[i]/au);
857 for(
size_t j=0; j<ncoef; j++)
858 phiprime0+=scales[j]*coefs[j]*std::exp(scales[j]*
r[0]/au);
869 static const size_t ncoef=3;
870 static G4double scales[ncoef]={-6.0, -1.2, -0.3};
871 static G4double coefs[ncoef]={0.10, 0.55, 0.35};
874 +std::pow(z2,0.6667));
875 std::vector<G4double>
r(npoints),
phi(npoints);
877 for(
size_t i=0; i<npoints; i++) {
878 G4double rr=(float)i/(
float)(npoints-1);
882 for(
size_t j=0; j<ncoef; j++)
883 sum+=coefs[j]*std::exp(scales[j]*
r[i]/au);
889 for(
size_t j=0; j<ncoef; j++)
890 phiprime0+=scales[j]*coefs[j]*std::exp(scales[j]*
r[0]/au);
904 +std::pow(z2,0.6667));
905 std::vector<G4double>
r(npoints),
phi(npoints);
907 for(
size_t i=0; i<npoints; i++) {
908 G4double rr=(float)i/(
float)(npoints-1);
914 G4double phipoly=1+y+0.3344*ysq+0.0485*y*ysq+0.002647*ysq*ysq;
915 phi[i]=phipoly*std::exp(-y);
920 G4double logphiprime0=(9.67/2.0)*(2*0.3344-1.0);
922 logphiprime0 *= (1.0/au);
968 std::vector<G4String>
970 std::vector<G4String> keys;
972 std::map<std::string, ScreeningFunc>::const_iterator sfunciter=
phiMap.begin();
973 for(; sfunciter !=
phiMap.end(); sfunciter++)
974 keys.push_back((*sfunciter).first);
987 return mc2*f/(std::sqrt(1.0+f)+1.0);
991 G4double s2th2=eratio*( (m1+mass2)*(m1+mass2)/(4.0*m1*mass2) );
993 return 2.0*std::asin(sth2);
1000 static const size_t sigLen=200;
1017 for (
G4int iEl=0; iEl<nMatElements; iEl++)
1019 G4Element* element = (*elementVector)[iEl];
1027 std::map<std::string, ScreeningFunc>::iterator sfunciter=
1028 phiMap.find(screeningKey);
1029 if(sfunciter==
phiMap.end()) {
1031 ed <<
"No such screening key <"
1032 << screeningKey <<
">";
1033 G4Exception(
"G4NativeScreenedCoulombCrossSection::LoadData",
1065 x0func->set_domain(1
e-6*
angstrom/au, 0.9999*screen->
xmax()/au);
1077 G4double eratkin=0.9999*(4*a1*a2)/((a1+a2)*(a1+a2));
1084 G4cout <<
"Native Screening: " << screeningKey <<
" "
1085 << z1 <<
" " << a1 <<
" " <<
1086 Z <<
" " << a2 <<
" " << recoilCutoff <<
G4endl;
1098 c2eps.
reset(0.0, 0.0, eps);
1112 x0=x0_solution(2*q-q*q);
1116 "failure in inverse solution to generate MFP tables");