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 "G4ChipsPionMinusInelasticXS.hh"
00042 #include "G4ChipsPionPlusInelasticXS.hh"
00043 #include "G4SystemOfUnits.hh"
00044 #include "G4DynamicParticle.hh"
00045 #include "G4ParticleDefinition.hh"
00046 #include "G4PionMinus.hh"
00047
00048
00049 #include "G4CrossSectionFactory.hh"
00050
00051 G4_DECLARE_XS_FACTORY(G4ChipsPionMinusInelasticXS);
00052
00053 G4ChipsPionMinusInelasticXS::G4ChipsPionMinusInelasticXS():G4VCrossSectionDataSet("ChipsPionMinusInelasticXS")
00054 {
00055
00056 lastLEN=0;
00057 lastHEN=0;
00058 lastN=0;
00059 lastZ=0;
00060 lastP=0.;
00061 lastTH=0.;
00062 lastCS=0.;
00063 lastI=0;
00064 LEN = new std::vector<G4double*>;
00065 HEN = new std::vector<G4double*>;
00066 }
00067
00068 G4ChipsPionMinusInelasticXS::~G4ChipsPionMinusInelasticXS()
00069 {
00070 G4int lens=LEN->size();
00071 for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
00072 delete LEN;
00073 G4int hens=HEN->size();
00074 for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
00075 delete HEN;
00076 }
00077
00078
00079 G4bool G4ChipsPionMinusInelasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00080 const G4Element*,
00081 const G4Material*)
00082 {
00083 G4ParticleDefinition* particle = Pt->GetDefinition();
00084 if (particle == G4PionMinus::PionMinus() ) return true;
00085 return false;
00086 }
00087
00088
00089
00090 G4double G4ChipsPionMinusInelasticXS::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, -211);
00099 }
00100
00101 G4double G4ChipsPionMinusInelasticXS::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,-211,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,-211,lastZ,lastN,pMom);
00148
00149
00150
00151 lastTH = 0;
00152 colN.push_back(tgN);
00153 colZ.push_back(tgZ);
00154 colP.push_back(pMom);
00155 colTH.push_back(lastTH);
00156 colCS.push_back(lastCS);
00157
00158 return lastCS*millibarn;
00159 }
00160 else
00161 {
00162 colP[lastI]=pMom;
00163 colCS[lastI]=lastCS;
00164 }
00165 }
00166 else if(pMom<=lastTH)
00167 {
00168 return 0.;
00169 }
00170 else
00171 {
00172 lastCS=CalculateCrossSection(1,j,-211,lastZ,lastN,pMom);
00173 lastP=pMom;
00174 }
00175 return lastCS*millibarn;
00176 }
00177
00178
00179 G4double G4ChipsPionMinusInelasticXS::CalculateCrossSection(G4int F, G4int I,
00180 G4int, G4int targZ, G4int targN, G4double Momentum)
00181 {
00182 static const G4double THmin=27.;
00183 static const G4double THmiG=THmin*.001;
00184 static const G4double dP=10.;
00185 static const G4double dPG=dP*.001;
00186 static const G4int nL=105;
00187 static const G4double Pmin=THmin+(nL-1)*dP;
00188 static const G4double Pmax=227000.;
00189 static const G4int nH=224;
00190 static const G4double milP=std::log(Pmin);
00191 static const G4double malP=std::log(Pmax);
00192 static const G4double dlP=(malP-milP)/(nH-1);
00193 static const G4double milPG=std::log(.001*Pmin);
00194 G4double sigma=0.;
00195 if(F&&I) sigma=0.;
00196
00197 if(F<=0)
00198 {
00199 if(F<0)
00200 {
00201 G4int sync=LEN->size();
00202 if(sync<=I) G4cerr<<"*!*G4ChipsPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
00203 lastLEN=(*LEN)[I];
00204 lastHEN=(*HEN)[I];
00205 }
00206 else
00207 {
00208 lastLEN = new G4double[nL];
00209 lastHEN = new G4double[nH];
00210
00211 G4double P=THmiG;
00212 for(G4int k=0; k<nL; k++)
00213 {
00214 lastLEN[k] = CrossSectionLin(targZ, targN, P);
00215 P+=dPG;
00216 }
00217 G4double lP=milPG;
00218 for(G4int n=0; n<nH; n++)
00219 {
00220 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
00221 lP+=dlP;
00222 }
00223
00224
00225 G4int sync=LEN->size();
00226 if(sync!=I)
00227 {
00228 G4cerr<<"***G4ChipsPiMinusNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
00229 <<", N="<<targN<<", F="<<F<<G4endl;
00230
00231 }
00232 LEN->push_back(lastLEN);
00233 HEN->push_back(lastHEN);
00234 }
00235 }
00236
00237 if (Momentum<lastTH) return 0.;
00238 else if (Momentum<Pmin)
00239 {
00240 sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
00241 }
00242 else if (Momentum<Pmax)
00243 {
00244 G4double lP=std::log(Momentum);
00245 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
00246 }
00247 else
00248 {
00249 G4double P=0.001*Momentum;
00250 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
00251 }
00252 if(sigma<0.) return 0.;
00253 return sigma;
00254 }
00255
00256
00257 G4double G4ChipsPionMinusInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00258 {
00259 G4double lP=std::log(P);
00260 return CrossSectionFormula(tZ, tN, P, lP);
00261 }
00262
00263
00264 G4double G4ChipsPionMinusInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00265 {
00266 G4double P=std::exp(lP);
00267 return CrossSectionFormula(tZ, tN, P, lP);
00268 }
00269
00270 G4double G4ChipsPionMinusInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
00271 G4double P, G4double lP)
00272 {
00273 G4double sigma=0.;
00274 if(tZ==1 && !tN)
00275 {
00276 G4double lr=lP+1.27;
00277 G4double LE=1.53/(lr*lr+.0676);
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 lm=lP+.36;
00284 G4double md=lm*lm+.04;
00285 G4double lh=lP-.017;
00286 G4double hd=lh*lh+.0025;
00287 G4double El=(.0557*ld2+2.4+7./sp)/(1.+.7/p4);
00288 G4double To=(.3*ld2+22.3+12./sp)/(1.+.4/p4);
00289 sigma=(To-El)+.4/md+.01/hd;
00290 sigma+=LE*2;
00291 }
00292 else if(tZ==1 && tN==1)
00293 {
00294 G4double p2=P*P;
00295 G4double d=lP-2.7;
00296 G4double f=lP+1.25;
00297 G4double gg=lP-.017;
00298 sigma=(.55*d*d+38.+23./std::sqrt(P))/(1.+.3/p2/p2)+18./(f*f+.1089)+.02/(gg*gg+.0025);
00299 }
00300 else if(tZ<97 && tN<152)
00301 {
00302 G4double d=lP-4.2;
00303 G4double p2=P*P;
00304 G4double p4=p2*p2;
00305 G4double a=tN+tZ;
00306 G4double al=std::log(a);
00307 G4double sa=std::sqrt(a);
00308 G4double ssa=std::sqrt(sa);
00309 G4double a2=a*a;
00310 G4double c=41.*std::exp(al*.68)*(1.+44./a2)/(1.+8./a)/(1.+200./a2/a2);
00311 G4double f=120*sa/(1.+24./a/ssa);
00312 G4double gg=-1.32-al*.043;
00313 G4double u=lP-gg;
00314 G4double h=al*(.388-.046*al);
00315 sigma=(c+d*d)/(1.+.17/p4)+f/(u*u+h*h);
00316 }
00317 else
00318 {
00319 G4cerr<<"-Warning-G4ChipsPiMinusNuclearCroSect::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 G4ChipsPionMinusInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
00327 {
00328 if(DX<=0. || N<2)
00329 {
00330 G4cerr<<"***G4ChipsPionMinusInelasticXS::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 }