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 "G4ChipsAntiBaryonInelasticXS.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4DynamicParticle.hh"
00044 #include "G4ParticleDefinition.hh"
00045 #include "G4AntiNeutron.hh"
00046 #include "G4AntiProton.hh"
00047 #include "G4AntiLambda.hh"
00048 #include "G4AntiSigmaPlus.hh"
00049 #include "G4AntiSigmaMinus.hh"
00050 #include "G4AntiSigmaZero.hh"
00051 #include "G4AntiXiMinus.hh"
00052 #include "G4AntiXiZero.hh"
00053 #include "G4AntiOmegaMinus.hh"
00054
00055
00056 #include "G4CrossSectionFactory.hh"
00057
00058 G4_DECLARE_XS_FACTORY(G4ChipsAntiBaryonInelasticXS);
00059
00060 G4ChipsAntiBaryonInelasticXS::G4ChipsAntiBaryonInelasticXS():G4VCrossSectionDataSet(Default_Name())
00061 {
00062 lastLEN=0;
00063 lastHEN=0;
00064 lastN=0;
00065 lastZ=0;
00066 lastP=0.;
00067 lastTH=0.;
00068 lastCS=0.;
00069 lastI=0;
00070 LEN = new std::vector<G4double*>;
00071 HEN = new std::vector<G4double*>;
00072 }
00073
00074 G4ChipsAntiBaryonInelasticXS::~G4ChipsAntiBaryonInelasticXS()
00075 {
00076 G4int lens=LEN->size();
00077 for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
00078 delete LEN;
00079 G4int hens=HEN->size();
00080 for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
00081 delete HEN;
00082 }
00083
00084 G4bool G4ChipsAntiBaryonInelasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00085 const G4Element*,
00086 const G4Material*)
00087 {
00088 G4ParticleDefinition* particle = Pt->GetDefinition();
00089
00090 if(particle == G4AntiNeutron::AntiNeutron())
00091 {
00092 return true;
00093 }
00094 else if(particle == G4AntiProton::AntiProton())
00095 {
00096 return true;
00097 }
00098 else if(particle == G4AntiLambda::AntiLambda())
00099 {
00100 return true;
00101 }
00102 else if(particle == G4AntiSigmaPlus::AntiSigmaPlus())
00103 {
00104 return true;
00105 }
00106 else if(particle == G4AntiSigmaMinus::AntiSigmaMinus())
00107 {
00108 return true;
00109 }
00110 else if(particle == G4AntiSigmaZero::AntiSigmaZero())
00111 {
00112 return true;
00113 }
00114 else if(particle == G4AntiXiMinus::AntiXiMinus())
00115 {
00116 return true;
00117 }
00118 else if(particle == G4AntiXiZero::AntiXiZero())
00119 {
00120 return true;
00121 }
00122 else if(particle == G4AntiOmegaMinus::AntiOmegaMinus())
00123 {
00124 return true;
00125 }
00126 return false;
00127 }
00128
00129
00130
00131 G4double G4ChipsAntiBaryonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
00132 const G4Isotope*,
00133 const G4Element*,
00134 const G4Material*)
00135 {
00136 G4double pMom=Pt->GetTotalMomentum();
00137 G4int tgN = A - tgZ;
00138 G4int pdg = Pt->GetDefinition()->GetPDGEncoding();
00139
00140 return GetChipsCrossSection(pMom, tgZ, tgN, pdg);
00141 }
00142
00143 G4double G4ChipsAntiBaryonInelasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int cPDG)
00144 {
00145 static G4int j;
00146 static std::vector <G4int> colN;
00147 static std::vector <G4int> colZ;
00148 static std::vector <G4double> colP;
00149 static std::vector <G4double> colTH;
00150 static std::vector <G4double> colCS;
00151
00152
00153 G4bool in=false;
00154 if(tgN!=lastN || tgZ!=lastZ)
00155 {
00156 in = false;
00157 lastP = 0.;
00158 lastN = tgN;
00159 lastZ = tgZ;
00160 lastI = colN.size();
00161 j = 0;
00162 if(lastI) for(G4int i=0; i<lastI; i++)
00163 {
00164 if(colN[i]==tgN && colZ[i]==tgZ)
00165 {
00166 lastI=i;
00167 lastTH =colTH[i];
00168 if(pMom<=lastTH)
00169 {
00170 return 0.;
00171 }
00172 lastP =colP [i];
00173 lastCS =colCS[i];
00174 in = true;
00175
00176 lastCS=CalculateCrossSection(-1,j,cPDG,lastZ,lastN,pMom);
00177 if(lastCS<=0. && pMom>lastTH)
00178 {
00179 lastCS=0.;
00180 lastTH=pMom;
00181 }
00182 break;
00183 }
00184 j++;
00185 }
00186 if(!in)
00187 {
00189 lastCS=CalculateCrossSection(0,j,cPDG,lastZ,lastN,pMom);
00190
00191
00192
00193 lastTH = 0;
00194 colN.push_back(tgN);
00195 colZ.push_back(tgZ);
00196 colP.push_back(pMom);
00197 colTH.push_back(lastTH);
00198 colCS.push_back(lastCS);
00199
00200 return lastCS*millibarn;
00201 }
00202 else
00203 {
00204 colP[lastI]=pMom;
00205 colCS[lastI]=lastCS;
00206 }
00207 }
00208 else if(pMom<=lastTH)
00209 {
00210 return 0.;
00211 }
00212 else
00213 {
00214 lastCS=CalculateCrossSection(1,j,cPDG,lastZ,lastN,pMom);
00215 lastP=pMom;
00216 }
00217 return lastCS*millibarn;
00218 }
00219
00220
00221 G4double G4ChipsAntiBaryonInelasticXS::CalculateCrossSection(G4int F, G4int I,
00222 G4int, G4int targZ, G4int targN, G4double Momentum)
00223 {
00224 static const G4double THmin=27.;
00225 static const G4double THmiG=THmin*.001;
00226 static const G4double dP=10.;
00227 static const G4double dPG=dP*.001;
00228 static const G4int nL=105;
00229 static const G4double Pmin=THmin+(nL-1)*dP;
00230 static const G4double Pmax=227000.;
00231 static const G4int nH=224;
00232 static const G4double milP=std::log(Pmin);
00233 static const G4double malP=std::log(Pmax);
00234 static const G4double dlP=(malP-milP)/(nH-1);
00235 static const G4double milPG=std::log(.001*Pmin);
00236 G4double sigma=0.;
00237 if(F&&I) sigma=0.;
00238
00239 if(F<=0)
00240 {
00241 if(F<0)
00242 {
00243 G4int sync=LEN->size();
00244 if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
00245 lastLEN=(*LEN)[I];
00246 lastHEN=(*HEN)[I];
00247 }
00248 else
00249 {
00250 lastLEN = new G4double[nL];
00251 lastHEN = new G4double[nH];
00252
00253 G4double P=THmiG;
00254 for(G4int k=0; k<nL; k++)
00255 {
00256 lastLEN[k] = CrossSectionLin(targZ, targN, P);
00257 P+=dPG;
00258 }
00259 G4double lP=milPG;
00260 for(G4int n=0; n<nH; n++)
00261 {
00262 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
00263 lP+=dlP;
00264 }
00265
00266
00267 G4int sync=LEN->size();
00268 if(sync!=I)
00269 {
00270 G4cerr<<"***G4QPiMinusNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
00271 <<", N="<<targN<<", F="<<F<<G4endl;
00272
00273 }
00274 LEN->push_back(lastLEN);
00275 HEN->push_back(lastHEN);
00276 }
00277 }
00278
00279 if (Momentum<lastTH) return 0.;
00280 else if (Momentum<Pmin)
00281 {
00282 sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
00283 }
00284 else if (Momentum<Pmax)
00285 {
00286 G4double lP=std::log(Momentum);
00287 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
00288 }
00289 else
00290 {
00291 G4double P=0.001*Momentum;
00292 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
00293 }
00294 if(sigma<0.) return 0.;
00295 return sigma;
00296 }
00297
00298
00299 G4double G4ChipsAntiBaryonInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00300 {
00301 G4double lP=std::log(P);
00302 return CrossSectionFormula(tZ, tN, P, lP);
00303 }
00304
00305
00306 G4double G4ChipsAntiBaryonInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00307 {
00308 G4double P=std::exp(lP);
00309 return CrossSectionFormula(tZ, tN, P, lP);
00310 }
00311
00312 G4double G4ChipsAntiBaryonInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
00313 G4double P, G4double lP)
00314 {
00315 G4double sigma=0.;
00316 if(tZ==1 && !tN)
00317 {
00318 G4double ld=lP-3.5;
00319 G4double ld2=ld*ld;
00320 G4double ye=std::exp(lP*1.25);
00321 G4double yt=std::exp(lP*0.35);
00322 G4double El=80./(ye+1.);
00323 G4double To=(80./yt+.3)/yt;
00324 sigma=(To-El)+.2443*ld2+31.48;
00325 }
00326 else if(tZ==1 && tN==1)
00327 {
00328 G4double r=lP-3.7;
00329 sigma=0.6*r*r+67.+90.*std::exp(-lP*.666);
00330 }
00331 else if(tZ<97 && tN<152)
00332 {
00333 G4double d=lP-4.2;
00334 G4double sp=std::sqrt(P);
00335 G4double a=tN+tZ;
00336 G4double sa=std::sqrt(a);
00337 G4double a2=a*a;
00338 G4double a3=a2*a;
00339 G4double a2s=a2*sa;
00340 G4double c=(170.+3600./a2s)/(1.+65./a2s)+40.*std::pow(a,0.712)/(1.+12.2/a)/(1.+34./a2);
00341 G4double r=(170.+0.01*a3)/(1.+a3/28000.);
00342 sigma=c+d*d+r/sp;
00343 }
00344 else
00345 {
00346 G4cerr<<"-Warning-G4QAntiBarNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00347 sigma=0.;
00348 }
00349 if(sigma<0.) return 0.;
00350 return sigma;
00351 }
00352
00353 G4double G4ChipsAntiBaryonInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
00354 {
00355 if(DX<=0. || N<2)
00356 {
00357 G4cerr<<"***G4ChipsAntiBaryonInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
00358 return Y[0];
00359 }
00360
00361 G4int N2=N-2;
00362 G4double d=(X-X0)/DX;
00363 G4int j=static_cast<int>(d);
00364 if (j<0) j=0;
00365 else if(j>N2) j=N2;
00366 d-=j;
00367 G4double yi=Y[j];
00368 G4double sigma=yi+(Y[j+1]-yi)*d;
00369
00370 return sigma;
00371 }