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 "G4QIonIonCrossSection.hh"
00048 #include "G4SystemOfUnits.hh"
00049
00050
00051 G4double* G4QIonIonCrossSection::lastLENI=0;
00052 G4double* G4QIonIonCrossSection::lastHENI=0;
00053 G4double* G4QIonIonCrossSection::lastLENE=0;
00054 G4double* G4QIonIonCrossSection::lastHENE=0;
00055 G4int G4QIonIonCrossSection::lastPDG=0;
00056 G4int G4QIonIonCrossSection::lastN=0;
00057 G4int G4QIonIonCrossSection::lastZ=0;
00058 G4double G4QIonIonCrossSection::lastP=0.;
00059 G4double G4QIonIonCrossSection::lastTH=0.;
00060 G4double G4QIonIonCrossSection::lastICS=0.;
00061 G4double G4QIonIonCrossSection::lastECS=0.;
00062 G4int G4QIonIonCrossSection::lastI=0;
00063
00064
00065 G4VQCrossSection* G4QIonIonCrossSection::GetPointer()
00066 {
00067 static G4QIonIonCrossSection theCrossSection;
00068 return &theCrossSection;
00069 }
00070
00071
00072
00073 G4double G4QIonIonCrossSection::GetCrossSection(G4bool fCS, G4double pMom, G4int tZ,
00074 G4int tN, G4int pPDG)
00075 {
00076 static G4int j;
00077 static std::vector <G4int> colPDG;
00078 static std::vector <G4int> colN;
00079 static std::vector <G4int> colZ;
00080 static std::vector <G4double> colP;
00081 static std::vector <G4double> colTH;
00082 static std::vector <G4double> colICS;
00083 static std::vector <G4double> colECS;
00084
00085 #ifdef pdebug
00086 G4cout<<"G4QIICS::GetCS:>>> f="<<fCS<<", Z="<<tZ<<"("<<lastZ<<"), N="<<tN<<"("<<lastN
00087 <<"),PDG="<<pPDG<<"("<<lastPDG<<"), p="<<pMom<<"("<<lastTH<<")"<<",Sz="
00088 <<colN.size()<<G4endl;
00089 #endif
00090 if(!pPDG)
00091 {
00092 #ifdef pdebug
00093 G4cout<<"G4QIonIonCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00094 #endif
00095 return 0.;
00096 }
00097 G4bool in=false;
00098 if(tN!=lastN || tZ!=lastZ || pPDG!=lastPDG)
00099 {
00100 in = false;
00101 lastP = 0.;
00102 lastPDG = pPDG;
00103 lastN = tN;
00104 lastZ = tZ;
00105 lastI = colN.size();
00106 j = 0;
00107 #ifdef pdebug
00108 G4cout<<"G4QIICS::GetCS:FindI="<<lastI<<",pPDG="<<pPDG<<",tN="<<tN<<",tZ="<<tZ<<G4endl;
00109 #endif
00110 if(lastI) for(G4int i=0; i<lastI; i++)
00111 {
00112 #ifdef pdebug
00113 G4cout<<"G4QII::GCS:P="<<colPDG[i]<<",N="<<colN[i]<<",Z="<<colZ[i]<<",j="<<j<<G4endl;
00114 #endif
00115 if(colPDG[i]==pPDG && colN[i]==tN && colZ[i]==tZ)
00116 {
00117 lastI=i;
00118 lastTH =colTH[i];
00119 #ifdef pdebug
00120 G4cout<<"G4QIICS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00121 #endif
00122 if(pMom<=lastTH)
00123 {
00124 #ifdef pdebug
00125 G4cout<<"G4QIICS::GetCS:Found P="<<pMom<<"<Threshold="<<lastTH<<"->XS=0"<<G4endl;
00126 #endif
00127 return 0.;
00128 }
00129 lastP =colP [i];
00130 lastICS=colICS[i];
00131 lastECS=colECS[i];
00132 if(std::fabs(lastP/pMom-1.)<tolerance)
00133 {
00134 #ifdef pdebug
00135 G4cout<<"G4QIonIonCS::GetCS:P="<<pMom<<",InXS="<<lastICS*millibarn<<",ElXS="
00136 <<lastECS*millibarn<<G4endl;
00137 #endif
00138 CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom);
00139 if(fCS) return lastICS*millibarn;
00140 return lastECS*millibarn;
00141 }
00142 in = true;
00143
00144 #ifdef pdebug
00145 G4cout<<"G4QIICS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00146 #endif
00147 lastICS=CalculateCrossSection( true,-1,j,lastPDG,lastZ,lastN,pMom);
00148 lastECS=CalculateCrossSection(false,-1,j,lastPDG,lastZ,lastN,pMom);
00149 #ifdef pdebug
00150 G4cout<<"G4QIonIonCS::GetCS:=>New(inDB) InCS="<<lastICS<<",ElCS="<<lastECS<<G4endl;
00151 #endif
00152 if((lastICS<=0. || lastECS<=0.) && pMom>lastTH)
00153 {
00154 #ifdef pdebug
00155 G4cout<<"G4QIonIonCS::GetCS:New,T="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00156 #endif
00157 lastTH=pMom;
00158 }
00159 break;
00160 }
00161 #ifdef pdebug
00162 G4cout<<"--->G4QIonIonCrossSec::GetCrosSec: pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
00163 <<",Z["<<i<<"]="<<colZ[i]<<",PDG="<<colPDG[i]<<G4endl;
00164 #endif
00165 j++;
00166 }
00167 if(!in)
00168 {
00169 #ifdef pdebug
00170 G4cout<<"G4QIICS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
00171 #endif
00173 lastICS=CalculateCrossSection(true ,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
00174 lastECS=CalculateCrossSection(false,0,j,lastPDG,lastZ,lastN,pMom);
00175 if(lastICS<=0. || lastECS<=0.)
00176 {
00177 lastTH = ThresholdEnergy(tZ, tN);
00178 #ifdef pdebug
00179 G4cout<<"G4QIonIonCrossSect::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
00180 #endif
00181 if(pMom>lastTH)
00182 {
00183 #ifdef pdebug
00184 G4cout<<"G4QIonIonCS::GetCS:1-st,P="<<pMom<<">Thresh="<<lastTH<<"->XS=0"<<G4endl;
00185 #endif
00186 lastTH=pMom;
00187 }
00188 }
00189 #ifdef pdebug
00190 G4cout<<"G4QIICS::GetCS: *New* ICS="<<lastICS<<", ECS="<<lastECS<<",N="<<lastN<<",Z="
00191 <<lastZ<<G4endl;
00192 #endif
00193 colN.push_back(tN);
00194 colZ.push_back(tZ);
00195 colPDG.push_back(pPDG);
00196 colP.push_back(pMom);
00197 colTH.push_back(lastTH);
00198 colICS.push_back(lastICS);
00199 colECS.push_back(lastECS);
00200 #ifdef pdebug
00201 G4cout<<"G4QIICS::GetCS:*1st*, P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn
00202 <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
00203 #endif
00204 if(fCS) return lastICS*millibarn;
00205 return lastECS*millibarn;
00206 }
00207 else
00208 {
00209 #ifdef pdebug
00210 G4cout<<"G4QIICS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00211 #endif
00212 colP[lastI]=pMom;
00213 colPDG[lastI]=pPDG;
00214 colICS[lastI]=lastICS;
00215 colECS[lastI]=lastECS;
00216 }
00217 }
00218 else if(pMom<=lastTH)
00219 {
00220 #ifdef pdebug
00221 G4cout<<"G4QIICS::GetCS: Current T="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00222 #endif
00223 return 0.;
00224 }
00225 else if(std::fabs(lastP/pMom-1.)<tolerance)
00226 {
00227 #ifdef pdebug
00228 G4cout<<"G4QIICS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", InCS="<<lastICS*millibarn
00229 <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
00230 #endif
00231 if(fCS) return lastICS*millibarn;
00232 return lastECS*millibarn;
00233 }
00234 else
00235 {
00236 #ifdef pdebug
00237 G4cout<<"G4QIICS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00238 #endif
00239 lastICS=CalculateCrossSection( true,1,j,lastPDG,lastZ,lastN,pMom);
00240 lastECS=CalculateCrossSection(false,1,j,lastPDG,lastZ,lastN,pMom);
00241 lastP=pMom;
00242 }
00243 #ifdef pdebug
00244 G4cout<<"G4QIICS::GetCroSec:*End*,P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn<<", ElCS="
00245 <<lastECS*millibarn<<"(mb)"<<G4endl;
00246 #endif
00247 if(fCS) return lastICS*millibarn;
00248 return lastECS*millibarn;
00249 }
00250
00251
00252 G4double G4QIonIonCrossSection::CalculateCrossSection(G4bool XS,G4int F,G4int I,G4int pPDG,
00253 G4int tZ,G4int tN, G4double TotMom)
00254 {
00255
00256
00257
00258 static const G4double THmin=0.;
00259 static const G4double dP=10.;
00260 static const G4int nL=100;
00261 static const G4double Pmin=THmin+(nL-1)*dP;
00262 static const G4double Pmax=300000.;
00263 static const G4int nH=100;
00264 static const G4double milP=std::log(Pmin);
00265 static const G4double malP=std::log(Pmax);
00266 static const G4double dlP=(malP-milP)/(nH-1);
00267
00268
00269 static std::vector <G4double*> LENI;
00270 static std::vector <G4double*> HENI;
00271 static std::vector <G4double*> LENE;
00272 static std::vector <G4double*> HENE;
00273 #ifdef debug
00274 G4cout<<"G4QIonIonCrossSection::CalcCS: Z="<<tZ<<", N="<<tN<<", P="<<TotMom<<G4endl;
00275 #endif
00276 G4int dPDG=pPDG/10;
00277 G4int zPDG=dPDG/1000;
00278 G4int zA=dPDG%1000;
00279 G4int pZ=zPDG%1000;
00280 G4int pN=zA-pZ;
00281 G4double Momentum=TotMom/zA;
00282 if (Momentum<THmin) return 0.;
00283 G4double sigma=0.;
00284 if(F&&I) sigma=0.;
00285 G4double tA=tN+tZ;
00286 G4double pA=zA;
00287 if(F<=0)
00288 {
00289 if(F<0 || !XS)
00290 {
00291 lastLENI=LENI[I];
00292 lastHENI=HENI[I];
00293 lastLENE=LENE[I];
00294 lastHENE=HENE[I];
00295 }
00296 else
00297 {
00298 lastLENI = new G4double[nL];
00299 lastHENI = new G4double[nH];
00300 lastLENE = new G4double[nL];
00301 lastHENE = new G4double[nH];
00302 G4int er=GetFunctions(pZ,pN,tZ,tN,lastLENI,lastHENI,lastLENE,lastHENE);
00303 if(er<1) G4cerr<<"*W*G4QIonIonCroSec::CalcCrossSection: pA="<<tA<<",tA="<<tA<<G4endl;
00304 #ifdef debug
00305 G4cout<<"G4QIonIonCrossSection::CalcCS: GetFunctions er="<<er<<",pA="<<pA<<",tA="<<tA
00306 <<G4endl;
00307 #endif
00308
00309 G4int sync=LENI.size();
00310 if(sync!=I) G4cout<<"*W*G4IonIonCrossSec::CalcCrossSect:Sync="<<sync<<"#"<<I<<G4endl;
00311 LENI.push_back(lastLENI);
00312 HENI.push_back(lastHENI);
00313 LENE.push_back(lastLENE);
00314 HENE.push_back(lastHENE);
00315 }
00316 }
00317
00318 if (Momentum<lastTH) return 0.;
00319 else if (Momentum<Pmin)
00320 {
00321 #ifdef debug
00322 G4cout<<"G4QIICS::CalCS:p="<<pA<<",t="<<tA<<",n="<<nL<<",T="<<THmin<<",d="<<dP<<G4endl;
00323 #endif
00324 if(tA<1. || pA<1.)
00325 {
00326 G4cout<<"-Warning-G4QIICS::CalcCS: pA="<<pA<<" or tA="<<tA<<" aren't nuclei"<<G4endl;
00327 sigma=0.;
00328 }
00329 else
00330 {
00331 G4double dPp=dP*pA;
00332 if(XS) sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENI);
00333 else sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENE);
00334 }
00335 #ifdef debugn
00336 if(sigma<0.) G4cout<<"-Warning-G4QIICS::CalcCS:pA="<<pA<<",tA="<<tA<<",XS="<<XS<<",P="
00337 <<Momentum<<", Th="<<THmin<<", dP="<<dP<<G4endl;
00338 #endif
00339 }
00340 else if (Momentum<Pmax*pA)
00341 {
00342 G4double lP=std::log(Momentum);
00343 #ifdef debug
00344 G4cout<<"G4QIonIonCS::CalcCS:before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl;
00345 #endif
00346 if(tA<1. || pA<1.)
00347 {
00348 G4cout<<"-Warning-G4QIICS::CalCS:pA="<<pA<<" or tA="<<tA<<" aren't composit"<<G4endl;
00349 sigma=0.;
00350 }
00351 else
00352 {
00353 G4double milPp=milP+std::log(pA);
00354 if(XS) sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENI);
00355 else sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENE);
00356 }
00357 }
00358 else
00359 {
00360 std::pair<G4double, G4double> inelel = CalculateXS(pZ, pN, tZ, tN, Momentum);
00361 if(XS) sigma=inelel.first;
00362 else sigma=inelel.second;
00363 }
00364 #ifdef debug
00365 G4cout<<"G4IonIonCrossSection::CalculateCrossSection: sigma="<<sigma<<G4endl;
00366 #endif
00367 if(sigma<0.) return 0.;
00368 return sigma;
00369 }
00370
00371
00372
00373
00374 G4int G4QIonIonCrossSection::GetFunctions(G4int pZ,G4int pN,G4int tZ,G4int tN,G4double* li,
00375 G4double* hi, G4double* le, G4double* he)
00376 {
00377
00378 static const G4double THmin=0.;
00379 static const G4double dP=10.;
00380 static const G4int nL=100;
00381 static const G4double Pmin=THmin+(nL-1)*dP;
00382 static const G4double Pmax=300000.;
00383 static const G4int nH=100;
00384 static const G4double milP=std::log(Pmin);
00385 static const G4double malP=std::log(Pmax);
00386 static const G4double dlP=(malP-milP)/(nH-1);
00387 static const G4double lP=std::exp(dlP);
00388
00389 if(pZ<1 || pN<0 || tZ<1 || tN<0)
00390 {
00391 G4cout<<"-W-G4QIonIonCS::GetFunct:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl;
00392 return -1;
00393 }
00394 G4int pA=pN+pZ;
00395 G4double dPp=dP*pA;
00396 G4double Mom=THmin;
00397 for(G4int k=0; k<nL; k++)
00398 {
00399 std::pair<G4double,G4double> len = CalculateXS(pZ, pN, tZ, tN, Mom);
00400 li[k]=len.first;
00401 le[k]=len.second;
00402 Mom+=dPp;
00403 }
00404 G4double lMom=Pmin*pA;
00405 for(G4int j=0; j<nH; j++)
00406 {
00407 std::pair<G4double,G4double> len = CalculateXS(pZ, pN, pZ, pN, lMom);
00408 hi[j]=len.first;
00409 he[j]=len.second;
00410 lMom*=lP;
00411 }
00412 #ifdef debug
00413 G4cout<<"G4QIonIonCS::GetFunctions: pZ="<<pZ<<", tZ="<<tZ<<" pair is calculated"<<G4endl;
00414 #endif
00415 return 1;
00416 }
00417
00418
00419 std::pair<G4double,G4double> G4QIonIonCrossSection::CalculateXS(G4int pZ,G4int pN,G4int tZ,
00420 G4int tN, G4double Mom)
00421 {
00422 static G4VQCrossSection* PElCSman = G4QProtonElasticCrossSection::GetPointer();
00423 static G4VQCrossSection* NElCSman = G4QNeutronElasticCrossSection::GetPointer();
00424 static G4VQCrossSection* InelPCSman = G4QProtonNuclearCrossSection::GetPointer();
00425 static G4VQCrossSection* InelNCSman = G4QNeutronNuclearCrossSection::GetPointer();
00426 G4double pA=pZ+pN;
00427 G4double tA=tZ+tN;
00428 if(pA<.9 || tA<.9 ||pA>239. || tA>239 || Mom < 0.) return std::make_pair(0.,0.);
00429 G4double inCS=0.;
00430 G4double elCS=0.;
00431 if(pA<1.1 )
00432 {
00433 if (pZ == 1 && !pN)
00434 {
00435 inCS=InelPCSman->GetCrossSection(true, Mom, tZ, tN, 2212);
00436 elCS=PElCSman->GetCrossSection(true, Mom, tZ, tN, 2212);
00437 }
00438 else if(pN == 1 && !pZ)
00439 {
00440 inCS=InelNCSman->GetCrossSection(true, Mom, tZ, tN, 2112);
00441 elCS=NElCSman->GetCrossSection(true, Mom, tZ, tN, 2112);
00442 }
00443 else G4cerr<<"-Warn-G4QIICS::CaCS:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl;
00444 }
00445 else
00446 {
00447 G4double T=ThresholdMomentum(pZ, pN, tZ, tN);
00448 if(Mom<=T)
00449 {
00450 elCS=0.;
00451 inCS=0.;
00452 }
00453 else
00454 {
00455 G4double P2=Mom*Mom;
00456 G4double T2=T*T;
00457 G4double R=1.-T2/P2;
00458
00459
00460
00461
00462 G4double tot=R*CalculateTotal(pA, tA, Mom);
00463 G4double rat=CalculateElTot(pA, tA, Mom);
00464 elCS=tot*rat;
00465 inCS=tot-elCS;
00466 }
00467 }
00468 return std::make_pair(inCS,elCS);
00469 }
00470
00471
00472 G4double G4QIonIonCrossSection::CalculateTotal(G4double pA, G4double tA, G4double Mom)
00473 {
00474 G4double y=std::log(Mom/1000.);
00475 G4double ab=pA+tA;
00476 G4double al=std::log(ab);
00477 G4double ap=std::log(pA*tA);
00478 G4double e=std::pow(pA,0.1)+std::pow(tA,0.1);
00479 G4double d=e-1.55/std::pow(al,0.2);
00480 G4double f=4.;
00481 if(pA>4. && tA>4.) f=3.3+130./ab/ab+2.25/e;
00482 G4double c=(385.+11.*ab/(1.+.02*ab*al)+(.5*ab+500./al/al)/al)*std::pow(d,f);
00483 G4double r=y-3.-4./ap/ap;
00484 #ifdef pdebug
00485 G4cout<<"G4QIonIonCS::CalcTot:P="<<Mom<<", stot="<<c+d*r*r<<", d="<<d<<", r="<<r<<G4endl;
00486 #endif
00487 return c+d*r*r;
00488 }
00489
00490
00491 G4double G4QIonIonCrossSection::CalculateElTot(G4double pA, G4double tA, G4double Mom)
00492 {
00493 G4double y=std::log(Mom/1000.);
00494 G4double ab=pA*tA;
00495 G4double ap=std::log(ab);
00496 G4double r=y-3.92-1.73/ap/ap;
00497 G4double d=.00166/(1.+.002*std::pow(ab,1.33333));
00498 G4double al=std::log(pA+tA);
00499 G4double e=1.+.235*(std::fabs(ap-5.78));
00500 if (std::fabs(pA-2.)<.1 && std::fabs(tA-2.)<.1) e=2.18;
00501 else if(std::fabs(pA-2.)<.1 && std::fabs(tA-3.)<.1) e=1.90;
00502 else if(std::fabs(pA-3.)<.1 && std::fabs(tA-2.)<.1) e=1.90;
00503 else if(std::fabs(pA-2.)<.1 && std::fabs(tA-4.)<.1) e=1.65;
00504 else if(std::fabs(pA-4.)<.1 && std::fabs(tA-2.)<.1) e=1.65;
00505 else if(std::fabs(pA-3.)<.1 && std::fabs(tA-4.)<.1) e=1.32;
00506 else if(std::fabs(pA-4.)<.1 && std::fabs(tA-3.)<.1) e=1.32;
00507 else if(std::fabs(pA-4.)<.1 && std::fabs(tA-4.)<.1) e=1.;
00508 G4double f=.37+.0188*al;
00509 G4double g_value=std::log(std::pow(pA,0.35)+std::pow(tA,0.35));
00510 G4double h=g_value*g_value;
00511 G4double c=f/(1.+e/h/h);
00512 #ifdef pdebug
00513 G4cout<<"G4QIonIonCS::CalcElT:P="<<Mom<<",el/tot="<<c+d*r*r<<",d="<<d<<", r="<<r<<G4endl;
00514 #endif
00515 return c+d*r*r;
00516 }
00517
00518
00519 G4double G4QIonIonCrossSection::ThresholdMomentum(G4int pZ, G4int pN, G4int tZ, G4int tN)
00520 {
00521 static const G4double third=1./3.;
00522 static const G4double pM = G4QPDGCode(2212).GetMass();
00523 static const G4double tpM= pM+pM;
00524 if(pZ<.99 || pN<0. || tZ<.99 || tN<0.) return 0.;
00525 G4double tA=tZ+tN;
00526 G4double pA=pZ+pN;
00527
00528 G4double dE=pZ*tZ/(std::pow(pA,third)+std::pow(tA,third))/pA;
00529 return std::sqrt(dE*(tpM+dE));
00530 }