00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "G4PreCompoundProton.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4Proton.hh"
00049
00050 G4PreCompoundProton::G4PreCompoundProton()
00051 : G4PreCompoundNucleon(G4Proton::Proton(), &theProtonCoulombBarrier)
00052 {
00053 ResidualA = GetRestA();
00054 ResidualZ = GetRestZ();
00055 theA = GetA();
00056 theZ = GetZ();
00057 ResidualAthrd = ResidualA13();
00058 FragmentAthrd = ResidualAthrd;
00059 FragmentA = theA + ResidualA;
00060 }
00061
00062 G4PreCompoundProton::~G4PreCompoundProton()
00063 {}
00064
00065 G4double G4PreCompoundProton::GetRj(G4int nParticles, G4int nCharged)
00066 {
00067 G4double rj = 0.0;
00068 if(nParticles > 0) {
00069 rj = static_cast<G4double>(nCharged)/static_cast<G4double>(nParticles);
00070 }
00071 return rj;
00072 }
00073
00075
00076
00077
00078
00079
00080
00081 G4double G4PreCompoundProton::CrossSection(G4double K)
00082 {
00083 ResidualA = GetRestA();
00084 ResidualZ = GetRestZ();
00085 theA = GetA();
00086 theZ = GetZ();
00087 ResidualAthrd = ResidualA13();
00088 FragmentA = theA + ResidualA;
00089 FragmentAthrd = g4pow->Z13(FragmentA);
00090
00091 if (OPTxs==0) { return GetOpt0(K); }
00092 else if( OPTxs==1) { return GetOpt1(K); }
00093 else if( OPTxs==2|| OPTxs==4) { return GetOpt2(K); }
00094 else if (OPTxs==3) { return GetOpt3(K); }
00095 else{
00096 std::ostringstream errOs;
00097 errOs << "BAD PROTON CROSS SECTION OPTION !!" <<G4endl;
00098 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00099 return 0.;
00100 }
00101 }
00102
00103 G4double G4PreCompoundProton::GetAlpha()
00104 {
00105 G4int aZ = ResidualZ;
00106 G4double C = 0.0;
00107 if (aZ >= 70)
00108 {
00109 C = 0.10;
00110 }
00111 else
00112 {
00113 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
00114 }
00115 return 1.0 + C;
00116 }
00117
00118 G4double G4PreCompoundProton::GetBeta()
00119 {
00120 return -GetCoulombBarrier();
00121 }
00122
00123
00124
00125
00126 G4double G4PreCompoundProton::GetOpt1(G4double K)
00127 {
00128 G4double Kc=K;
00129
00130
00131 if (K > 50*MeV) { Kc = 50*MeV; }
00132
00133 G4double landa, landa0, landa1, mu, mm0, mu1,nu, nu0, nu1, nu2,xs;
00134 G4double p, p0, p1, p2,Ec,delta,q,r,ji;
00135
00136 p0 = 15.72;
00137 p1 = 9.65;
00138 p2 = -449.0;
00139 landa0 = 0.00437;
00140 landa1 = -16.58;
00141 mm0 = 244.7;
00142 mu1 = 0.503;
00143 nu0 = 273.1;
00144 nu1 = -182.4;
00145 nu2 = -1.872;
00146 delta=0.;
00147
00148 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
00149 p = p0 + p1/Ec + p2/(Ec*Ec);
00150 landa = landa0*ResidualA + landa1;
00151
00152 G4double resmu1 = g4pow->powZ(ResidualA,mu1);
00153 mu = mm0*resmu1;
00154 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
00155 q = landa - nu/(Ec*Ec) - 2*p*Ec;
00156 r = mu + 2*nu/Ec + p*(Ec*Ec);
00157
00158 ji=std::max(Kc,Ec);
00159 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
00160 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
00161 if (xs <0.0) {xs=0.0;}
00162
00163 return xs;
00164 }
00165
00166
00167
00168 G4double G4PreCompoundProton::GetOpt2(G4double K)
00169 {
00170
00171 G4double eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
00172
00173
00174
00175
00176
00177 if(!useSICB && K<=theCoulombBarrier) { return 0.0; }
00178
00179 eekin=K;
00180 G4int rnneu=ResidualA-ResidualZ;
00181 ekin=eekin/1000;
00182 r0=1.36*1.e-15;
00183 fac=pi*r0*r0;
00184 b0=2.247-0.915*(1.-1./ResidualAthrd);
00185 fac1=b0*(1.-1./ResidualAthrd);
00186 fac2=1.;
00187 if(rnneu > 1.5) { fac2 = g4pow->logZ(rnneu); }
00188 xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1);
00189 xine_th=(1.-0.15*std::exp(-ekin))*xine_th/(1.00-0.0007*ResidualA);
00190 ff1=0.70-0.0020*ResidualA;
00191 ff2=1.00+1/G4double(ResidualA);
00192 ff3=0.8+18/G4double(ResidualA)-0.002*ResidualA;
00193 fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2))));
00194 xine_th=xine_th*(1.+ff3*fac);
00195 ff1=1.-1/G4double(ResidualA)-0.001*ResidualA;
00196 ff2=1.17-2.7/G4double(ResidualA)-0.0014*ResidualA;
00197 fac=-8.*ff1*(std::log10(ekin)+2.0*ff2);
00198 fac=1./(1.+std::exp(fac));
00199 xine_th=xine_th*fac;
00200 if (xine_th < 0.0){
00201 std::ostringstream errOs;
00202 G4cout<<"WARNING: negative Wellisch cross section "<<G4endl;
00203 errOs << "RESIDUAL: A=" << ResidualA << " Z=" << ResidualZ <<G4endl;
00204 errOs <<" xsec("<<ekin<<" MeV) ="<<xine_th <<G4endl;
00205 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00206 }
00207 return xine_th;
00208 }
00209
00210
00211 G4double G4PreCompoundProton::GetOpt3(const G4double K)
00212 {
00213
00214
00215
00216 G4double landa, landa0, landa1, mu, mm0, mu1,nu, nu0, nu1, nu2;
00217 G4double p, p0, p1, p2;
00218 p0 = 15.72;
00219 p1 = 9.65;
00220 p2 = -300.;
00221 landa0 = 0.00437;
00222 landa1 = -16.58;
00223 mm0 = 244.7;
00224 mu1 = 0.503;
00225 nu0 = 273.1;
00226 nu1 = -182.4;
00227 nu2 = -1.872;
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig;
00239 G4double b,ecut,cut,ecut2,geom,elab;
00240
00241 G4double flow = 1.e-18;
00242 G4double spill= 1.e+18;
00243
00244 if (ResidualA <= 60) { signor = 0.92; }
00245 else if (ResidualA < 100) { signor = 0.8 + ResidualA*0.002; }
00246
00247 ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
00248 ecsq = ec * ec;
00249 p = p0 + p1/ec + p2/ecsq;
00250 landa = landa0*ResidualA + landa1;
00251 a = g4pow->powZ(ResidualA,mu1);
00252 mu = mm0 * a;
00253 nu = a* (nu0+nu1*ec+nu2*ecsq);
00254
00255 c =std::min(3.15,ec*0.5);
00256 w = 0.7 * c / 3.15;
00257
00258 xnulam = nu / landa;
00259 if (xnulam > spill) { xnulam=0.; }
00260 if (xnulam >= flow) { etest =std::sqrt(xnulam) + 7.; }
00261
00262 a = -2.*p*ec + landa - nu/ecsq;
00263 b = p*ecsq + mu + 2.*nu/ec;
00264 ecut = 0.;
00265 cut = a*a - 4.*p*b;
00266 if (cut > 0.) { ecut = std::sqrt(cut); }
00267 ecut = (ecut-a) / (p+p);
00268 ecut2 = ecut;
00269
00270
00271
00272
00273 if (cut < 0.) { ecut2 = ecut; }
00274 elab = K * FragmentA /G4double(ResidualA);
00275 sig = 0.;
00276 if (elab <= ec) {
00277 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
00278
00279 signor2 = (ec-elab-c) / w;
00280 signor2 = 1. + std::exp(signor2);
00281 sig = sig / signor2;
00282 }
00283 else{
00284 sig = (landa*elab+mu+nu/elab) * signor;
00285 geom = 0.;
00286
00287 if (xnulam < flow || elab < etest)
00288 {
00289 if (sig <0.0) {sig=0.0;}
00290 return sig;
00291 }
00292 geom = std::sqrt(theA*K);
00293 geom = 1.23*ResidualAthrd + ra + 4.573/geom;
00294 geom = 31.416 * geom * geom;
00295 sig = std::max(geom,sig);
00296
00297 }
00298 return sig;
00299 }