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 #include "G4ChipsKaonMinusInelasticXS.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4DynamicParticle.hh"
00043 #include "G4ParticleDefinition.hh"
00044 #include "G4KaonMinus.hh"
00045
00046
00047 #include "G4CrossSectionFactory.hh"
00048
00049 G4_DECLARE_XS_FACTORY(G4ChipsKaonMinusInelasticXS);
00050
00051
00052
00053 G4ChipsKaonMinusInelasticXS::G4ChipsKaonMinusInelasticXS():G4VCrossSectionDataSet(Default_Name())
00054 {
00055 lastLEN=0;
00056 lastHEN=0;
00057 lastN=0;
00058 lastZ=0;
00059 lastP=0.;
00060 lastTH=0.;
00061 lastCS=0.;
00062 lastI=0;
00063 LEN = new std::vector<G4double*>;
00064 HEN = new std::vector<G4double*>;
00065 }
00066
00067 G4ChipsKaonMinusInelasticXS::~G4ChipsKaonMinusInelasticXS()
00068 {
00069 G4int lens=LEN->size();
00070 for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
00071 delete LEN;
00072
00073 G4int hens=HEN->size();
00074 for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
00075 delete HEN;
00076 }
00077
00078 G4bool G4ChipsKaonMinusInelasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00079 const G4Element*,
00080 const G4Material*)
00081 {
00082 G4ParticleDefinition* particle = Pt->GetDefinition();
00083 if (particle == G4KaonMinus::KaonMinus() ) return true;
00084 return false;
00085 }
00086
00087
00088
00089
00090 G4double G4ChipsKaonMinusInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
00091 const G4Isotope*,
00092 const G4Element*,
00093 const G4Material*)
00094 {
00095 G4double pMom=Pt->GetTotalMomentum();
00096 G4int tgN = A - tgZ;
00097
00098 return GetChipsCrossSection(pMom, tgZ, tgN, -321);
00099 }
00100
00101 G4double G4ChipsKaonMinusInelasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int)
00102 {
00103 static G4int j;
00104 static std::vector <G4int> colN;
00105 static std::vector <G4int> colZ;
00106 static std::vector <G4double> colP;
00107 static std::vector <G4double> colTH;
00108 static std::vector <G4double> colCS;
00109
00110
00111 G4bool in=false;
00112 if(tgN!=lastN || tgZ!=lastZ)
00113 {
00114 in = false;
00115 lastP = 0.;
00116 lastN = tgN;
00117 lastZ = tgZ;
00118 lastI = colN.size();
00119 j = 0;
00120 if(lastI) for(G4int i=0; i<lastI; i++)
00121 {
00122 if(colN[i]==tgN && colZ[i]==tgZ)
00123 {
00124 lastI=i;
00125 lastTH =colTH[i];
00126 if(pMom<=lastTH)
00127 {
00128 return 0.;
00129 }
00130 lastP =colP [i];
00131 lastCS =colCS[i];
00132 in = true;
00133
00134 lastCS=CalculateCrossSection(-1,j,-321,lastZ,lastN,pMom);
00135 if(lastCS<=0. && pMom>lastTH)
00136 {
00137 lastCS=0.;
00138 lastTH=pMom;
00139 }
00140 break;
00141 }
00142 j++;
00143 }
00144 if(!in)
00145 {
00147 lastCS=CalculateCrossSection(0,j,-321,lastZ,lastN,pMom);
00148
00149
00150
00151
00152
00153 lastTH = 0;
00154 colN.push_back(tgN);
00155 colZ.push_back(tgZ);
00156 colP.push_back(pMom);
00157 colTH.push_back(lastTH);
00158 colCS.push_back(lastCS);
00159
00160 return lastCS*millibarn;
00161 }
00162 else
00163 {
00164 colP[lastI]=pMom;
00165 colCS[lastI]=lastCS;
00166 }
00167 }
00168 else if(pMom<=lastTH)
00169 {
00170 return 0.;
00171 }
00172 else
00173 {
00174 lastCS=CalculateCrossSection(1,j,-321,lastZ,lastN,pMom);
00175 lastP=pMom;
00176 }
00177 return lastCS*millibarn;
00178 }
00179
00180
00181 G4double G4ChipsKaonMinusInelasticXS::CalculateCrossSection(G4int F, G4int I,
00182 G4int, G4int targZ, G4int targN, G4double Momentum)
00183 {
00184 static const G4double THmin=27.;
00185 static const G4double THmiG=THmin*.001;
00186 static const G4double dP=10.;
00187 static const G4double dPG=dP*.001;
00188 static const G4int nL=105;
00189 static const G4double Pmin=THmin+(nL-1)*dP;
00190 static const G4double Pmax=227000.;
00191 static const G4int nH=224;
00192 static const G4double milP=std::log(Pmin);
00193 static const G4double malP=std::log(Pmax);
00194 static const G4double dlP=(malP-milP)/(nH-1);
00195 static const G4double milPG=std::log(.001*Pmin);
00196 G4double sigma=0.;
00197 if(F&&I) sigma=0.;
00198
00199 if(F<=0)
00200 {
00201 if(F<0)
00202 {
00203 G4int sync=LEN->size();
00204 if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
00205 lastLEN=(*LEN)[I];
00206 lastHEN=(*HEN)[I];
00207 }
00208 else
00209 {
00210 lastLEN = new G4double[nL];
00211 lastHEN = new G4double[nH];
00212
00213 G4double P=THmiG;
00214 for(G4int k=0; k<nL; k++)
00215 {
00216 lastLEN[k] = CrossSectionLin(targZ, targN, P);
00217 P+=dPG;
00218 }
00219 G4double lP=milPG;
00220 for(G4int n=0; n<nH; n++)
00221 {
00222 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
00223 lP+=dlP;
00224 }
00225
00226
00227 G4int sync=LEN->size();
00228 if(sync!=I)
00229 {
00230 G4cerr<<"***G4ChipsKaonMinusCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
00231 <<", N="<<targN<<", F="<<F<<G4endl;
00232
00233 }
00234 LEN->push_back(lastLEN);
00235 HEN->push_back(lastHEN);
00236 }
00237 }
00238
00239 if (Momentum<lastTH) return 0.;
00240 else if (Momentum<Pmin)
00241 {
00242 sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
00243 }
00244 else if (Momentum<Pmax)
00245 {
00246 G4double lP=std::log(Momentum);
00247 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
00248 }
00249 else
00250 {
00251 G4double P=0.001*Momentum;
00252 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
00253 }
00254 if(sigma<0.) return 0.;
00255 return sigma;
00256 }
00257
00258
00259 G4double G4ChipsKaonMinusInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00260 {
00261 G4double lP=std::log(P);
00262 return CrossSectionFormula(tZ, tN, P, lP);
00263 }
00264
00265
00266 G4double G4ChipsKaonMinusInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00267 {
00268 G4double P=std::exp(lP);
00269 return CrossSectionFormula(tZ, tN, P, lP);
00270 }
00271
00272 G4double G4ChipsKaonMinusInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
00273 G4double P, G4double lP)
00274 {
00275 G4double sigma=0.;
00276 if(tZ==1 && !tN)
00277 {
00278 G4double ld=lP-3.5;
00279 G4double ld2=ld*ld;
00280 G4double p2=P*P;
00281 G4double p4=p2*p2;
00282 G4double sp=std::sqrt(P);
00283 G4double psp=P*sp;
00284 G4double lm=P-.39;
00285 G4double md=lm*lm+.000156;
00286 G4double lh=P-1.;
00287 G4double hd=lh*lh+.0156;
00288 G4double El=(.0557*ld2+2.23)/(1.-.7/sp+.075/p4);
00289 G4double To=(.3*ld2+19.5)/(1.-.21/sp+.52/p4);
00290 sigma=8.8/psp+(To-El)+.002/md+.15/hd;
00291 }
00292 else if(tZ==1 && tN==1)
00293 {
00294 G4double p2=P*P;
00295 G4double dX=lP-3.7;
00296 G4double dR=P-.94;
00297 G4double sp=std::sqrt(P);
00298 sigma=(.6*dX*dX+36.)/(1.-.11/sp+.52/p2/p2)+.7/(dR*dR+.0256)+18./P/sp;
00299 }
00300 else if(tZ<97 && tN<152)
00301 {
00302 G4double d=lP-4.2;
00303 G4double sp=std::sqrt(P);
00304 G4double p2=P*P;
00305 G4double a=tN+tZ;
00306 G4double sa=std::sqrt(a);
00307 G4double al=std::log(a);
00308 G4double a2=a*a;
00309 G4double c=52.*std::exp(al*0.6)*(1.+97./a2)/(1.+9.8/a)/(1.+47./a2);
00310 G4double gg=-.2-.003*a;
00311 G4double h=.5+.07*a;
00312 G4double v=P-1.;
00313 G4double f=.6*a*sa/(1.+.00002*a2);
00314 G4double u=.125+.127*al;
00315 sigma=(c+d*d)/(1.+gg/sp+h/p2/p2)+f/(v*v+u*u)+20.*sa/P/sp;
00316 }
00317 else
00318 {
00319 G4cerr<<"-Warning-G4ChipsKMinusNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00320 sigma=0.;
00321 }
00322 if(sigma<0.) return 0.;
00323 return sigma;
00324 }
00325
00326 G4double G4ChipsKaonMinusInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
00327 {
00328 if(DX<=0. || N<2)
00329 {
00330 G4cerr<<"***G4ChipsKaonMinusInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
00331 return Y[0];
00332 }
00333
00334 G4int N2=N-2;
00335 G4double d=(X-X0)/DX;
00336 G4int j=static_cast<int>(d);
00337 if (j<0) j=0;
00338 else if(j>N2) j=N2;
00339 d-=j;
00340 G4double yi=Y[j];
00341 G4double sigma=yi+(Y[j+1]-yi)*d;
00342
00343 return sigma;
00344 }