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 #include "G4ChipsKaonPlusInelasticXS.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4DynamicParticle.hh"
00044 #include "G4ParticleDefinition.hh"
00045 #include "G4KaonPlus.hh"
00046 #include "G4Proton.hh"
00047 #include "G4PionPlus.hh"
00048
00049
00050 #include "G4CrossSectionFactory.hh"
00051
00052 G4_DECLARE_XS_FACTORY(G4ChipsKaonPlusInelasticXS);
00053
00054
00055 G4ChipsKaonPlusInelasticXS::G4ChipsKaonPlusInelasticXS():G4VCrossSectionDataSet(Default_Name())
00056 {
00057
00058 lastLEN=0;
00059 lastHEN=0;
00060 lastN=0;
00061 lastZ=0;
00062 lastP=0.;
00063 lastTH=0.;
00064 lastCS=0.;
00065 lastI=0;
00066 LEN = new std::vector<G4double*>;
00067 HEN = new std::vector<G4double*>;
00068 }
00069
00070 G4ChipsKaonPlusInelasticXS::~G4ChipsKaonPlusInelasticXS()
00071 {
00072 G4int lens=LEN->size();
00073 for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
00074 delete LEN;
00075
00076 G4int hens=HEN->size();
00077 for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
00078 delete HEN;
00079 }
00080
00081
00082 G4bool G4ChipsKaonPlusInelasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00083 const G4Element*,
00084 const G4Material*)
00085 {
00086 G4ParticleDefinition* particle = Pt->GetDefinition();
00087 if (particle == G4KaonPlus::KaonPlus() ) return true;
00088 return false;
00089 }
00090
00091
00092
00093
00094 G4double G4ChipsKaonPlusInelasticXS::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, 321);
00103 }
00104
00105 G4double G4ChipsKaonPlusInelasticXS::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
00125 if(lastI) for(G4int i=0; i<lastI; i++)
00126 {
00127 if(colN[i]==tgN && colZ[i]==tgZ)
00128 {
00129 lastI=i;
00130 lastTH =colTH[i];
00131
00132 if(pMom<=lastTH)
00133 {
00134 return 0.;
00135 }
00136 lastP =colP [i];
00137 lastCS =colCS[i];
00138 in = true;
00139
00140 lastCS=CalculateCrossSection(-1,j,321,lastZ,lastN,pMom);
00141
00142 if(lastCS<=0. && pMom>lastTH)
00143 {
00144 lastCS=0.;
00145 lastTH=pMom;
00146 }
00147 break;
00148 }
00149 j++;
00150 }
00151 if(!in)
00152 {
00154 lastCS=CalculateCrossSection(0,j,321,lastZ,lastN,pMom);
00155
00156
00157
00158
00159 lastTH = 0;
00160 colN.push_back(tgN);
00161 colZ.push_back(tgZ);
00162 colP.push_back(pMom);
00163 colTH.push_back(lastTH);
00164 colCS.push_back(lastCS);
00165
00166 return lastCS*millibarn;
00167 }
00168 else
00169 {
00170 colP[lastI]=pMom;
00171 colCS[lastI]=lastCS;
00172 }
00173 }
00174 else if(pMom<=lastTH)
00175 {
00176 return 0.;
00177 }
00178 else
00179 {
00180 lastCS=CalculateCrossSection(1,j,321,lastZ,lastN,pMom);
00181 lastP=pMom;
00182 }
00183 return lastCS*millibarn;
00184 }
00185
00186
00187 G4double G4ChipsKaonPlusInelasticXS::CalculateCrossSection(G4int F, G4int I,
00188 G4int, G4int targZ, G4int targN, G4double Momentum)
00189 {
00190 static const G4double THmin=27.;
00191 static const G4double THmiG=THmin*.001;
00192 static const G4double dP=10.;
00193 static const G4double dPG=dP*.001;
00194 static const G4int nL=105;
00195 static const G4double Pmin=THmin+(nL-1)*dP;
00196 static const G4double Pmax=227000.;
00197 static const G4int nH=224;
00198 static const G4double milP=std::log(Pmin);
00199 static const G4double malP=std::log(Pmax);
00200 static const G4double dlP=(malP-milP)/(nH-1);
00201 static const G4double milPG=std::log(.001*Pmin);
00202
00203 G4double sigma=0.;
00204 if(F&&I) sigma=0.;
00205 G4double A=targN+targZ;
00206
00207 if(F<=0)
00208 {
00209 if(F<0)
00210 {
00211 G4int sync=LEN->size();
00212 if(sync<=I) G4cerr<<"*!*G4ChipsKPlusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
00213 lastLEN=(*LEN)[I];
00214 lastHEN=(*HEN)[I];
00215 }
00216 else
00217 {
00218 lastLEN = new G4double[nL];
00219 lastHEN = new G4double[nH];
00220
00221 G4double P=THmiG;
00222 for(G4int k=0; k<nL; k++)
00223 {
00224 lastLEN[k] = CrossSectionLin(targZ, targN, P);
00225 P+=dPG;
00226 }
00227 G4double lP=milPG;
00228 for(G4int n=0; n<nH; n++)
00229 {
00230 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
00231 lP+=dlP;
00232 }
00233
00234
00235 G4int sync=LEN->size();
00236 if(sync!=I)
00237 {
00238 G4cerr<<"***G4ChipsKPlusNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
00239 <<", N="<<targN<<", F="<<F<<G4endl;
00240
00241 }
00242 LEN->push_back(lastLEN);
00243 HEN->push_back(lastHEN);
00244 }
00245 }
00246
00247
00248 if (Momentum<lastTH) return 0.;
00249 else if (Momentum<Pmin)
00250 {
00251 if(A<=1. && Momentum < 600.) sigma=0.;
00252 else sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
00253 }
00254 else if (Momentum<Pmax)
00255 {
00256 G4double lP=std::log(Momentum);
00257 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
00258 }
00259 else
00260 {
00261 G4double P=0.001*Momentum;
00262 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
00263 }
00264 if(sigma<0.) return 0.;
00265 return sigma;
00266 }
00267
00268
00269 G4double G4ChipsKaonPlusInelasticXS::ThresholdMomentum(G4int tZ, G4int tN)
00270 {
00271 static const G4double third=1./3.;
00272 static const G4double prM = G4Proton::Proton()->GetPDGMass();
00273 static const G4double piM = G4PionPlus::PionPlus()->GetPDGMass()+.1;
00274 static const G4double pM = G4KaonPlus::KaonPlus()->GetPDGMass();
00275 static const G4double tpM= pM+pM;
00276 G4double tA=tZ+tN;
00277 if(tZ<.99 || tN<0.) return 0.;
00278 G4double tM=931.5*tA;
00279 G4double dE=piM;
00280 if(tZ==1 && tN==0) tM=prM;
00281 else dE=tZ/(1.+std::pow(tA,third));
00282
00283 G4double T=dE+dE*(dE/2+pM)/tM;
00284 return std::sqrt(T*(tpM+T));
00285 }
00286
00287
00288 G4double G4ChipsKaonPlusInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00289 {
00290 G4double lP=std::log(P);
00291 return CrossSectionFormula(tZ, tN, P, lP);
00292 }
00293
00294
00295 G4double G4ChipsKaonPlusInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00296 {
00297 G4double P=std::exp(lP);
00298 return CrossSectionFormula(tZ, tN, P, lP);
00299 }
00300
00301 G4double G4ChipsKaonPlusInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
00302 G4double P, G4double lP)
00303 {
00304 G4double sigma=0.;
00305 if(tZ==1 && !tN)
00306 {
00307 G4double ld=lP-3.5;
00308 G4double ld2=ld*ld;
00309 G4double sp=std::sqrt(P);
00310 G4double p2=P*P;
00311 G4double p4=p2*p2;
00312 G4double lm=P-1.;
00313 G4double md=lm*lm+.372;
00314 G4double El=(.0557*ld2+2.23)/(1.-.7/sp+.1/p4);
00315 G4double To=(.3*ld2+19.5)/(1.+.46/sp+1.6/p4);
00316 sigma=(To-El)+.6/md;
00317 }
00318 else if(tZ<97 && tN<152)
00319 {
00320 G4double p2=P*P;
00321 G4double p4=p2*p2;
00322 G4double a=tN+tZ;
00323 G4double al=std::log(a);
00324 G4double sa=std::sqrt(a);
00325 G4double asa=a*sa;
00326 G4double a2=a*a;
00327 G4double a3=a2*a;
00328 G4double a4=a2*a2;
00329 G4double a8=a4*a4;
00330 G4double a12=a8*a4;
00331 G4double f=.6;
00332 G4double r=.5;
00333 G4double gg=3.7;
00334 G4double c=36.;
00335 G4double ss=3.5;
00336 G4double t=3.;
00337 G4double u=.44;
00338 G4double v=5.E-9;
00339 if(tZ>1 && tN>1)
00340 {
00341 f=1.;
00342 r=1./(1.+.007*a2);
00343 gg=4.2;
00344 c=52.*std::exp(al*.6)*(1.+95./a2)/(1.+9./a)/(1.+46./a2);
00345 ss=(40.+.14*a)/(1.+12./a);
00346 G4double y=std::exp(al*1.7);
00347 t=.185*y/(1.+.00012*y);
00348 u=(1.+80./asa)/(1.+200./asa);
00349 v=(1.+3.E-6*a4*(1.+6.E-7*a3+4.E10/a12))/a3/20000.;
00350 }
00351 G4double d=lP-gg;
00352 G4double w=P-1.;
00353 G4double rD=ss/(w*w+.36);
00354 G4double h=P-.44;
00355 G4double rR=t/(h*h+u*u);
00356 sigma=(f*d*d+c)/(1.+r/std::sqrt(P)+1./p4)+(rD+rR)/(1+v/p4/p4);
00357 }
00358 else
00359 {
00360 G4cerr<<"-Warning-G4ChipsKaonPlusNuclearCroSect::CSForm:Bad A, Z="<<tZ<<", N="<<tN<<G4endl;
00361 sigma=0.;
00362 }
00363 if(sigma<0.) return 0.;
00364 return sigma;
00365 }
00366
00367 G4double G4ChipsKaonPlusInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
00368 {
00369 if(DX<=0. || N<2)
00370 {
00371 G4cerr<<"***G4ChipsKaonPlusInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
00372 return Y[0];
00373 }
00374
00375 G4int N2=N-2;
00376 G4double d=(X-X0)/DX;
00377 G4int j=static_cast<int>(d);
00378 if (j<0) j=0;
00379 else if(j>N2) j=N2;
00380 d-=j;
00381 G4double yi=Y[j];
00382 G4double sigma=yi+(Y[j+1]-yi)*d;
00383
00384 return sigma;
00385 }