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 "G4QPionPlusNuclearCrossSection.hh"
00048 #include "G4SystemOfUnits.hh"
00049
00050
00051 G4double* G4QPionPlusNuclearCrossSection::lastLEN=0;
00052 G4double* G4QPionPlusNuclearCrossSection::lastHEN=0;
00053 G4int G4QPionPlusNuclearCrossSection::lastN=0;
00054 G4int G4QPionPlusNuclearCrossSection::lastZ=0;
00055 G4double G4QPionPlusNuclearCrossSection::lastP=0.;
00056 G4double G4QPionPlusNuclearCrossSection::lastTH=0.;
00057 G4double G4QPionPlusNuclearCrossSection::lastCS=0.;
00058 G4int G4QPionPlusNuclearCrossSection::lastI=0;
00059 std::vector<G4double*>* G4QPionPlusNuclearCrossSection::LEN = new std::vector<G4double*>;
00060 std::vector<G4double*>* G4QPionPlusNuclearCrossSection::HEN = new std::vector<G4double*>;
00061
00062
00063 G4VQCrossSection* G4QPionPlusNuclearCrossSection::GetPointer()
00064 {
00065 static G4QPionPlusNuclearCrossSection theCrossSection;
00066 return &theCrossSection;
00067 }
00068
00069 G4QPionPlusNuclearCrossSection::~G4QPionPlusNuclearCrossSection()
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 G4QPionPlusNuclearCrossSection::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<<"G4QPipCS::GetCS:>>>f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
00094 <<"("<<lastN<<"),PDG=211, thresh="<<lastTH<<",Sz="<<colN.size()<<G4endl;
00095 #endif
00096 if(PDG!=211) G4cout<<"-Warning-G4QPionPlusCS::GetCS:**Not a PiPlus**,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<<"G4QPipCS::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<<"G4QPipCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00117 #endif
00118 if(pMom<=lastTH)
00119 {
00120 #ifdef debug
00121 G4cout<<"G4QPipCS::GCS: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<<"..G4QPipCS::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<<"G4QPipCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00140 #endif
00141 lastCS=CalculateCrossSection(fCS,-1,j,211,lastZ,lastN,pMom);
00142 #ifdef debug
00143 G4cout<<"G4QPipNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00144 #endif
00145 if(lastCS<=0. && pMom>lastTH)
00146 {
00147 #ifdef debug
00148 G4cout<<"G4QPipNCS::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<<"-->G4QPipNucCrossSec::GetCrosSec: pPDG=211, j="<<j<<", N="<<colN[i]
00157 <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
00158 #endif
00159 j++;
00160 }
00161 #ifdef debug
00162 G4cout<<"-?-G4QPipCS::GetCS:R,Z="<<tgZ<<",N="<<tgN<<",in="<<in<<",j="<<j<<" ?"<<G4endl;
00163 #endif
00164 if(!in)
00165 {
00166 #ifdef debug
00167 G4cout<<"^^^G4QPipCS::GeCS:CalcNew P="<<pMom<<", f="<<fCS<<", lastI="<<lastI<<G4endl;
00168 #endif
00170 lastCS=CalculateCrossSection(fCS,0,j,211,lastZ,lastN,pMom); //calculate & create
00171
00172
00173
00174 lastTH = ThresholdEnergy(tgZ, tgN);
00175 #ifdef debug
00176 G4cout<<"G4QPipCrossSection::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<<"G4QPipNCS::GetCrosSec:lCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00185 #endif
00186
00187 #ifdef pdebug
00188 G4cout<<"G4QPipCS::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<<"G4QPipNucCS::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<<"G4QPipCS::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,211,lastZ,lastN,pMom);
00222 lastP=pMom;
00223 }
00224 #ifdef debug
00225 G4cout<<"==>G4QPipCS::GetCroSec:P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00226 #endif
00227 return lastCS*millibarn;
00228 }
00229
00230
00231 G4double G4QPionPlusNuclearCrossSection::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<<"G4QPipNuCS::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<<"G4QPipNucCS::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<<"-*->G4QPipNucCS::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<<"G4QPipNCS::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<<"G4QPipNCS::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<<"G4QPipNuCS::CalCS: 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<<"G4QPipNucCS::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<<"G4QPionPlusNuclearCrossSection::CalcCS: CS="<<sigma<<G4endl;
00328 #endif
00329 if(sigma<0.) return 0.;
00330 return sigma;
00331 }
00332
00333
00334 G4double G4QPionPlusNuclearCrossSection::ThresholdMomentum(G4int tZ, G4int tN)
00335 {
00336 static const G4double third=1./3.;
00337 static const G4double pM = G4QPDGCode(211).GetMass();
00338 static const G4double tpM= pM+pM;
00339 G4double tA=tZ+tN;
00340 if(tZ<.99 || tN<0.) return 0.;
00341 else if(tZ==1 && tN==0) return 300.;
00342
00343 G4double dE=tZ/(1.+std::pow(tA,third));
00344 G4double tM=931.5*tA;
00345 G4double T=dE+dE*(dE/2+pM)/tM;
00346 return std::sqrt(T*(tpM+T));
00347 }
00348
00349
00350 G4double G4QPionPlusNuclearCrossSection::CrossSectionLin(G4int tZ, G4int tN, G4double P)
00351 {
00352 G4double lP=std::log(P);
00353 return CrossSectionFormula(tZ, tN, P, lP);
00354 }
00355
00356
00357 G4double G4QPionPlusNuclearCrossSection::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
00358 {
00359 G4double P=std::exp(lP);
00360 return CrossSectionFormula(tZ, tN, P, lP);
00361 }
00362
00363 G4double G4QPionPlusNuclearCrossSection::CrossSectionFormula(G4int tZ, G4int tN,
00364 G4double P, G4double lP)
00365 {
00366 G4double sigma=0.;
00367 if(tZ==1 && !tN)
00368 {
00369 G4double ld=lP-3.5;
00370 G4double ld2=ld*ld;
00371 G4double p2=P*P;
00372 G4double p4=p2*p2;
00373 G4double sp=std::sqrt(P);
00374 G4double lm=lP-.32;
00375 G4double md=lm*lm+.04;
00376 G4double El=(.0557*ld2+2.4+6./sp)/(1.+3./p4);
00377 G4double To=(.3*ld2+22.3+5./sp)/(1.+1./p4);
00378 sigma=(To-El)+.1/md;
00379 }
00380 else if(tZ==1 && tN==1)
00381 {
00382 G4double p2=P*P;
00383 G4double d=lP-2.7;
00384 G4double f=lP+1.25;
00385 G4double g_value=lP-.017;
00386 sigma=(.55*d*d+38.+23./std::sqrt(P))/(1.+.3/p2/p2)+18./(f*f+.1089)+.02/(g_value*g_value+.0025);
00387 }
00388 else if(tZ<97 && tN<152)
00389 {
00390 G4double d=lP-4.2;
00391 G4double p2=P*P;
00392 G4double p4=p2*p2;
00393 G4double a=tN+tZ;
00394 G4double al=std::log(a);
00395 G4double sa=std::sqrt(a);
00396 G4double ssa=std::sqrt(sa);
00397 G4double a2=a*a;
00398 G4double c=41.*std::exp(al*.68)*(1.+44./a2)/(1.+8./a)/(1.+200./a2/a2);
00399 G4double f=290.*ssa/(1.+34./a/ssa);
00400 G4double g_value=-1.32-al*.043;
00401 G4double u=lP-g_value;
00402 G4double h=al*(.4-.055*al);
00403 G4double r=.01+a2*5.E-8;
00404 sigma=(c+d*d)/(1.+(.2-.009*sa)/p4)+f/(u*u+h*h)/(1.+r/p2);
00405 #ifdef pdebug
00406 G4cout<<"G4QPiPlusNucCS::CSForm: A="<<a<<",P="<<P<<",CS="<<sigma<<",c="<<c<<",g="<<g_value
00407 <<",d="<<d<<",r="<<r<<",e="<<e<<",h="<<h<<G4endl;
00408 #endif
00409 }
00410 else
00411 {
00412 G4cerr<<"-Warning-G4QPiPlusNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
00413 sigma=0.;
00414 }
00415 if(sigma<0.) return 0.;
00416 return sigma;
00417 }