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 "G4ChipsPionPlusInelasticXS.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4DynamicParticle.hh"
00044 #include "G4ParticleDefinition.hh"
00045 #include "G4PionPlus.hh"
00046
00047
00048 #include "G4CrossSectionFactory.hh"
00049
00050 G4_DECLARE_XS_FACTORY(G4ChipsPionPlusInelasticXS);
00051
00052 G4ChipsPionPlusInelasticXS::G4ChipsPionPlusInelasticXS():G4VCrossSectionDataSet(Default_Name())
00053 {
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
00068 G4ChipsPionPlusInelasticXS::~G4ChipsPionPlusInelasticXS()
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 G4bool G4ChipsPionPlusInelasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00079 const G4Element*,
00080 const G4Material*)
00081 {
00082 G4ParticleDefinition* particle = Pt->GetDefinition();
00083 if (particle == G4PionPlus::PionPlus() ) return true;
00084 return false;
00085 }
00086
00087
00088
00089 G4double G4ChipsPionPlusInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
00090 const G4Isotope*,
00091 const G4Element*,
00092 const G4Material*)
00093 {
00094 G4double pMom=Pt->GetTotalMomentum();
00095 G4int tgN = A - tgZ;
00096
00097 return GetChipsCrossSection(pMom, tgZ, tgN, 211);
00098 }
00099
00100
00101 G4double G4ChipsPionPlusInelasticXS::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 G4ChipsPionPlusInelasticXS::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 G4ChipsPionPlusInelasticXS::ThresholdMomentum(G4int tZ, G4int tN)
00258 {
00259 static const G4double third=1./3.;
00260 static const G4double pM = G4PionPlus::PionPlus()->Definition()->GetPDGMass();
00261 static const G4double tpM= pM+pM;
00262 G4double tA=tZ+tN;
00263 if(tZ<.99 || tN<0.) return 0.;
00264 else if(tZ==1 && tN==0) return 300.;
00265
00266 G4double dE=tZ/(1.+std::pow(tA,third));
00267 G4double tM=931.5*tA;
00268 G4double T=dE+dE*(dE/2+pM)/tM;
00269 return std::sqrt(T*(tpM+T));
00270 }
00271
00272
00273 G4double G4ChipsPionPlusInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00274 {
00275 G4double lP=std::log(P);
00276 return CrossSectionFormula(tZ, tN, P, lP);
00277 }
00278
00279
00280 G4double G4ChipsPionPlusInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00281 {
00282 G4double P=std::exp(lP);
00283 return CrossSectionFormula(tZ, tN, P, lP);
00284 }
00285
00286 G4double G4ChipsPionPlusInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
00287 G4double P, G4double lP)
00288 {
00289 G4double sigma=0.;
00290 if(tZ==1 && !tN)
00291 {
00292 G4double ld=lP-3.5;
00293 G4double ld2=ld*ld;
00294 G4double p2=P*P;
00295 G4double p4=p2*p2;
00296 G4double sp=std::sqrt(P);
00297 G4double lm=lP-.32;
00298 G4double md=lm*lm+.04;
00299 G4double El=(.0557*ld2+2.4+6./sp)/(1.+3./p4);
00300 G4double To=(.3*ld2+22.3+5./sp)/(1.+1./p4);
00301 sigma=(To-El)+.1/md;
00302 }
00303 else if(tZ==1 && tN==1)
00304 {
00305 G4double p2=P*P;
00306 G4double d=lP-2.7;
00307 G4double f=lP+1.25;
00308 G4double gg=lP-.017;
00309 sigma=(.55*d*d+38.+23./std::sqrt(P))/(1.+.3/p2/p2)+18./(f*f+.1089)+.02/(gg*gg+.0025);
00310 }
00311 else if(tZ<97 && tN<152)
00312 {
00313 G4double d=lP-4.2;
00314 G4double p2=P*P;
00315 G4double p4=p2*p2;
00316 G4double a=tN+tZ;
00317 G4double al=std::log(a);
00318 G4double sa=std::sqrt(a);
00319 G4double ssa=std::sqrt(sa);
00320 G4double a2=a*a;
00321 G4double c=41.*std::exp(al*.68)*(1.+44./a2)/(1.+8./a)/(1.+200./a2/a2);
00322 G4double f=290.*ssa/(1.+34./a/ssa);
00323 G4double gg=-1.32-al*.043;
00324 G4double u=lP-gg;
00325 G4double h=al*(.4-.055*al);
00326 G4double r=.01+a2*5.E-8;
00327 sigma=(c+d*d)/(1.+(.2-.009*sa)/p4)+f/(u*u+h*h)/(1.+r/p2);
00328 }
00329 else
00330 {
00331 G4cerr<<"-Warning-G4ChipsPiPlusNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00332 sigma=0.;
00333 }
00334 if(sigma<0.) return 0.;
00335 return sigma;
00336 }
00337
00338 G4double G4ChipsPionPlusInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
00339 {
00340 if(DX<=0. || N<2)
00341 {
00342 G4cerr<<"***G4ChipsPionPlusInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
00343 return Y[0];
00344 }
00345
00346 G4int N2=N-2;
00347 G4double d=(X-X0)/DX;
00348 G4int j=static_cast<int>(d);
00349 if (j<0) j=0;
00350 else if(j>N2) j=N2;
00351 d-=j;
00352 G4double yi=Y[j];
00353 G4double sigma=yi+(Y[j+1]-yi)*d;
00354
00355 return sigma;
00356 }