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