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
00042
00043
00044
00045
00046
00047 #include "G4QHyperonPlusNuclearCrossSection.hh"
00048 #include "G4SystemOfUnits.hh"
00049
00050
00051 G4double* G4QHyperonPlusNuclearCrossSection::lastLEN=0;
00052 G4double* G4QHyperonPlusNuclearCrossSection::lastHEN=0;
00053 G4int G4QHyperonPlusNuclearCrossSection::lastN=0;
00054 G4int G4QHyperonPlusNuclearCrossSection::lastZ=0;
00055 G4double G4QHyperonPlusNuclearCrossSection::lastP=0.;
00056 G4double G4QHyperonPlusNuclearCrossSection::lastTH=0.;
00057 G4double G4QHyperonPlusNuclearCrossSection::lastCS=0.;
00058 G4int G4QHyperonPlusNuclearCrossSection::lastI=0;
00059 std::vector<G4double*>* G4QHyperonPlusNuclearCrossSection::LEN= new std::vector<G4double*>;
00060 std::vector<G4double*>* G4QHyperonPlusNuclearCrossSection::HEN= new std::vector<G4double*>;
00061
00062
00063 G4VQCrossSection* G4QHyperonPlusNuclearCrossSection::GetPointer()
00064 {
00065 static G4QHyperonPlusNuclearCrossSection theCrossSection;
00066 return &theCrossSection;
00067 }
00068
00069 G4QHyperonPlusNuclearCrossSection::~G4QHyperonPlusNuclearCrossSection()
00070 {
00071 G4int lens=LEN->size();
00072 for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
00073 delete LEN;
00074 G4int hens=HEN->size();
00075 for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
00076 delete HEN;
00077 }
00078
00079
00080
00081 G4double G4QHyperonPlusNuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
00082 G4int tgZ, G4int tgN, G4int PDG)
00083 {
00084
00085 static G4int j;
00086 static std::vector <G4int> colN;
00087 static std::vector <G4int> colZ;
00088 static std::vector <G4double> colP;
00089 static std::vector <G4double> colTH;
00090 static std::vector <G4double> colCS;
00091
00092 #ifdef debug
00093 G4cout<<"G4QHyPNCS::GetCS:>>>f="<<fCS<<",p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
00094 <<"("<<lastN<<"), PDG="<<PDG<<", thresh="<<lastTH<<",Sz="<<colN.size()<<G4endl;
00095 #endif
00096 if(PDG!=3222) G4cout<<"-Warning-G4QPiMinusCS::GetCS:**Not a PiMinus**,PDG="<<PDG<<G4endl;
00097 G4bool in=false;
00098 if(tgN!=lastN || tgZ!=lastZ)
00099 {
00100 in = false;
00101 lastP = 0.;
00102 lastN = tgN;
00103 lastZ = tgZ;
00104 lastI = colN.size();
00105 j = 0;
00106 #ifdef debug
00107 G4cout<<"G4QHyPNuclCS::GetCS: the amount of records in the AMDB lastI="<<lastI<<G4endl;
00108 #endif
00109 if(lastI) for(G4int i=0; i<lastI; i++)
00110 {
00111 if(colN[i]==tgN && colZ[i]==tgZ)
00112 {
00113 lastI=i;
00114 lastTH =colTH[i];
00115 #ifdef debug
00116 G4cout<<"G4QHyPCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00117 #endif
00118 if(pMom<=lastTH)
00119 {
00120 #ifdef debug
00121 G4cout<<"G4QPCS::GetCS:Found,P="<<pMom<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
00122 #endif
00123 return 0.;
00124 }
00125 lastP =colP [i];
00126 lastCS =colCS[i];
00127 if(std::fabs(lastP-pMom)<tolerance*pMom)
00128
00129 {
00130 #ifdef debug
00131 G4cout<<"G4QHyPNCS::GetCS:.DoNothing.P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00132 #endif
00133
00134 return lastCS*millibarn;
00135 }
00136 in = true;
00137
00138 #ifdef debug
00139 G4cout<<"G4QHyPNCS::G:UpdDB,P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00140 #endif
00141 lastCS=CalculateCrossSection(fCS,-1,j,PDG,lastZ,lastN,pMom);
00142 #ifdef debug
00143 G4cout<<"G4QHyPNuCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00144 #endif
00145 if(lastCS<=0. && pMom>lastTH)
00146 {
00147 #ifdef debug
00148 G4cout<<"G4QHyPNuCS::GetCS: New P="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00149 #endif
00150 lastCS=0.;
00151 lastTH=pMom;
00152 }
00153 break;
00154 }
00155 #ifdef debug
00156 G4cout<<"-->G4QHyPNucCrossSec::GetCrosSec: pPDG="<<PDG<<", j="<<j<<", N="<<colN[i]
00157 <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
00158 #endif
00159 j++;
00160 }
00161 #ifdef debug
00162 G4cout<<"-?-G4QHyPNCS::GetCS:R,Z="<<tgZ<<",N="<<tgN<<",i="<<in<<",j="<<j<<" ?"<<G4endl;
00163 #endif
00164 if(!in)
00165 {
00166 #ifdef debug
00167 G4cout<<"^G4QHyPNCS::GetCS:CalcNew P="<<pMom<<", f="<<fCS<<", lastI="<<lastI<<G4endl;
00168 #endif
00170 lastCS=CalculateCrossSection(fCS,0,j,PDG,lastZ,lastN,pMom); //calculate & create
00171
00172
00173
00174 lastTH = ThresholdEnergy(tgZ, tgN);
00175 #ifdef debug
00176 G4cout<<"G4QHypNucCrossSect::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
00177 #endif
00178 colN.push_back(tgN);
00179 colZ.push_back(tgZ);
00180 colP.push_back(pMom);
00181 colTH.push_back(lastTH);
00182 colCS.push_back(lastCS);
00183 #ifdef debug
00184 G4cout<<"G4QHyPNCS::GetCrosSec:lCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00185 #endif
00186
00187 #ifdef pdebug
00188 G4cout<<"G4QHyPCS::GeCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00189 #endif
00190 return lastCS*millibarn;
00191 }
00192 else
00193 {
00194 #ifdef debug
00195 G4cout<<"G4QHyPNucCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00196 #endif
00197 colP[lastI]=pMom;
00198 colCS[lastI]=lastCS;
00199 }
00200 }
00201 else if(pMom<=lastTH)
00202 {
00203 #ifdef debug
00204 G4cout<<"G4QHyPNCS::GetCS: CurrentP="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00205 #endif
00206 return 0.;
00207 }
00208 else if(std::fabs(lastP-pMom)<tolerance*pMom)
00209
00210 {
00211 #ifdef debug
00212 G4cout<<"..G4QPCS::GetCS:OldNZ&P="<<lastP<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00213 #endif
00214 return lastCS*millibarn;
00215 }
00216 else
00217 {
00218 #ifdef debug
00219 G4cout<<"-!-G4QPCS::GetCS:UseCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00220 #endif
00221 lastCS=CalculateCrossSection(fCS,1,j,PDG,lastZ,lastN,pMom);
00222 lastP=pMom;
00223 }
00224 #ifdef debug
00225 G4cout<<"==>G4QHyPNCS::GetCrosS:P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00226 #endif
00227 return lastCS*millibarn;
00228 }
00229
00230
00231 G4double G4QHyperonPlusNuclearCrossSection::CalculateCrossSection(G4bool, G4int F, G4int I,
00232 G4int, G4int targZ, G4int targN, G4double Momentum)
00233 {
00234 static const G4double THmin=27.;
00235 static const G4double THmiG=THmin*.001;
00236 static const G4double dP=10.;
00237 static const G4double dPG=dP*.001;
00238 static const G4int nL=105;
00239 static const G4double Pmin=THmin+(nL-1)*dP;
00240 static const G4double Pmax=227000.;
00241 static const G4int nH=224;
00242 static const G4double milP=std::log(Pmin);
00243 static const G4double malP=std::log(Pmax);
00244 static const G4double dlP=(malP-milP)/(nH-1);
00245 static const G4double milPG=std::log(.001*Pmin);
00246 #ifdef debug
00247 G4cout<<"G4QHyPNuCS::CalCS:N="<<targN<<",Z="<<targZ<<",P="<<Momentum<<">"<<THmin<<G4endl;
00248 #endif
00249 G4double sigma=0.;
00250 if(F&&I) sigma=0.;
00251
00252 #ifdef debug
00253 G4cout<<"G4QHyPNucCS::CalCS: A="<<A<<", F="<<F<<",I="<<I<<",nL="<<nL<<",nH="<<nH<<G4endl;
00254 #endif
00255 if(F<=0)
00256 {
00257 if(F<0)
00258 {
00259 G4int sync=LEN->size();
00260 if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
00261 lastLEN=(*LEN)[I];
00262 lastHEN=(*HEN)[I];
00263 }
00264 else
00265 {
00266 lastLEN = new G4double[nL];
00267 lastHEN = new G4double[nH];
00268
00269 G4double P=THmiG;
00270 for(G4int n=0; n<nL; n++)
00271 {
00272 lastLEN[n] = CrossSectionLin(targZ, targN, P);
00273 P+=dPG;
00274 }
00275 G4double lP=milPG;
00276 for(G4int n=0; n<nH; n++)
00277 {
00278 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
00279 lP+=dlP;
00280 }
00281 #ifdef debug
00282 G4cout<<"-*->G4QHyPNucCS::CalcCS: Tab for Z="<<targZ<<",N="<<targN<<",I="<<I<<G4endl;
00283 #endif
00284
00285
00286 G4int sync=LEN->size();
00287 if(sync!=I)
00288 {
00289 G4cerr<<"***G4QHyPNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
00290 <<", N="<<targN<<", F="<<F<<G4endl;
00291
00292 }
00293 LEN->push_back(lastLEN);
00294 HEN->push_back(lastHEN);
00295 }
00296 }
00297
00298 #ifdef debug
00299 G4cout<<"G4QHyPNCS::CalcCS:lTH="<<lastTH<<",Pi="<<Pmin<<",dP="<<dP<<",dlP="<<dlP<<G4endl;
00300 #endif
00301 if (Momentum<lastTH) return 0.;
00302 else if (Momentum<Pmin)
00303 {
00304 #ifdef debug
00305 G4cout<<"G4QHyPNCS::CalcCS:bLEN nL="<<nL<<",TH="<<THmin<<",dP="<<dP<<G4endl;
00306 #endif
00307 sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
00308 #ifdef debugn
00309 if(sigma<0.)
00310 G4cout<<"G4QHyPNuCS::CalcCS: E="<<Momentum<<",T="<<THmin<<",d="<<dP<<G4endl;
00311 #endif
00312 }
00313 else if (Momentum<Pmax)
00314 {
00315 G4double lP=std::log(Momentum);
00316 #ifdef debug
00317 G4cout<<"G4QHyPNucCS::CalcCS: before HEN nH="<<nH<<", iE="<<milP<<",dlP="<<dlP<<G4endl;
00318 #endif
00319 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
00320 }
00321 else
00322 {
00323 G4double P=0.001*Momentum;
00324 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
00325 }
00326 #ifdef debug
00327 G4cout<<"G4QHyperonPlusNuclearCrossSection::CalcCS: CS="<<sigma<<G4endl;
00328 #endif
00329 if(sigma<0.) return 0.;
00330 return sigma;
00331 }
00332
00333
00334 G4double G4QHyperonPlusNuclearCrossSection::ThresholdMomentum(G4int tZ, G4int tN)
00335 {
00336 static const G4double third=1./3.;
00337 static const G4double prM = G4QPDGCode(2212).GetMass();
00338 static const G4double pM = G4QPDGCode(3112).GetMass();
00339 static const G4double tpM= pM+pM;
00340 G4double tA=tZ+tN;
00341 if(tZ<.99 || tN<0.) return 0.;
00342 G4double tM=931.5*tA;
00343 G4double dE=65.88;
00344 if(tZ==1 && tN==0) tM=prM;
00345 else dE=tZ/(1.+std::pow(tA,third));
00346
00347 G4double T=dE+dE*(dE/2+pM)/tM;
00348 return std::sqrt(T*(tpM+T));
00349 }
00350
00351
00352 G4double G4QHyperonPlusNuclearCrossSection::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00353 {
00354 G4double lP=std::log(P);
00355 return CrossSectionFormula(tZ, tN, P, lP);
00356 }
00357
00358
00359 G4double G4QHyperonPlusNuclearCrossSection::CrossSectionLog(G4int tZ, G4int tN,G4double lP)
00360 {
00361 G4double P=std::exp(lP);
00362 return CrossSectionFormula(tZ, tN, P, lP);
00363 }
00364
00365 G4double G4QHyperonPlusNuclearCrossSection::CrossSectionFormula(G4int tZ, G4int tN,
00366 G4double P, G4double lP)
00367 {
00368 G4double sigma=0.;
00369 if(tZ==1 && !tN)
00370 {
00371 G4double ld=lP-3.5;
00372 G4double ld2=ld*ld;
00373 G4double p2=P*P;
00374 G4double p4=p2*p2;
00375 G4double sp=std::sqrt(P);
00376 G4double El=(.0557*ld2+6.72+99./p2)/(1.+2./sp+2./p4);
00377 G4double To=(.3*ld2+38.2+900./sp)/(1.+27./sp+3./p4);
00378 sigma=To-El;
00379 }
00380 else if(tZ<97 && tN<152)
00381 {
00382 G4double d=lP-4.2;
00383 G4double p2=P*P;
00384 G4double p4=p2*p2;
00385 G4double sp=std::sqrt(P);
00386 G4double ssp=std::sqrt(sp);
00387 G4double a=tN+tZ;
00388 G4double al=std::log(a);
00389 G4double sa=std::sqrt(a);
00390 G4double a2=a*a;
00391 G4double a2s=a2*sa;
00392 G4double a4=a2*a2;
00393 G4double a8=a4*a4;
00394 G4double c=(170.+3600./a2s)/(1.+65./a2s);
00395 G4double g_value=40.*std::exp(al*0.712)/(1.+12.2/a)/(1.+34./a2);
00396 G4double e=700.;
00397 G4double r=0.27;
00398 G4double h=3.E-7;
00399 G4double t=0.3;
00400 if(tZ>1 || tN>1)
00401 {
00402 e=99.+146.*a2/(1.+a2/12.)/(1.+2.E-19*a8);
00403 r=0.15;
00404 h=(.0017+.2E-9*a4*sa)/a2/a;
00405 t=(.2+.00056*a2)/(1.+a2*.0006);
00406 }
00407 sigma=(c+d*d)/(1.+t/ssp+r/p4)+(g_value+e*std::exp(-6.*P))/(1.+h/p4/p4);
00408 #ifdef pdebug
00409 G4cout<<"G4QHyperPlusNucCS::CSForm: A="<<a<<",P="<<P<<",CS="<<sigma<<",c="<<c<<",g="<<g_value
00410 <<",d="<<d<<",r="<<r<<",e="<<e<<",h="<<h<<G4endl;
00411 #endif
00412 }
00413 else
00414 {
00415 G4cerr<<"-Warning-G4QHyperPlNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00416 sigma=0.;
00417 }
00418 if(sigma<0.) return 0.;
00419 return sigma;
00420 }