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 #include "G4ChipsHyperonInelasticXS.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4DynamicParticle.hh"
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4Lambda.hh"
00041 #include "G4SigmaPlus.hh"
00042 #include "G4SigmaMinus.hh"
00043 #include "G4SigmaZero.hh"
00044 #include "G4XiMinus.hh"
00045 #include "G4XiZero.hh"
00046 #include "G4OmegaMinus.hh"
00047
00048
00049 #include "G4CrossSectionFactory.hh"
00050
00051 G4_DECLARE_XS_FACTORY(G4ChipsHyperonInelasticXS);
00052
00053 G4ChipsHyperonInelasticXS::G4ChipsHyperonInelasticXS():G4VCrossSectionDataSet(Default_Name())
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 G4ChipsHyperonInelasticXS::~G4ChipsHyperonInelasticXS()
00069 {
00070 G4int lens=LEN->size();
00071 for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
00072 delete LEN;
00073
00074 G4int hens=HEN->size();
00075 for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
00076 delete HEN;
00077 }
00078
00079 G4bool G4ChipsHyperonInelasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00080 const G4Element*,
00081 const G4Material*)
00082 {
00083 G4ParticleDefinition* particle = Pt->GetDefinition();
00084 if (particle == G4Lambda::Lambda())
00085 {
00086 return true;
00087 }
00088 else if(particle == G4SigmaPlus::SigmaPlus())
00089 {
00090 return true;
00091 }
00092 else if(particle == G4SigmaMinus::SigmaMinus())
00093 {
00094 return true;
00095 }
00096 else if(particle == G4SigmaZero::SigmaZero())
00097 {
00098 return true;
00099 }
00100 else if(particle == G4XiMinus::XiMinus())
00101 {
00102 return true;
00103 }
00104 else if(particle == G4XiZero::XiZero())
00105 {
00106 return true;
00107 }
00108 else if(particle == G4OmegaMinus::OmegaMinus())
00109 {
00110 return true;
00111 }
00112 return false;
00113 }
00114
00115
00116
00117 G4double G4ChipsHyperonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
00118 const G4Isotope*,
00119 const G4Element*,
00120 const G4Material*)
00121 {
00122 G4double pMom=Pt->GetTotalMomentum();
00123 G4int tgN = A - tgZ;
00124 G4int pdg = Pt->GetDefinition()->GetPDGEncoding();
00125
00126 return GetChipsCrossSection(pMom, tgZ, tgN, pdg);
00127 }
00128
00129 G4double G4ChipsHyperonInelasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int PDG)
00130 {
00131 static G4int j;
00132 static std::vector <G4int> colN;
00133 static std::vector <G4int> colZ;
00134 static std::vector <G4double> colP;
00135 static std::vector <G4double> colTH;
00136 static std::vector <G4double> colCS;
00137
00138
00139 G4bool in=false;
00140 if(tgN!=lastN || tgZ!=lastZ)
00141 {
00142 in = false;
00143 lastP = 0.;
00144 lastN = tgN;
00145 lastZ = tgZ;
00146 lastI = colN.size();
00147 j = 0;
00148
00149 if(lastI) for(G4int i=0; i<lastI; i++)
00150 {
00151 if(colN[i]==tgN && colZ[i]==tgZ)
00152 {
00153 lastI=i;
00154 lastTH =colTH[i];
00155
00156 if(pMom<=lastTH)
00157 {
00158 return 0.;
00159 }
00160 lastP =colP [i];
00161 lastCS =colCS[i];
00162 in = true;
00163
00164 lastCS=CalculateCrossSection(-1,j,PDG,lastZ,lastN,pMom);
00165
00166 if(lastCS<=0. && pMom>lastTH)
00167 {
00168 lastCS=0.;
00169 lastTH=pMom;
00170 }
00171 break;
00172 }
00173 j++;
00174 }
00175 if(!in)
00176 {
00178 lastCS=CalculateCrossSection(0,j,PDG,lastZ,lastN,pMom);
00179
00180
00181
00182 lastTH = 0;
00183 colN.push_back(tgN);
00184 colZ.push_back(tgZ);
00185 colP.push_back(pMom);
00186 colTH.push_back(lastTH);
00187 colCS.push_back(lastCS);
00188
00189 return lastCS*millibarn;
00190 }
00191 else
00192 {
00193 colP[lastI]=pMom;
00194 colCS[lastI]=lastCS;
00195 }
00196 }
00197 else if(pMom<=lastTH)
00198 {
00199 return 0.;
00200 }
00201 else
00202 {
00203 lastCS=CalculateCrossSection(1,j,PDG,lastZ,lastN,pMom);
00204 lastP=pMom;
00205 }
00206 return lastCS*millibarn;
00207 }
00208
00209
00210 G4double G4ChipsHyperonInelasticXS::CalculateCrossSection(G4int F, G4int I,
00211 G4int, G4int targZ, G4int targN, G4double Momentum)
00212 {
00213 static const G4double THmin=27.;
00214 static const G4double THmiG=THmin*.001;
00215 static const G4double dP=10.;
00216 static const G4double dPG=dP*.001;
00217 static const G4int nL=105;
00218 static const G4double Pmin=THmin+(nL-1)*dP;
00219 static const G4double Pmax=227000.;
00220 static const G4int nH=224;
00221 static const G4double milP=std::log(Pmin);
00222 static const G4double malP=std::log(Pmax);
00223 static const G4double dlP=(malP-milP)/(nH-1);
00224 static const G4double milPG=std::log(.001*Pmin);
00225 G4double sigma=0.;
00226 if(F&&I) sigma=0.;
00227
00228 if(F<=0)
00229 {
00230 if(F<0)
00231 {
00232 G4int sync=LEN->size();
00233 if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
00234 lastLEN=(*LEN)[I];
00235 lastHEN=(*HEN)[I];
00236 }
00237 else
00238 {
00239 lastLEN = new G4double[nL];
00240 lastHEN = new G4double[nH];
00241
00242 G4double P=THmiG;
00243 for(G4int k=0; k<nL; k++)
00244 {
00245 lastLEN[k] = CrossSectionLin(targZ, targN, P);
00246 P+=dPG;
00247 }
00248 G4double lP=milPG;
00249 for(G4int n=0; n<nH; n++)
00250 {
00251 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
00252 lP+=dlP;
00253 }
00254
00255
00256 G4int sync=LEN->size();
00257 if(sync!=I)
00258 {
00259 G4cerr<<"***G4QHyperNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
00260 <<", N="<<targN<<", F="<<F<<G4endl;
00261
00262 }
00263 LEN->push_back(lastLEN);
00264 HEN->push_back(lastHEN);
00265 }
00266 }
00267
00268 if (Momentum<lastTH) return 0.;
00269 else if (Momentum<Pmin)
00270 {
00271 sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
00272 }
00273 else if (Momentum<Pmax)
00274 {
00275 G4double lP=std::log(Momentum);
00276 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
00277 }
00278 else
00279 {
00280 G4double P=0.001*Momentum;
00281 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
00282 }
00283 if(sigma<0.) return 0.;
00284 return sigma;
00285 }
00286
00287
00288 G4double G4ChipsHyperonInelasticXS::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 G4ChipsHyperonInelasticXS::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 G4ChipsHyperonInelasticXS::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 p2=P*P;
00310 G4double p4=p2*p2;
00311 G4double sp=std::sqrt(P);
00312 G4double El=(.0557*ld2+6.72+99./p2)/(1.+2./sp+2./p4);
00313 G4double To=(.3*ld2+38.2+900./sp)/(1.+27./sp+3./p4);
00314 sigma=To-El;
00315 }
00316 else if(tZ<97 && tN<152)
00317 {
00318 G4double d=lP-4.2;
00319 G4double p2=P*P;
00320 G4double p4=p2*p2;
00321 G4double sp=std::sqrt(P);
00322 G4double ssp=std::sqrt(sp);
00323 G4double a=tN+tZ;
00324 G4double al=std::log(a);
00325 G4double sa=std::sqrt(a);
00326 G4double a2=a*a;
00327 G4double a2s=a2*sa;
00328 G4double a4=a2*a2;
00329 G4double a8=a4*a4;
00330 G4double c=(170.+3600./a2s)/(1.+65./a2s);
00331 G4double gg=42.*(std::exp(al*0.8)+4.E-8*a4)/(1.+28./a)/(1.+5.E-5*a2);
00332 G4double e=390.;
00333 G4double r=0.27;
00334 G4double h=2.E-7;
00335 G4double t=0.3;
00336 if(tZ>1 || tN>1)
00337 {
00338 e=380.+18.*a2/(1.+a2/60.)/(1.+2.E-19*a8);
00339 r=0.15;
00340 h=1.E-8*a2/(1.+a2/17.)/(1.+3.E-20*a8);
00341 t=(.2+.00056*a2)/(1.+a2*.0006);
00342 }
00343 sigma=(c+d*d)/(1.+t/ssp+r/p4)+(gg+e*std::exp(-6.*P))/(1.+h/p4/p4);
00344 #ifdef pdebug
00345 G4cout<<"G4QHyperonNucCS::CSForm: A="<<a<<",P="<<P<<",CS="<<sigma<<",c="<<c<<",g="<<gg
00346 <<",d="<<d<<",r="<<r<<",e="<<e<<",h="<<h<<G4endl;
00347 #endif
00348 }
00349 else
00350 {
00351 G4cerr<<"-Warning-G4QHyperonNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00352 sigma=0.;
00353 }
00354 if(sigma<0.) return 0.;
00355 return sigma;
00356 }
00357
00358 G4double G4ChipsHyperonInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
00359 {
00360 if(DX<=0. || N<2)
00361 {
00362 G4cerr<<"***G4ChipsHyperonInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
00363 return Y[0];
00364 }
00365
00366 G4int N2=N-2;
00367 G4double d=(X-X0)/DX;
00368 G4int j=static_cast<int>(d);
00369 if (j<0) j=0;
00370 else if(j>N2) j=N2;
00371 d-=j;
00372 G4double yi=Y[j];
00373 G4double sigma=yi+(Y[j+1]-yi)*d;
00374
00375 return sigma;
00376 }