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 #include "G4ChipsProtonInelasticXS.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4DynamicParticle.hh"
00046 #include "G4ParticleDefinition.hh"
00047 #include "G4Proton.hh"
00048
00049
00050 #include "G4CrossSectionFactory.hh"
00051
00052 G4_DECLARE_XS_FACTORY(G4ChipsProtonInelasticXS);
00053
00054 G4ChipsProtonInelasticXS::G4ChipsProtonInelasticXS():G4VCrossSectionDataSet(Default_Name())
00055 {
00056
00057 lastLEN=0;
00058 lastHEN=0;
00059 lastN=0;
00060 lastZ=0;
00061 lastP=0.;
00062 lastTH=0.;
00063 lastCS=0.;
00064 lastI=0;
00065
00066 LEN = new std::vector<G4double*>;
00067 HEN = new std::vector<G4double*>;
00068 }
00069
00070 G4ChipsProtonInelasticXS::~G4ChipsProtonInelasticXS()
00071 {
00072
00073
00074
00075
00076
00077
00078
00079
00080 }
00081
00082 G4bool G4ChipsProtonInelasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00083 const G4Element*,
00084 const G4Material*)
00085 {
00086 G4ParticleDefinition* particle = Pt->GetDefinition();
00087 if (particle == G4Proton::Proton() ) return true;
00088 return false;
00089 }
00090
00091
00092
00093
00094 G4double G4ChipsProtonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
00095 const G4Isotope*,
00096 const G4Element*,
00097 const G4Material*)
00098 {
00099 G4double pMom=Pt->GetTotalMomentum();
00100 G4int tgN = A - tgZ;
00101
00102 return GetChipsCrossSection(pMom, tgZ, tgN, 2212);
00103 }
00104
00105 G4double G4ChipsProtonInelasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int)
00106 {
00107 static G4int j;
00108 static std::vector <G4int> colN;
00109 static std::vector <G4int> colZ;
00110 static std::vector <G4double> colP;
00111 static std::vector <G4double> colTH;
00112 static std::vector <G4double> colCS;
00113
00114
00115 G4bool in=false;
00116 if(tgN!=lastN || tgZ!=lastZ)
00117 {
00118 in = false;
00119 lastP = 0.;
00120 lastN = tgN;
00121 lastZ = tgZ;
00122 lastI = colN.size();
00123 j = 0;
00124 if(lastI) for(G4int i=0; i<lastI; i++)
00125 {
00126 if(colN[i]==tgN && colZ[i]==tgZ)
00127 {
00128 lastI=i;
00129 lastTH =colTH[i];
00130 if(pMom<=lastTH)
00131 {
00132 return 0.;
00133 }
00134 lastP =colP [i];
00135 lastCS =colCS[i];
00136 in = true;
00137
00138 lastCS=CalculateCrossSection(-1,j,2212,lastZ,lastN,pMom);
00139 if(lastCS<=0. && pMom>lastTH)
00140 {
00141 lastCS=0.;
00142 lastTH=pMom;
00143 }
00144 break;
00145 }
00146 j++;
00147 }
00148 if(!in)
00149 {
00151 lastCS=CalculateCrossSection(0,j,2212,lastZ,lastN,pMom);
00152
00153
00154
00155 lastTH = 0;
00156 colN.push_back(tgN);
00157 colZ.push_back(tgZ);
00158 colP.push_back(pMom);
00159 colTH.push_back(lastTH);
00160 colCS.push_back(lastCS);
00161
00162 return lastCS*millibarn;
00163 }
00164 else
00165 {
00166 colP[lastI]=pMom;
00167 colCS[lastI]=lastCS;
00168 }
00169 }
00170 else if(pMom<=lastTH)
00171 {
00172 return 0.;
00173 }
00174 else
00175 {
00176 lastCS=CalculateCrossSection(1,j,2212,lastZ,lastN,pMom);
00177 lastP=pMom;
00178 }
00179 return lastCS*millibarn;
00180 }
00181
00182
00183 G4double G4ChipsProtonInelasticXS::CalculateCrossSection(G4int F, G4int I,
00184 G4int, G4int targZ, G4int targN, G4double Momentum)
00185 {
00186 static const G4double THmin=27.;
00187 static const G4double THmiG=THmin*.001;
00188 static const G4double dP=10.;
00189 static const G4double dPG=dP*.001;
00190 static const G4int nL=105;
00191 static const G4double Pmin=THmin+(nL-1)*dP;
00192 static const G4double Pmax=227000.;
00193 static const G4int nH=224;
00194 static const G4double milP=std::log(Pmin);
00195 static const G4double malP=std::log(Pmax);
00196 static const G4double dlP=(malP-milP)/(nH-1);
00197 static const G4double milPG=std::log(.001*Pmin);
00198 G4double sigma=0.;
00199 if(F&&I) sigma=0.;
00200
00201 if(F<=0)
00202 {
00203 if(F<0)
00204 {
00205 G4int sync=LEN->size();
00206 if(sync<=I) G4cout<<"*!*G4QProtonNuclCS::CalcCrossSect:Sync="<<sync<<"<="<<I<<G4endl;
00207 lastLEN=(*LEN)[I];
00208 lastHEN=(*HEN)[I];
00209 }
00210 else
00211 {
00212 lastLEN = new G4double[nL];
00213 lastHEN = new G4double[nH];
00214
00215 G4double P=THmiG;
00216 for(G4int k=0; k<nL; k++)
00217 {
00218 lastLEN[k] = CrossSectionLin(targZ, targN, P);
00219 P+=dPG;
00220 }
00221 G4double lP=milPG;
00222 for(G4int n=0; n<nH; n++)
00223 {
00224 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
00225 lP+=dlP;
00226 }
00227
00228
00229 G4int sync=LEN->size();
00230 if(sync!=I)
00231 {
00232 G4cout<<"***G4ChipsProtonNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
00233 <<", N="<<targN<<", F="<<F<<G4endl;
00234
00235 }
00236 LEN->push_back(lastLEN);
00237 HEN->push_back(lastHEN);
00238 }
00239 }
00240
00241 if (Momentum<lastTH) return 0.;
00242 else if (Momentum<Pmin)
00243 {
00244 sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
00245 }
00246 else if (Momentum<Pmax)
00247 {
00248 G4double lP=std::log(Momentum);
00249 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
00250 }
00251 else
00252 {
00253 G4double P=0.001*Momentum;
00254 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
00255 }
00256 if(sigma<0.) return 0.;
00257 return sigma;
00258 }
00259
00260
00261 G4double G4ChipsProtonInelasticXS::ThresholdMomentum(G4int tZ, G4int tN)
00262 {
00263 static const G4double third=1./3.;
00264 static const G4double pM = G4Proton::Proton()->Definition()->GetPDGMass();
00265 static const G4double tpM= pM+pM;
00266
00267 G4double tA=tZ+tN;
00268 if(tZ<.99 || tN<0.) return 0.;
00269 else if(tZ==1 && tN==0) return 800.;
00270
00271 G4double dE=tZ/(1.+std::pow(tA,third));
00272 G4double tM=931.5*tA;
00273 G4double T=dE+dE*(dE/2+pM)/tM;
00274 return std::sqrt(T*(tpM+T));
00275 }
00276
00277
00278 G4double G4ChipsProtonInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00279 {
00280 G4double sigma=0.;
00281 if(P<ThresholdMomentum(tZ,tN)*.001) return sigma;
00282 G4double lP=std::log(P);
00283 if(tZ==1&&!tN){if(P>.35) sigma=CrossSectionFormula(tZ,tN,P,lP);}
00284 else if(tZ<97 && tN<152)
00285 {
00286 G4double pex=0.;
00287 G4double pos=0.;
00288 G4double wid=1.;
00289 if(tZ==13 && tN==14)
00290 {
00291 pex=230.;
00292 pos=.13;
00293 wid=8.e-5;
00294 }
00295 else if(tZ<7)
00296 {
00297 if(tZ==6 && tN==6)
00298 {
00299 pex=320.;
00300 pos=.14;
00301 wid=7.e-6;
00302 }
00303 else if(tZ==5 && tN==6)
00304 {
00305 pex=270.;
00306 pos=.17;
00307 wid=.002;
00308 }
00309 else if(tZ==4 && tN==5)
00310 {
00311 pex=600.;
00312 pos=.132;
00313 wid=.005;
00314 }
00315 else if(tZ==3 && tN==4)
00316 {
00317 pex=280.;
00318 pos=.19;
00319 wid=.0025;
00320 }
00321 else if(tZ==3 && tN==3)
00322 {
00323 pex=370.;
00324 pos=.171;
00325 wid=.006;
00326 }
00327 else if(tZ==2 && tN==1)
00328 {
00329 pex=30.;
00330 pos=.22;
00331 wid=.0005;
00332 }
00333 }
00334 sigma=CrossSectionFormula(tZ,tN,P,lP);
00335 if(pex>0.)
00336 {
00337 G4double dp=P-pos;
00338 sigma+=pex*std::exp(-dp*dp/wid);
00339 }
00340 }
00341 else
00342 {
00343 G4cerr<<"-Warning-G4ChipsProtonNuclearXS::CSLin:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00344 sigma=0.;
00345 }
00346 if(sigma<0.) return 0.;
00347 return sigma;
00348 }
00349
00350
00351 G4double G4ChipsProtonInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00352 {
00353 G4double P=std::exp(lP);
00354 return CrossSectionFormula(tZ, tN, P, lP);
00355 }
00356
00357 G4double G4ChipsProtonInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
00358 G4double P, G4double lP)
00359 {
00360 G4double sigma=0.;
00361 if(tZ==1 && !tN)
00362 {
00363 G4double p2=P*P;
00364 G4double lp=lP-3.5;
00365 G4double lp2=lp*lp;
00366 G4double rp2=1./p2;
00367 G4double El=(.0557*lp2+6.72+30./P)/(1.+.49*rp2/P);
00368 G4double To=(.3*lp2+38.2)/(1.+.54*rp2*rp2);
00369 sigma=To-El;
00370 }
00371 else if(tZ<97 && tN<152)
00372 {
00373
00374 G4double d=lP-4.2;
00375 G4double p2=P*P;
00376 G4double p4=p2*p2;
00377 G4double a=tN+tZ;
00378 G4double al=std::log(a);
00379 G4double sa=std::sqrt(a);
00380 G4double a2=a*a;
00381 G4double a2s=a2*sa;
00382 G4double a4=a2*a2;
00383 G4double a8=a4*a4;
00384 G4double a12=a8*a4;
00385 G4double a16=a8*a8;
00386 G4double c=(170.+3600./a2s)/(1.+65./a2s);
00387 G4double dl=al-3.;
00388 G4double dl2=dl*dl;
00389 G4double r=.21+.62*dl2/(1.+.5*dl2);
00390 G4double gg=40.*std::exp(al*0.712)/(1.+12.2/a)/(1.+34./a2);
00391 G4double e=318.+a4/(1.+.0015*a4/std::exp(al*0.09))/(1.+4.e-28*a12)+
00392 8.e-18/(1./a16+1.3e-20)/(1.+1.e-21*a12);
00393 G4double ss=3.57+.009*a2/(1.+.0001*a2*a);
00394 G4double h=(.01/a4+2.5e-6/a)*(1.+6.e-6*a2*a)/(1.+6.e7/a12/a2);
00395 sigma=(c+d*d)/(1.+r/p4)+(gg+e*std::exp(-ss*P))/(1.+h/p4/p4);
00396 }
00397 else
00398 {
00399 G4cerr<<"-Warning-G4QProtonNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00400 sigma=0.;
00401 }
00402 if(sigma<0.) return 0.;
00403 return sigma;
00404 }
00405
00406 G4double G4ChipsProtonInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
00407 {
00408 if(DX<=0. || N<2)
00409 {
00410 G4cerr<<"***G4ChipsProtonInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
00411 return Y[0];
00412 }
00413
00414 G4int N2=N-2;
00415 G4double d=(X-X0)/DX;
00416 G4int j=static_cast<int>(d);
00417 if (j<0) j=0;
00418 else if(j>N2) j=N2;
00419 d-=j;
00420 G4double yi=Y[j];
00421 G4double sigma=yi+(Y[j+1]-yi)*d;
00422
00423 return sigma;
00424 }