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
00048
00049 #include "G4QANuENuclearCrossSection.hh"
00050 #include "G4SystemOfUnits.hh"
00051
00052
00053 G4bool G4QANuENuclearCrossSection::onlyCS=true;
00054 G4double G4QANuENuclearCrossSection::lastSig=0.;
00055 G4double G4QANuENuclearCrossSection::lastQEL=0.;
00056 G4int G4QANuENuclearCrossSection::lastL=0;
00057 G4double G4QANuENuclearCrossSection::lastE=0.;
00058 G4double* G4QANuENuclearCrossSection::lastEN=0;
00059 G4double* G4QANuENuclearCrossSection::lastTX=0;
00060 G4double* G4QANuENuclearCrossSection::lastQE=0;
00061 G4int G4QANuENuclearCrossSection::lastPDG=0;
00062 G4int G4QANuENuclearCrossSection::lastN=0;
00063 G4int G4QANuENuclearCrossSection::lastZ=0;
00064 G4double G4QANuENuclearCrossSection::lastP=0.;
00065 G4double G4QANuENuclearCrossSection::lastTH=0.;
00066 G4double G4QANuENuclearCrossSection::lastCS=0.;
00067 G4int G4QANuENuclearCrossSection::lastI=0;
00068 std::vector<G4double*>* G4QANuENuclearCrossSection::TX = new std::vector<G4double*>;
00069 std::vector<G4double*>* G4QANuENuclearCrossSection::QE = new std::vector<G4double*>;
00070
00071
00072 G4VQCrossSection* G4QANuENuclearCrossSection::GetPointer()
00073 {
00074 static G4QANuENuclearCrossSection theCrossSection;
00075 return &theCrossSection;
00076 }
00077
00078 G4QANuENuclearCrossSection::~G4QANuENuclearCrossSection()
00079 {
00080 G4int lens=TX->size();
00081 for(G4int i=0; i<lens; ++i) delete[] (*TX)[i];
00082 delete TX;
00083 G4int hens=QE->size();
00084 for(G4int i=0; i<hens; ++i) delete[] (*QE)[i];
00085 delete QE;
00086 }
00087
00088
00089
00090 G4double G4QANuENuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
00091 G4int tgZ, G4int tgN, G4int pPDG)
00092 {
00093 static G4int j;
00094 static std::vector <G4int> colPDG;
00095 static std::vector <G4int> colN;
00096 static std::vector <G4int> colZ;
00097 static std::vector <G4double> colP;
00098 static std::vector <G4double> colTH;
00099 static std::vector <G4double> colCS;
00100
00101 G4double pEn=pMom;
00102 #ifdef debug
00103 G4cout<<"G4QAENCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
00104 <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
00105 <<colN.size()<<G4endl;
00106
00107 #endif
00108 if(pPDG!=-12)
00109 {
00110 #ifdef debug
00111 G4cout<<"G4QAENCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00112
00113 #endif
00114 return 0.;
00115 }
00116 G4bool in=false;
00117 if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)
00118 {
00119 in = false;
00120 lastP = 0.;
00121 lastPDG = pPDG;
00122 lastN = tgN;
00123 lastZ = tgZ;
00124 lastI = colN.size();
00125 j = 0;
00126 if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG)
00127 {
00128 if(colN[i]==tgN && colZ[i]==tgZ)
00129 {
00130 lastI=i;
00131 lastTH =colTH[i];
00132 #ifdef pdebug
00133 G4cout<<"G4QAENCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00134
00135 #endif
00136 if(pEn<=lastTH)
00137 {
00138 #ifdef pdebug
00139 G4cout<<"G4QAENCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl;
00140
00141 #endif
00142 return 0.;
00143 }
00144 lastP =colP [i];
00145 lastCS =colCS[i];
00146 if(std::fabs(lastP/pMom-1.)<tolerance)
00147 {
00148 #ifdef pdebug
00149 G4cout<<"G4QAENCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00150
00151 #endif
00152 return lastCS*millibarn;
00153 }
00154 in = true;
00155
00156 #ifdef pdebug
00157 G4cout<<"G4QAENCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00158 #endif
00159 lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom);
00160 #ifdef pdebug
00161 G4cout<<"G4QAENCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00162
00163 #endif
00164 if(lastCS<=0. && pEn>lastTH)
00165 {
00166 #ifdef pdebug
00167 G4cout<<"G4QAENCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00168 #endif
00169 lastTH=pEn;
00170 }
00171 break;
00172 }
00173 #ifdef pdebug
00174 G4cout<<"---G4QAENCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
00175 <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl;
00176
00177 #endif
00178 j++;
00179 }
00180 if(!in)
00181 {
00182 #ifdef pdebug
00183 G4cout<<"G4QAENCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<lastI<<G4endl;
00184 #endif
00186 lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
00187 if(lastCS<=0.)
00188 {
00189 lastTH = ThresholdEnergy(tgZ, tgN);
00190 #ifdef pdebug
00191 G4cout<<"G4QAENCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
00192 #endif
00193 if(pEn>lastTH)
00194 {
00195 #ifdef pdebug
00196 G4cout<<"G4QAENCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00197 #endif
00198 lastTH=pEn;
00199 }
00200 }
00201 #ifdef pdebug
00202 G4cout<<"G4QAENCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00203
00204 #endif
00205 colN.push_back(tgN);
00206 colZ.push_back(tgZ);
00207 colPDG.push_back(pPDG);
00208 colP.push_back(pMom);
00209 colTH.push_back(lastTH);
00210 colCS.push_back(lastCS);
00211 #ifdef pdebug
00212 G4cout<<"G4QAENCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl;
00213
00214 #endif
00215 return lastCS*millibarn;
00216 }
00217 else
00218 {
00219 #ifdef pdebug
00220 G4cout<<"G4QAENCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00221 #endif
00222 colP[lastI]=pMom;
00223 colPDG[lastI]=pPDG;
00224 colCS[lastI]=lastCS;
00225 }
00226 }
00227 else if(pEn<=lastTH)
00228 {
00229 #ifdef pdebug
00230 G4cout<<"G4QAENCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00231
00232 #endif
00233 return 0.;
00234 }
00235 else if(std::fabs(lastP/pMom-1.)<tolerance)
00236 {
00237 #ifdef pdebug
00238 G4cout<<"G4QAENCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00239
00240 #endif
00241 return lastCS*millibarn;
00242 }
00243 else
00244 {
00245 #ifdef pdebug
00246 G4cout<<"G4QAENCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00247 #endif
00248 lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom);
00249 lastP=pMom;
00250 }
00251 #ifdef pdebug
00252 G4cout<<"G4QAENCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00253
00254 #endif
00255 return lastCS*millibarn;
00256 }
00257
00258
00259 G4double G4QANuENuclearCrossSection::ThresholdEnergy(G4int Z, G4int N, G4int)
00260 {
00261
00262
00263
00264 static const G4double mN=.931494043;
00265 static const G4double dmN=mN+mN;
00266 static const G4double me=.00051099892;
00267 static const G4double me2=me*me;
00268 static const G4double thresh=me+me2/dmN;
00269
00270
00271 G4double dN=0.;
00272 if(Z>0||N>0) dN=thresh*GeV;
00273
00274 return dN;
00275 }
00276
00277
00278 G4double G4QANuENuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I,
00279 G4int , G4int targZ, G4int targN, G4double Momentum)
00280 {
00281 static const G4double mb38=1.E-11;
00282 static const G4int nE=65;
00283 static const G4int mL=nE-1;
00284 static const G4double mN=.931494043;
00285 static const G4double dmN=mN+mN;
00286 static const G4double me=.00051099892;
00287 static const G4double me2=me*me;
00288 static const G4double EMi=me+me2/dmN;
00289 static const G4double EMa=300.;
00290
00291 static std::vector <G4double> colH;
00292 static G4bool first=true;
00293
00294
00295
00296 onlyCS=CS;
00297 lastE=Momentum/GeV;
00298 if (lastE<=EMi)
00299 {
00300 lastE=0.;
00301 lastSig=0.;
00302 return 0.;
00303 }
00304 G4int Z=targZ;
00305 if(F<=0)
00306 {
00307 if(F<0)
00308 {
00309 lastTX =(*TX)[I];
00310 lastQE =(*QE)[I];
00311 }
00312 else
00313 {
00314 if(first)
00315 {
00316 lastEN = new G4double[nE];
00317 Z=-Z;
00318 first=false;
00319 }
00320 lastTX = new G4double[nE];
00321 lastQE = new G4double[nE];
00322 G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);
00323 if(res<0) G4cerr<<"*W*G4NuENuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
00324
00325 G4int sync=TX->size();
00326 if(sync!=I) G4cerr<<"***G4NuENuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
00327 TX->push_back(lastTX);
00328 QE->push_back(lastQE);
00329 }
00330 }
00331
00332 if (lastE<=EMi)
00333 {
00334 lastE=0.;
00335 lastSig=0.;
00336 return 0.;
00337 }
00338 if(lastE<EMa)
00339 {
00340 G4int chk=1;
00341 G4int ran=mL/2;
00342 G4int sep=ran;
00343 while(ran>=2)
00344 {
00345 G4int newran=ran/2;
00346 G4double oldL=lastEN[sep];
00347 if(lastE<=oldL) sep-=newran;
00348 else sep+=newran;
00349 #ifdef pdebug
00350 G4cout<<"G4ANuE::CCS:n="<<newran<<",s="<<sep<<",E="<<lastE<<",oE="<<oldL<<G4endl;
00351 #endif
00352 ran=newran;
00353 chk=chk+chk;
00354 }
00355 if(chk+chk!=mL) G4cerr<<"*Warn*G4NuENuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
00356 G4double lowE=lastEN[sep];
00357 G4double highE=lastEN[sep+1];
00358 G4double lowTX=lastTX[sep];
00359 if(lastE<lowE||sep>=mL||lastE>highE)
00360 G4cerr<<"*Warn*G4ANuENuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
00361 <<", sep="<<sep<<", mL="<<mL<<G4endl;
00362 lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX;
00363 if(!onlyCS)
00364 {
00365 G4double lowQE=lastQE[sep];
00366 lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
00367 #ifdef pdebug
00368 G4cout<<"G4ANuENuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
00369 #endif
00370 }
00371 }
00372 else
00373 {
00374 lastSig=lastTX[mL];
00375 lastQEL=lastQE[mL];
00376 }
00377 if(lastQEL<0.) lastQEL = 0.;
00378 if(lastSig<0.) lastSig = 0.;
00379
00380 lastSig*=mb38;
00381 if(!onlyCS) lastQEL*=mb38;
00382 return lastSig;
00383 }
00384
00385
00386
00387
00388
00389
00390 G4int G4QANuENuclearCrossSection::GetFunctions (G4int z, G4int n,
00391 G4double* t, G4double* q, G4double* e)
00392 {
00393 static const G4double mN=.931494043;
00394 static const G4double dmN=mN+mN;
00395 static const G4double me=.00051099892;
00396 static const G4double me2=me*me;
00397 static const G4double thresh=me+me2/dmN;
00398 static const G4int nE=65;
00399 static const G4double nuEn[nE]={thresh,
00400 .00051331,.00053602,.00056078,.00058783,.00061743,.00064990,.00068559,.00072492,
00401 .00076834,.00081641,.00086975,.00092912,.00099536,.00106950,.00115273,.00124646,
00402 .00135235,.00147241,.00160901,.00176503,.00194392,.00214986,.00238797,.00266448,
00403 .00298709,.00336531,.00381094,.00433879,.00496745,.00572047,.00662785,.00772806,
00404 .00907075,.01072050,.01276190,.01530660,.01850330,.02255110,.02771990,.03437780,
00405 .04303240,.05438970,.06944210,.08959920,.11688400,.15423600,.20597200,.27851200,
00406 .38153100,.52979600,.74616300,1.0665200,1.5480900,2.2834800,3.4251100,5.2281000,
00407 8.1270200,12.875900,20.808500,34.331200,57.877800,99.796200,176.16300,318.68200};
00408 static const G4double TOTX[nE]={0.,
00409 .00046923,.00160693,.00229560,.00288772,.00344344,.00398831,.00453725,.00510087,
00410 .00568797,.00630663,.00696489,.00767121,.00843477,.00926586,.01017610,.01117910,
00411 .01229040,.01352820,.01491430,.01647420,.01823820,.02024280,.02253130,.02515600,
00412 .02818010,.03167950,.03574640,.04049260,.04605330,.05739330,.07028190,.08396290,
00413 .09969400,.11756200,.13812500,.16207000,.19099600,.22447700,.26434900,.31051900,
00414 .36357400,.42242900,.48500000,.54692800,.60140500,.63996800,.65519300,.64339200,
00415 .60784300,.55734400,.50212600,.45116800,.40962600,.37854800,.35702100,.34330900,
00416 .33576200,.33300100,.33387100,.33737700,.34235200,.34880100,.35507500,.35961200};
00417 static const G4double QELX[nE]={0.,
00418 2.40856e-7,8.61338e-7,1.28732e-6,1.69747e-6,2.12608e-6,2.59201e-6,3.11071e-6,3.69770e-6,
00419 4.37028e-6,5.14877e-6,6.05773e-6,7.12746e-6,8.39567e-6,9.90987e-6,1.17304e-5,1.39343e-5,
00420 1.66209e-5,1.99191e-5,2.39974e-5,2.90775e-5,3.54536e-5,4.35191e-5,5.38039e-5,6.70278e-5,
00421 8.41765e-5,.000106611,.000136228,.000175689,.000228768,.000328317,.000465818,.000648870,
00422 .000904299,.001260330,.001762740,.002480740,.003534040,.005062210,.007327720,.010675000,
00423 .015645500,.022975800,.033679400,.049004300,.070294800,.098706100,.134951000,.179193000,
00424 .231297000,.287287000,.344760000,.404397000,.465915000,.527060000,.584085000,.632748000,
00425 .671346000,.699744000,.719355000,.732541000,.740996000,.746249000,.749265000,.751160000};
00426
00427 G4int first=0;
00428 if(z<0.)
00429 {
00430 first=1;
00431 z=-z;
00432 }
00433 if(z<1 || z>92)
00434 {
00435 G4cout<<"***G4QANuENuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
00436 return -1;
00437 }
00438 for(G4int k=0; k<nE; k++)
00439 {
00440 G4double a=n+z;
00441 G4double za=z+a;
00442 G4double dz=z+z;
00443 G4double da=a+a;
00444 G4double ta=da+a;
00445 if(first) e[k]=nuEn[k];
00446 t[k]=TOTX[k]*nuEn[k]*(za+za)/ta+QELX[k]*(dz+dz-da)/ta;
00447 q[k]=QELX[k]*dz/a;
00448 }
00449 return first;
00450 }
00451
00452
00453 G4double G4QANuENuclearCrossSection::GetQEL_ExchangeQ2()
00454 {
00455 static const G4double me=.00051099892;
00456 static const G4double me2=me*me;
00457 static const G4double hme2=me2/2;
00458 static const double MN=.931494043;
00459 static const double MN2=MN*MN;
00460 static const G4double power=-3.5;
00461 static const G4double pconv=1./power;
00462 static const G4int nQ2=101;
00463 static const G4int lQ2=nQ2-1;
00464 static const G4int bQ2=lQ2-1;
00465
00466 static const G4double Xl[nQ2]={5.20224e-16,
00467 .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898,
00468 .093878, .103332, .112844, .122410, .132023, .141680, .151376, .161109, .170875, .180672,
00469 .190499, .200352, .210230, .220131, .230055, .239999, .249963, .259945, .269944, .279960,
00470 .289992, .300039, .310099, .320173, .330260, .340359, .350470, .360592, .370724, .380867,
00471 .391019, .401181, .411352, .421531, .431719, .441915, .452118, .462329, .472547, .482771,
00472 .493003, .503240, .513484, .523734, .533989, .544250, .554517, .564788, .575065, .585346,
00473 .595632, .605923, .616218, .626517, .636820, .647127, .657438, .667753, .678072, .688394,
00474 .698719, .709048, .719380, .729715, .740053, .750394, .760738, .771085, .781434, .791786,
00475 .802140, .812497, .822857, .833219, .843582, .853949, .864317, .874687, .885060, .895434,
00476 .905810, .916188, .926568, .936950, .947333, .957719, .968105, .978493, .988883, .999275};
00477
00478 static const G4double Xmax=Xl[lQ2];
00479 static const G4double Xmin=Xl[0];
00480 static const G4double dX=(Xmax-Xmin)/lQ2;
00481 static const G4double inl[nQ2]={0,
00482 1.52225, 2.77846, 3.96651, 5.11612, 6.23990, 7.34467, 8.43466, 9.51272, 10.5809, 11.6406,
00483 12.6932, 13.7394, 14.7801, 15.8158, 16.8471, 17.8743, 18.8979, 19.9181, 20.9353, 21.9496,
00484 22.9614, 23.9707, 24.9777, 25.9826, 26.9855, 27.9866, 28.9860, 29.9837, 30.9798, 31.9745,
00485 32.9678, 33.9598, 34.9505, 35.9400, 36.9284, 37.9158, 38.9021, 39.8874, 40.8718, 41.8553,
00486 42.8379, 43.8197, 44.8007, 45.7810, 46.7605, 47.7393, 48.7174, 49.6950, 50.6718, 51.6481,
00487 52.6238, 53.5990, 54.5736, 55.5476, 56.5212, 57.4943, 58.4670, 59.4391, 60.4109, 61.3822,
00488 62.3531, 63.3236, 64.2937, 65.2635, 66.2329, 67.2019, 68.1707, 69.1390, 70.1071, 71.0748,
00489 72.0423, 73.0095, 73.9763, 74.9429, 75.9093, 76.8754, 77.8412, 78.8068, 79.7721, 80.7373,
00490 81.7022, 82.6668, 83.6313, 84.5956, 85.5596, 86.5235, 87.4872, 88.4507, 89.4140, 90.3771,
00491 91.3401, 92.3029, 93.2656, 94.2281, 95.1904, 96.1526, 97.1147, 98.0766, 99.0384, 100.000};
00492 G4double Enu=lastE;
00493 G4double dEnu=Enu+Enu;
00494 G4double Enu2=Enu*Enu;
00495 G4double ME=Enu*MN;
00496 G4double dME=ME+ME;
00497 G4double dEMN=(dEnu+MN)*ME;
00498 G4double MEm=ME-hme2;
00499 G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
00500 G4double E2M=MN*Enu2-(Enu+MN)*hme2;
00501 G4double ymax=(E2M+sqE)/dEMN;
00502 G4double ymin=(E2M-sqE)/dEMN;
00503 G4double rmin=1.-ymin;
00504 G4double rhm2E=hme2/Enu2;
00505 G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
00506 G4double Q2ma=dME*ymax;
00507 G4double Xma=std::pow((1.+Q2mi),power);
00508 G4double Xmi=std::pow((1.+Q2ma),power);
00509
00510 G4double rXi=(Xmi-Xmin)/dX;
00511 G4int iXi=static_cast<int>(rXi);
00512 if(iXi<0) iXi=0;
00513 if(iXi>bQ2) iXi=bQ2;
00514 G4double dXi=rXi-iXi;
00515 G4double bnti=inl[iXi];
00516 G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
00517
00518 G4double rXa=(Xma-Xmin)/dX;
00519 G4int iXa=static_cast<int>(rXa);
00520 if(iXa<0) iXa=0;
00521 if(iXa>bQ2) iXa=bQ2;
00522 G4double dXa=rXa-iXa;
00523 G4double bnta=inl[iXa];
00524 G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
00525
00526 G4double intx=inti+(inta-inti)*G4UniformRand();
00527 G4int intc=static_cast<int>(intx);
00528 if(intc<0) intc=0;
00529 if(intc>bQ2) intc=bQ2;
00530 G4double dint=intx-intc;
00531 G4double mX=Xl[intc];
00532 G4double X=mX+dint*(Xl[intc+1]-mX);
00533 G4double Q2=std::pow(X,pconv)-1.;
00534 return Q2*GeV*GeV;
00535 }
00536
00537
00538 G4double G4QANuENuclearCrossSection::GetNQE_ExchangeQ2()
00539 {
00540 static const double mpi=.13957018;
00541 static const G4double me=.00051099892;
00542 static const G4double me2=me*me;
00543 static const G4double hme2=me2/2;
00544 static const double MN=.931494043;
00545 static const double MN2=MN*MN;
00546 static const double dMN=MN+MN;
00547 static const double mcV=(dMN+mpi)*mpi;
00548 static const G4int power=7;
00549 static const G4double pconv=1./power;
00550 static const G4int nX=21;
00551 static const G4int lX=nX-1;
00552 static const G4int bX=lX-1;
00553 static const G4int nE=20;
00554 static const G4int bE=nE-1;
00555 static const G4int pE=bE-1;
00556
00557 static const G4double X0[nX]={5.21412e-05,
00558 .437860, .681908, .891529, 1.08434, 1.26751, 1.44494, 1.61915, 1.79198, 1.96493, 2.13937,
00559 2.31664, 2.49816, 2.68559, 2.88097, 3.08705, 3.30774, 3.54917, 3.82233, 4.15131, 4.62182};
00560 static const G4double X1[nX]={.00102591,
00561 1.00443, 1.55828, 2.03126, 2.46406, 2.87311, 3.26723, 3.65199, 4.03134, 4.40835, 4.78561,
00562 5.16549, 5.55031, 5.94252, 6.34484, 6.76049, 7.19349, 7.64917, 8.13502, 8.66246, 9.25086};
00563 static const G4double X2[nX]={.0120304,
00564 2.59903, 3.98637, 5.15131, 6.20159, 7.18024, 8.10986, 9.00426, 9.87265, 10.7217, 11.5564,
00565 12.3808, 13.1983, 14.0116, 14.8234, 15.6359, 16.4515, 17.2723, 18.1006, 18.9386, 19.7892};
00566 static const G4double X3[nX]={.060124,
00567 5.73857, 8.62595, 10.9849, 13.0644, 14.9636, 16.7340, 18.4066, 20.0019, 21.5342, 23.0142,
00568 24.4497, 25.8471, 27.2114, 28.5467, 29.8564, 31.1434, 32.4102, 33.6589, 34.8912, 36.1095};
00569 static const G4double X4[nX]={.0992363,
00570 8.23746, 12.1036, 15.1740, 17.8231, 20.1992, 22.3792, 24.4092, 26.3198, 28.1320, 29.8615,
00571 31.5200, 33.1169, 34.6594, 36.1536, 37.6044, 39.0160, 40.3920, 41.7353, 43.0485, 44.3354};
00572 static const G4double X5[nX]={.0561127,
00573 7.33661, 10.5694, 13.0778, 15.2061, 17.0893, 18.7973, 20.3717, 21.8400, 23.2211, 24.5291,
00574 25.7745, 26.9655, 28.1087, 29.2094, 30.2721, 31.3003, 32.2972, 33.2656, 34.2076, 35.1265};
00575 static const G4double X6[nX]={.0145859,
00576 4.81774, 6.83565, 8.37399, 9.66291, 10.7920, 11.8075, 12.7366, 13.5975, 14.4025, 15.1608,
00577 15.8791, 16.5628, 17.2162, 17.8427, 18.4451, 19.0259, 19.5869, 20.1300, 20.6566, 21.1706};
00578 static const G4double X7[nX]={.00241155,
00579 2.87095, 4.02492, 4.89243, 5.61207, 6.23747, 6.79613, 7.30433, 7.77270, 8.20858, 8.61732,
00580 9.00296, 9.36863, 9.71682, 10.0495, 10.3684, 10.6749, 10.9701, 11.2550, 11.5306, 11.7982};
00581 static const G4double X8[nX]={.000316863,
00582 1.76189, 2.44632, 2.95477, 3.37292, 3.73378, 4.05420, 4.34415, 4.61009, 4.85651, 5.08666,
00583 5.30299, 5.50738, 5.70134, 5.88609, 6.06262, 6.23178, 6.39425, 6.55065, 6.70149, 6.84742};
00584 static const G4double X9[nX]={3.73544e-05,
00585 1.17106, 1.61289, 1.93763, 2.20259, 2.42976, 2.63034, 2.81094, 2.97582, 3.12796, 3.26949,
00586 3.40202, 3.52680, 3.64482, 3.75687, 3.86360, 3.96557, 4.06323, 4.15697, 4.24713, 4.33413};
00587 static const G4double XA[nX]={4.19131e-06,
00588 .849573, 1.16208, 1.38955, 1.57379, 1.73079, 1.86867, 1.99221, 2.10451, 2.20770, 2.30332,
00589 2.39252, 2.47622, 2.55511, 2.62977, 2.70066, 2.76818, 2.83265, 2.89437, 2.95355, 3.01051};
00590 static const G4double XB[nX]={4.59981e-07,
00591 .666131, .905836, 1.07880, 1.21796, 1.33587, 1.43890, 1.53080, 1.61399, 1.69011, 1.76040,
00592 1.82573, 1.88682, 1.94421, 1.99834, 2.04959, 2.09824, 2.14457, 2.18878, 2.23107, 2.27162};
00593 static const G4double XC[nX]={4.99861e-08,
00594 .556280, .752730, .893387, 1.00587, 1.10070, 1.18317, 1.25643, 1.32247, 1.38269, 1.43809,
00595 1.48941, 1.53724, 1.58203, 1.62416, 1.66391, 1.70155, 1.73728, 1.77128, 1.80371, 1.83473};
00596 static const G4double XD[nX]={5.40832e-09,
00597 .488069, .657650, .778236, .874148, .954621, 1.02432, 1.08599, 1.14138, 1.19172, 1.23787,
00598 1.28049, 1.32008, 1.35705, 1.39172, 1.42434, 1.45514, 1.48429, 1.51197, 1.53829, 1.56339};
00599 static const G4double XE[nX]={5.84029e-10,
00600 .445057, .597434, .705099, .790298, .861468, .922865, .976982, 1.02542, 1.06930, 1.10939,
00601 1.14630, 1.18050, 1.21233, 1.24208, 1.27001, 1.29630, 1.32113, 1.34462, 1.36691, 1.38812};
00602 static const G4double XF[nX]={6.30137e-11,
00603 .418735, .560003, .659168, .737230, .802138, .857898, .906854, .950515, .989915, 1.02580,
00604 1.05873, 1.08913, 1.11734, 1.14364, 1.16824, 1.19133, 1.21306, 1.23358, 1.25298, 1.27139};
00605 static const G4double XG[nX]={6.79627e-12,
00606 .405286, .539651, .633227, .706417, .766929, .818642, .863824, .903931, .939963, .972639,
00607 1.00250, 1.02995, 1.05532, 1.07887, 1.10082, 1.12134, 1.14058, 1.15867, 1.17572, 1.19183};
00608 static const G4double XH[nX]={7.32882e-13,
00609 .404391, .535199, .625259, .695036, .752243, .800752, .842823, .879906, .912994, .942802,
00610 .969862, .994583, 1.01729, 1.03823, 1.05763, 1.07566, 1.09246, 1.10816, 1.12286, 1.13667};
00611 static const G4double XI[nX]={7.90251e-14,
00612 .418084, .548382, .636489, .703728, .758106, .803630, .842633, .876608, .906576, .933269,
00613 .957233, .978886, .998556, 1.01651, 1.03295, 1.04807, 1.06201, 1.07489, 1.08683, 1.09792};
00614 static const G4double XJ[nX]={8.52083e-15,
00615 .447299, .579635, .666780, .731788, .783268, .825512, .861013, .891356, .917626, .940597,
00616 .960842, .978802, .994820, 1.00917, 1.02208, 1.03373, 1.04427, 1.05383, 1.06253, 1.07046};
00617
00618 static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
00619 X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
00620 static const G4double dX[nE]={
00621 (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
00622 (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
00623 (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
00624 (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
00625 (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
00626 static const G4double* Xl[nE]=
00627 {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
00628 static const G4double I0[nX]={0,
00629 .354631, 1.08972, 2.05138, 3.16564, 4.38343, 5.66828, 6.99127, 8.32858, 9.65998, 10.9680,
00630 12.2371, 13.4536, 14.6050, 15.6802, 16.6686, 17.5609, 18.3482, 19.0221, 19.5752, 20.0000};
00631 static const G4double I1[nX]={0,
00632 .281625, .877354, 1.67084, 2.60566, 3.64420, 4.75838, 5.92589, 7.12829, 8.34989, 9.57708,
00633 10.7978, 12.0014, 13.1781, 14.3190, 15.4162, 16.4620, 17.4496, 18.3724, 19.2245, 20.0000};
00634 static const G4double I2[nX]={0,
00635 .201909, .642991, 1.24946, 1.98463, 2.82370, 3.74802, 4.74263, 5.79509, 6.89474, 8.03228,
00636 9.19947, 10.3889, 11.5938, 12.8082, 14.0262, 15.2427, 16.4527, 17.6518, 18.8356, 20.0000};
00637 static const G4double I3[nX]={0,
00638 .140937, .461189, .920216, 1.49706, 2.17728, 2.94985, 3.80580, 4.73758, 5.73867, 6.80331,
00639 7.92637, 9.10316, 10.3294, 11.6013, 12.9150, 14.2672, 15.6548, 17.0746, 18.5239, 20.0000};
00640 static const G4double I4[nX]={0,
00641 .099161, .337358, .694560, 1.16037, 1.72761, 2.39078, 3.14540, 3.98768, 4.91433, 5.92245,
00642 7.00942, 8.17287, 9.41060, 10.7206, 12.1010, 13.5500, 15.0659, 16.6472, 18.2924, 20.0000};
00643 static const G4double I5[nX]={0,
00644 .071131, .255084, .543312, .932025, 1.41892, 2.00243, 2.68144, 3.45512, 4.32283, 5.28411,
00645 6.33859, 7.48602, 8.72621, 10.0590, 11.4844, 13.0023, 14.6128, 16.3158, 18.1115, 20.0000};
00646 static const G4double I6[nX]={0,
00647 .053692, .202354, .443946, .778765, 1.20774, 1.73208, 2.35319, 3.07256, 3.89177, 4.81249,
00648 5.83641, 6.96528, 8.20092, 9.54516, 10.9999, 12.5670, 14.2486, 16.0466, 17.9630, 20.0000};
00649 static const G4double I7[nX]={0,
00650 .043065, .168099, .376879, .672273, 1.05738, 1.53543, 2.10973, 2.78364, 3.56065, 4.44429,
00651 5.43819, 6.54610, 7.77186, 9.11940, 10.5928, 12.1963, 13.9342, 15.8110, 17.8313, 20.0000};
00652 static const G4double I8[nX]={0,
00653 .036051, .143997, .327877, .592202, .941572, 1.38068, 1.91433, 2.54746, 3.28517, 4.13277,
00654 5.09574, 6.17984, 7.39106, 8.73568, 10.2203, 11.8519, 13.6377, 15.5854, 17.7033, 20.0000};
00655 static const G4double I9[nX]={0,
00656 .030977, .125727, .289605, .528146, .846967, 1.25183, 1.74871, 2.34384, 3.04376, 3.85535,
00657 4.78594, 5.84329, 7.03567, 8.37194, 9.86163, 11.5150, 13.3430, 15.3576, 17.5719, 20.0000};
00658 static const G4double IA[nX]={0,
00659 .027129, .111420, .258935, .475812, .768320, 1.14297, 1.60661, 2.16648, 2.83034, 3.60650,
00660 4.50394, 5.53238, 6.70244, 8.02569, 9.51488, 11.1841, 13.0488, 15.1264, 17.4362, 20.0000};
00661 static const G4double IB[nX]={0,
00662 .024170, .100153, .234345, .433198, .703363, 1.05184, 1.48607, 2.01409, 2.64459, 3.38708,
00663 4.25198, 5.25084, 6.39647, 7.70319, 9.18708, 10.8663, 12.7617, 14.8968, 17.2990, 20.0000};
00664 static const G4double IC[nX]={0,
00665 .021877, .091263, .214670, .398677, .650133, .976322, 1.38510, 1.88504, 2.48555, 3.19709,
00666 4.03129, 5.00127, 6.12184, 7.40989, 8.88482, 10.5690, 12.4888, 14.6748, 17.1638, 20.0000};
00667 static const G4double ID[nX]={0,
00668 .020062, .084127, .198702, .370384, .606100, .913288, 1.30006, 1.77535, 2.34912, 3.03253,
00669 3.83822, 4.78063, 5.87634, 7.14459, 8.60791, 10.2929, 12.2315, 14.4621, 17.0320, 20.0000};
00670 static const G4double IE[nX]={0,
00671 .018547, .078104, .185102, .346090, .567998, .858331, 1.22535, 1.67824, 2.22735, 2.88443,
00672 3.66294, 4.57845, 5.64911, 6.89637, 8.34578, 10.0282, 11.9812, 14.2519, 16.8993, 20.0000};
00673 static const G4double IF[nX]={0,
00674 .017143, .072466, .172271, .323007, .531545, .805393, 1.15288, 1.58338, 2.10754, 2.73758,
00675 3.48769, 4.37450, 5.41770, 6.64092, 8.07288, 9.74894, 11.7135, 14.0232, 16.7522, 20.0000};
00676 static const G4double IG[nX]={0,
00677 .015618, .066285, .158094, .297316, .490692, .745653, 1.07053, 1.47479, 1.96931, 2.56677,
00678 3.28205, 4.13289, 5.14068, 6.33158, 7.73808, 9.40133, 11.3745, 13.7279, 16.5577, 20.0000};
00679 static const G4double IH[nX]={0,
00680 .013702, .058434, .139923, .264115, .437466, .667179, .961433, 1.32965, 1.78283, 2.33399,
00681 2.99871, 3.79596, 4.74916, 5.88771, 7.24937, 8.88367, 10.8576, 13.2646, 16.2417, 20.0000};
00682 static const G4double II[nX]={0,
00683 .011264, .048311, .116235, .220381, .366634, .561656, .813132, 1.13008, 1.52322, 2.00554,
00684 2.59296, 3.30542, 4.16834, 5.21490, 6.48964, 8.05434, 9.99835, 12.4580, 15.6567, 20.0000};
00685 static const G4double IJ[nX]={0,
00686 .008628, .037206, .089928, .171242, .286114, .440251, .640343, .894382, 1.21208, 1.60544,
00687 2.08962, 2.68414, 3.41486, 4.31700, 5.44048, 6.85936, 8.69067, 11.1358, 14.5885, 20.0000};
00688 static const G4double* Il[nE]=
00689 {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
00690 static const G4double lE[nE]={
00691 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
00692 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
00693 static const G4double lEmi=lE[0];
00694 static const G4double lEma=lE[nE-1];
00695 static const G4double dlE=(lEma-lEmi)/bE;
00696
00697 G4double Enu=lastE;
00698 G4double lEn=std::log(Enu);
00699 G4double rE=(lEn-lEmi)/dlE;
00700 G4int fE=static_cast<int>(rE);
00701 if(fE<0) fE=0;
00702 if(fE>pE)fE=pE;
00703 G4int sE=fE+1;
00704 G4double dE=rE-fE;
00705 G4double dEnu=Enu+Enu;
00706 G4double Enu2=Enu*Enu;
00707 G4double Ee=Enu-me;
00708 G4double ME=Enu*MN;
00709 G4double dME=ME+ME;
00710 G4double dEMN=(dEnu+MN)*ME;
00711 G4double MEm=ME-hme2;
00712 G4double sqE=Enu*std::sqrt(MEm*MEm-me2*MN2);
00713 G4double E2M=MN*Enu2-(Enu+MN)*hme2;
00714 G4double ymax=(E2M+sqE)/dEMN;
00715 G4double ymin=(E2M-sqE)/dEMN;
00716 G4double rmin=1.-ymin;
00717 G4double rhm2E=hme2/Enu2;
00718 G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
00719 G4double Q2ma=dME*ymax;
00720 G4double Q2nq=Ee*dMN-mcV;
00721 if(Q2ma>Q2nq) Q2ma=Q2nq;
00722
00723 G4double Rmi=Q2mi/Q2ma;
00724 G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78);
00725
00726 G4double Xmi=std::pow((shift-Rmi),power);
00727 G4double Xma=std::pow((shift-1.),power);
00728
00729 G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]);
00730 G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]);
00731 G4double rXi=(Xmi-iXmi)/idX;
00732 G4int iXi=static_cast<int>(rXi);
00733 if(iXi<0) iXi=0;
00734 if(iXi>bX) iXi=bX;
00735 G4double dXi=rXi-iXi;
00736 G4double bntil=Il[fE][iXi];
00737 G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
00738 G4double bntir=Il[sE][iXi];
00739 G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
00740 G4double inti=intil+dE*(intir-intil);
00741
00742 G4double rXa=(Xma-iXmi)/idX;
00743 G4int iXa=static_cast<int>(rXa);
00744 if(iXa<0) iXa=0;
00745 if(iXa>bX) iXa=bX;
00746 G4double dXa=rXa-iXa;
00747 G4double bntal=Il[fE][iXa];
00748 G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
00749 G4double bntar=Il[sE][iXa];
00750 G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
00751 G4double inta=intal+dE*(intar-intal);
00752
00753
00754 G4double intx=inti+(inta-inti)*G4UniformRand();
00755 G4int intc=static_cast<int>(intx);
00756 if(intc<0) intc=0;
00757 if(intc>bX) intc=bX;
00758 G4double dint=intx-intc;
00759 G4double mXl=Xl[fE][intc];
00760 G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
00761 G4double mXr=Xl[sE][intc];
00762 G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
00763 G4double X=Xlb+dE*(Xrb-Xlb);
00764 G4double R=shift-std::pow(X,pconv);
00765 G4double Q2=R*Q2ma;
00766 return Q2*GeV*GeV;
00767 }
00768
00769
00770 G4double G4QANuENuclearCrossSection::GetDirectPart(G4double Q2)
00771 {
00772 G4double f=Q2/4.62;
00773 G4double ff=f*f;
00774 G4double r=ff*ff;
00775 G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
00776
00777 return 1.-s_value*(1.-s_value/2);
00778 }
00779
00780
00781 G4double G4QANuENuclearCrossSection::GetNPartons(G4double Q2)
00782 {
00783 return 3.+.3581*std::log(1.+Q2/.04);
00784 }
00785
00786
00787 G4int G4QANuENuclearCrossSection::GetExchangePDGCode() {return -211;}