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 "G4QKaonMinusNuclearCrossSection.hh"
00048 #include "G4SystemOfUnits.hh"
00049
00050
00051 G4double* G4QKaonMinusNuclearCrossSection::lastLEN=0;
00052 G4double* G4QKaonMinusNuclearCrossSection::lastHEN=0;
00053 G4int G4QKaonMinusNuclearCrossSection::lastN=0;
00054 G4int G4QKaonMinusNuclearCrossSection::lastZ=0;
00055 G4double G4QKaonMinusNuclearCrossSection::lastP=0.;
00056 G4double G4QKaonMinusNuclearCrossSection::lastTH=0.;
00057 G4double G4QKaonMinusNuclearCrossSection::lastCS=0.;
00058 G4int G4QKaonMinusNuclearCrossSection::lastI=0;
00059 std::vector<G4double*>* G4QKaonMinusNuclearCrossSection::LEN = new std::vector<G4double*>;
00060 std::vector<G4double*>* G4QKaonMinusNuclearCrossSection::HEN = new std::vector<G4double*>;
00061
00062
00063 G4VQCrossSection* G4QKaonMinusNuclearCrossSection::GetPointer()
00064 {
00065 static G4QKaonMinusNuclearCrossSection theCrossSection;
00066 return &theCrossSection;
00067 }
00068
00069 G4QKaonMinusNuclearCrossSection::~G4QKaonMinusNuclearCrossSection()
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 G4QKaonMinusNuclearCrossSection::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<<"G4QKmCS::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!=-321) G4cout<<"-Warning-G4QPiMinusCS::GetCS:**Not a KMinus**,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<<"G4QKmCS::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<<"G4QKmCS::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<<"..G4QKmCS::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<<"G4QKmCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00140 #endif
00141 lastCS=CalculateCrossSection(fCS,-1,j,-321,lastZ,lastN,pMom);
00142 #ifdef debug
00143 G4cout<<"G4QKmCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00144 #endif
00145 if(lastCS<=0. && pMom>lastTH)
00146 {
00147 #ifdef debug
00148 G4cout<<"G4QKmCS::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<<"-->G4QKmCrossSec::GetCrosSec: pPDG=-321, j="<<j<<", N="<<colN[i]
00157 <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
00158 #endif
00159 j++;
00160 }
00161 #ifdef debug
00162 G4cout<<"-?-G4QKmCS::GetCS:RC Z="<<tgZ<<",N="<<tgN<<",in="<<in<<",j="<<j<<" ?"<<G4endl;
00163 #endif
00164 if(!in)
00165 {
00166 #ifdef debug
00167 G4cout<<"^^^G4QKmCS::GetCS:CalcNew P="<<pMom<<", f="<<fCS<<", lastI="<<lastI<<G4endl;
00168 #endif
00170 lastCS=CalculateCrossSection(fCS,0,j,-321,lastZ,lastN,pMom); //calculate & create
00171
00172
00173
00174 lastTH = ThresholdEnergy(tgZ, tgN);
00175 #ifdef debug
00176 G4cout<<"G4QKmCrossSection::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<<"G4QKmCS::GetCrosSec:recCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00185 #endif
00186
00187 #ifdef pdebug
00188 G4cout<<"G4QKmCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00189 #endif
00190 return lastCS*millibarn;
00191 }
00192 else
00193 {
00194 #ifdef debug
00195 G4cout<<"G4QKmCS::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<<"G4QKmCS::GetCS: Current P="<<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,-321,lastZ,lastN,pMom);
00222 lastP=pMom;
00223 }
00224 #ifdef debug
00225 G4cout<<"==>G4QKmCS::GetCroSec: P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00226 #endif
00227 return lastCS*millibarn;
00228 }
00229
00230
00231 G4double G4QKaonMinusNuclearCrossSection::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<<"G4QKmNCS::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<<"G4QKmNucCS::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<<"-*->G4QKmNucCS::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<<"***G4QPiMinusNuclCS::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<<"G4QKmNCS::CalcCS:lTH="<<lastTH<<",Pmi="<<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<<"G4QKmNCS::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<<"G4QKmNuCS::CalcCS: E="<<Momentum<<",T="<<THmin<<",dP="<<dP<<G4endl;
00311 #endif
00312 }
00313 else if (Momentum<Pmax)
00314 {
00315 G4double lP=std::log(Momentum);
00316 #ifdef debug
00317 G4cout<<"G4QKmNucCS::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<<"G4QKaonMinusNuclearCrossSection::CalcCS: CS="<<sigma<<G4endl;
00328 #endif
00329 if(sigma<0.) return 0.;
00330 return sigma;
00331 }
00332
00333
00334 G4double G4QKaonMinusNuclearCrossSection::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00335 {
00336 G4double lP=std::log(P);
00337 return CrossSectionFormula(tZ, tN, P, lP);
00338 }
00339
00340
00341 G4double G4QKaonMinusNuclearCrossSection::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00342 {
00343 G4double P=std::exp(lP);
00344 return CrossSectionFormula(tZ, tN, P, lP);
00345 }
00346
00347 G4double G4QKaonMinusNuclearCrossSection::CrossSectionFormula(G4int tZ, G4int tN,
00348 G4double P, G4double lP)
00349 {
00350 G4double sigma=0.;
00351 if(tZ==1 && !tN)
00352 {
00353 G4double ld=lP-3.5;
00354 G4double ld2=ld*ld;
00355 G4double p2=P*P;
00356 G4double p4=p2*p2;
00357 G4double sp=std::sqrt(P);
00358 G4double psp=P*sp;
00359 G4double lm=P-.39;
00360 G4double md=lm*lm+.000156;
00361 G4double lh=P-1.;
00362 G4double hd=lh*lh+.0156;
00363 G4double El=(.0557*ld2+2.23)/(1.-.7/sp+.075/p4);
00364 G4double To=(.3*ld2+19.5)/(1.-.21/sp+.52/p4);
00365 sigma=8.8/psp+(To-El)+.002/md+.15/hd;
00366 }
00367 else if(tZ==1 && tN==1)
00368 {
00369 G4double p2=P*P;
00370 G4double dX=lP-3.7;
00371 G4double dR=P-.94;
00372 G4double sp=std::sqrt(P);
00373 sigma=(.6*dX*dX+36.)/(1.-.11/sp+.52/p2/p2)+.7/(dR*dR+.0256)+18./P/sp;
00374 }
00375 else if(tZ<97 && tN<152)
00376 {
00377 G4double d=lP-4.2;
00378 G4double sp=std::sqrt(P);
00379 G4double p2=P*P;
00380 G4double a=tN+tZ;
00381 G4double sa=std::sqrt(a);
00382 G4double al=std::log(a);
00383 G4double a2=a*a;
00384 G4double c=52.*std::exp(al*0.6)*(1.+97./a2)/(1.+9.8/a)/(1.+47./a2);
00385 G4double g_value=-.2-.003*a;
00386 G4double h=.5+.07*a;
00387 G4double v=P-1.;
00388 G4double f=.6*a*sa/(1.+.00002*a2);
00389 G4double u=.125+.127*al;
00390 sigma=(c+d*d)/(1.+g_value/sp+h/p2/p2)+f/(v*v+u*u)+20.*sa/P/sp;
00391 #ifdef pdebug
00392 G4cout<<"G4QKMinNucCS::CSForm: A="<<a<<",P="<<P<<",CS="<<sigma<<",c="<<c<<",g="<<g_value
00393 <<",d="<<d<<",r="<<r<<",e="<<e<<",h="<<h<<G4endl;
00394 #endif
00395 }
00396 else
00397 {
00398 G4cerr<<"-Warning-G4QKMinusNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00399 sigma=0.;
00400 }
00401 if(sigma<0.) return 0.;
00402 return sigma;
00403 }