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 "G4QAntiBaryonNuclearCrossSection.hh"
00048 #include "G4SystemOfUnits.hh"
00049
00050
00051 G4double* G4QAntiBaryonNuclearCrossSection::lastLEN=0;
00052 G4double* G4QAntiBaryonNuclearCrossSection::lastHEN=0;
00053 G4int G4QAntiBaryonNuclearCrossSection::lastN=0;
00054 G4int G4QAntiBaryonNuclearCrossSection::lastZ=0;
00055 G4double G4QAntiBaryonNuclearCrossSection::lastP=0.;
00056 G4double G4QAntiBaryonNuclearCrossSection::lastTH=0.;
00057 G4double G4QAntiBaryonNuclearCrossSection::lastCS=0.;
00058 G4int G4QAntiBaryonNuclearCrossSection::lastI=0;
00059 std::vector<G4double*>* G4QAntiBaryonNuclearCrossSection::LEN = new std::vector<G4double*>;
00060 std::vector<G4double*>* G4QAntiBaryonNuclearCrossSection::HEN = new std::vector<G4double*>;
00061
00062
00063 G4VQCrossSection* G4QAntiBaryonNuclearCrossSection::GetPointer()
00064 {
00065 static G4QAntiBaryonNuclearCrossSection theCrossSection;
00066 return &theCrossSection;
00067 }
00068
00069 G4QAntiBaryonNuclearCrossSection::~G4QAntiBaryonNuclearCrossSection()
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 G4QAntiBaryonNuclearCrossSection::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<<"G4QaBCS::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>-2111 || PDG==-3112 || PDG==-3312 || PDG==-3334)
00097 G4cout<<"-Warning-G4QAntiBaryonCS::GetCS: Not Negat/Zero AntiBaryon,PDG="<<PDG<<G4endl;
00098 G4bool in=false;
00099 if(tgN!=lastN || tgZ!=lastZ)
00100 {
00101 in = false;
00102 lastP = 0.;
00103 lastN = tgN;
00104 lastZ = tgZ;
00105 lastI = colN.size();
00106 j = 0;
00107 #ifdef debug
00108 G4cout<<"G4QaBarNucCS::GetCS: the amount of records in the AMDB lastI="<<lastI<<G4endl;
00109 #endif
00110 if(lastI) for(G4int i=0; i<lastI; i++)
00111 {
00112 if(colN[i]==tgN && colZ[i]==tgZ)
00113 {
00114 lastI=i;
00115 lastTH =colTH[i];
00116 #ifdef debug
00117 G4cout<<"G4QaBCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00118 #endif
00119 if(pMom<=lastTH)
00120 {
00121 #ifdef debug
00122 G4cout<<"G4QPCS::GetCS:Found,P="<<pMom<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
00123 #endif
00124 return 0.;
00125 }
00126 lastP =colP [i];
00127 lastCS =colCS[i];
00128 if(std::fabs(lastP-pMom)<tolerance*pMom)
00129
00130 {
00131 #ifdef debug
00132 G4cout<<"..G4QaBCS::GetCS:.DoNothing.P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00133 #endif
00134
00135 return lastCS*millibarn;
00136 }
00137 in = true;
00138
00139 #ifdef debug
00140 G4cout<<"G4QaBCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00141 #endif
00142 lastCS=CalculateCrossSection(fCS,-1,j,PDG,lastZ,lastN,pMom);
00143 #ifdef debug
00144 G4cout<<"G4QaBCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00145 #endif
00146 if(lastCS<=0. && pMom>lastTH)
00147 {
00148 #ifdef debug
00149 G4cout<<"G4QaBNucCS::GetCS: New P="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00150 #endif
00151 lastCS=0.;
00152 lastTH=pMom;
00153 }
00154 break;
00155 }
00156 #ifdef debug
00157 G4cout<<"-->G4QaBNucCrossSec::GetCrosSec: pPDG="<<PDG<<", j="<<j<<", N="<<colN[i]
00158 <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
00159 #endif
00160 j++;
00161 }
00162 #ifdef debug
00163 G4cout<<"-?-G4QaBCS::GetCS:RC Z="<<tgZ<<",N="<<tgN<<",in="<<in<<",j="<<j<<" ?"<<G4endl;
00164 #endif
00165 if(!in)
00166 {
00167 #ifdef debug
00168 G4cout<<"^^^G4QaBCS::GetCS:CalcNew P="<<pMom<<", f="<<fCS<<", lastI="<<lastI<<G4endl;
00169 #endif
00171 lastCS=CalculateCrossSection(fCS,0,j,PDG,lastZ,lastN,pMom); //calculate & create
00172
00173
00174
00175 lastTH = ThresholdEnergy(tgZ, tgN);
00176 #ifdef debug
00177 G4cout<<"G4QaBCrossSection::GetCrossSect: NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
00178 #endif
00179 colN.push_back(tgN);
00180 colZ.push_back(tgZ);
00181 colP.push_back(pMom);
00182 colTH.push_back(lastTH);
00183 colCS.push_back(lastCS);
00184 #ifdef debug
00185 G4cout<<"G4QaBCS::GetCrosSec:recCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00186 #endif
00187
00188 #ifdef pdebug
00189 G4cout<<"G4QaBCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00190 #endif
00191 return lastCS*millibarn;
00192 }
00193 else
00194 {
00195 #ifdef debug
00196 G4cout<<"G4QAntiBaryNucCrossSection::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00197 #endif
00198 colP[lastI]=pMom;
00199 colCS[lastI]=lastCS;
00200 }
00201 }
00202 else if(pMom<=lastTH)
00203 {
00204 #ifdef debug
00205 G4cout<<"G4QaBNCS::GetCS: Current P="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00206 #endif
00207 return 0.;
00208 }
00209 else if(std::fabs(lastP-pMom)<tolerance*pMom)
00210
00211 {
00212 #ifdef debug
00213 G4cout<<"..G4QPCS::GetCS:OldNZ&P="<<lastP<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00214 #endif
00215 return lastCS*millibarn;
00216 }
00217 else
00218 {
00219 #ifdef debug
00220 G4cout<<"-!-G4QPCS::GetCS:UseCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00221 #endif
00222 lastCS=CalculateCrossSection(fCS,1,j,PDG,lastZ,lastN,pMom);
00223 lastP=pMom;
00224 }
00225 #ifdef debug
00226 G4cout<<"==>G4QaBCS::GetCroSec: P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00227 #endif
00228 return lastCS*millibarn;
00229 }
00230
00231
00232 G4double G4QAntiBaryonNuclearCrossSection::CalculateCrossSection(G4bool, G4int F, G4int I,
00233 G4int, G4int targZ, G4int targN, G4double Momentum)
00234 {
00235 static const G4double THmin=27.;
00236 static const G4double THmiG=THmin*.001;
00237 static const G4double dP=10.;
00238 static const G4double dPG=dP*.001;
00239 static const G4int nL=105;
00240 static const G4double Pmin=THmin+(nL-1)*dP;
00241 static const G4double Pmax=227000.;
00242 static const G4int nH=224;
00243 static const G4double milP=std::log(Pmin);
00244 static const G4double malP=std::log(Pmax);
00245 static const G4double dlP=(malP-milP)/(nH-1);
00246 static const G4double milPG=std::log(.001*Pmin);
00247 #ifdef debug
00248 G4cout<<"G4QaBarNCS::CalCS:N="<<targN<<",Z="<<targZ<<",P="<<Momentum<<">"<<THmin<<G4endl;
00249 #endif
00250 G4double sigma=0.;
00251 if(F&&I) sigma=0.;
00252
00253 #ifdef debug
00254 G4cout<<"G4QaBarNucCS::CalCS: A="<<A<<",F="<<F<<",I="<<I<<",nL="<<nL<<",nH="<<nH<<G4endl;
00255 #endif
00256 if(F<=0)
00257 {
00258 if(F<0)
00259 {
00260 G4int sync=LEN->size();
00261 if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
00262 lastLEN=(*LEN)[I];
00263 lastHEN=(*HEN)[I];
00264 }
00265 else
00266 {
00267 lastLEN = new G4double[nL];
00268 lastHEN = new G4double[nH];
00269
00270 G4double P=THmiG;
00271 for(G4int n=0; n<nL; n++)
00272 {
00273 lastLEN[n] = CrossSectionLin(targZ, targN, P);
00274 P+=dPG;
00275 }
00276 G4double lP=milPG;
00277 for(G4int n=0; n<nH; n++)
00278 {
00279 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
00280 lP+=dlP;
00281 }
00282 #ifdef debug
00283 G4cout<<"-*->G4QaBarNucCS::CalcCS:Tab for Z="<<targZ<<",N="<<targN<<",I="<<I<<G4endl;
00284 #endif
00285
00286
00287 G4int sync=LEN->size();
00288 if(sync!=I)
00289 {
00290 G4cerr<<"***G4QPiMinusNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
00291 <<", N="<<targN<<", F="<<F<<G4endl;
00292
00293 }
00294 LEN->push_back(lastLEN);
00295 HEN->push_back(lastHEN);
00296 }
00297 }
00298
00299 #ifdef debug
00300 G4cout<<"G4QaBNCS::CalcCS:lTH="<<lastTH<<",Pmi="<<Pmin<<",dP="<<dP<<",dlP="<<dlP<<G4endl;
00301 #endif
00302 if (Momentum<lastTH) return 0.;
00303 else if (Momentum<Pmin)
00304 {
00305 #ifdef debug
00306 G4cout<<"G4QaBNCS::CalcCS:bLEN nL="<<nL<<",TH="<<THmin<<",dP="<<dP<<G4endl;
00307 #endif
00308 sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
00309 #ifdef debugn
00310 if(sigma<0.)
00311 G4cout<<"G4QaBNuCS::CalcCS: E="<<Momentum<<",T="<<THmin<<",dP="<<dP<<G4endl;
00312 #endif
00313 }
00314 else if (Momentum<Pmax)
00315 {
00316 G4double lP=std::log(Momentum);
00317 #ifdef debug
00318 G4cout<<"G4QaBarNucCS::CalcCS: before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl;
00319 #endif
00320 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
00321 }
00322 else
00323 {
00324 G4double P=0.001*Momentum;
00325 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
00326 }
00327 #ifdef debug
00328 G4cout<<"G4QAntiBaryonNuclearCrossSection::CalcCS: CS="<<sigma<<G4endl;
00329 #endif
00330 if(sigma<0.) return 0.;
00331 return sigma;
00332 }
00333
00334
00335 G4double G4QAntiBaryonNuclearCrossSection::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00336 {
00337 G4double lP=std::log(P);
00338 return CrossSectionFormula(tZ, tN, P, lP);
00339 }
00340
00341
00342 G4double G4QAntiBaryonNuclearCrossSection::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00343 {
00344 G4double P=std::exp(lP);
00345 return CrossSectionFormula(tZ, tN, P, lP);
00346 }
00347
00348 G4double G4QAntiBaryonNuclearCrossSection::CrossSectionFormula(G4int tZ, G4int tN,
00349 G4double P, G4double lP)
00350 {
00351 G4double sigma=0.;
00352 if(tZ==1 && !tN)
00353 {
00354 G4double ld=lP-3.5;
00355 G4double ld2=ld*ld;
00356 G4double ye=std::exp(lP*1.25);
00357 G4double yt=std::exp(lP*0.35);
00358 G4double El=80./(ye+1.);
00359 G4double To=(80./yt+.3)/yt;
00360 sigma=(To-El)+.2443*ld2+31.48;
00361 }
00362 else if(tZ==1 && tN==1)
00363 {
00364 G4double r=lP-3.7;
00365 sigma=0.6*r*r+67.+90.*std::exp(-lP*.666);
00366 }
00367 else if(tZ<97 && tN<152)
00368 {
00369 G4double d=lP-4.2;
00370 G4double sp=std::sqrt(P);
00371 G4double a=tN+tZ;
00372 G4double sa=std::sqrt(a);
00373 G4double a2=a*a;
00374 G4double a3=a2*a;
00375 G4double a2s=a2*sa;
00376 G4double c=(170.+3600./a2s)/(1.+65./a2s)+40.*std::pow(a,0.712)/(1.+12.2/a)/(1.+34./a2);
00377 G4double r=(170.+0.01*a3)/(1.+a3/28000.);
00378 sigma=c+d*d+r/sp;
00379 #ifdef pdebug
00380 G4cout<<"G4QAntiBarNucCS::CSForm: A="<<a<<",P="<<P<<",CS="<<sigma<<",c="<<c<<",g="<<g
00381 <<",d="<<d<<",r="<<r<<",e="<<e<<",h="<<h<<G4endl;
00382 #endif
00383 }
00384 else
00385 {
00386 G4cerr<<"-Warning-G4QAntiBarNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00387 sigma=0.;
00388 }
00389 if(sigma<0.) return 0.;
00390 return sigma;
00391 }