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 #include "G4QuasiElRatios.hh"
00043 #include "G4PhysicalConstants.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4Proton.hh"
00046 #include "G4Neutron.hh"
00047 #include "G4Deuteron.hh"
00048 #include "G4Triton.hh"
00049 #include "G4He3.hh"
00050 #include "G4Alpha.hh"
00051 #include "G4ThreeVector.hh"
00052 #include "G4CrossSectionDataSetRegistry.hh"
00053
00054
00055
00056 std::vector<G4double*> G4QuasiElRatios::vT;
00057 std::vector<G4double*> G4QuasiElRatios::vL;
00058 std::vector<std::pair<G4double,G4double>*> G4QuasiElRatios::vX;
00059
00060 G4QuasiElRatios::G4QuasiElRatios()
00061 {
00062
00063 PCSmanager=(G4ChipsProtonElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name());
00064
00065 NCSmanager=(G4ChipsNeutronElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name());
00066 }
00067
00068 G4QuasiElRatios::~G4QuasiElRatios()
00069 {
00070 std::vector<G4double*>::iterator pos;
00071 for(pos=vT.begin(); pos<vT.end(); pos++)
00072 { delete [] *pos; }
00073 vT.clear();
00074 for(pos=vL.begin(); pos<vL.end(); pos++)
00075 { delete [] *pos; }
00076 vL.clear();
00077
00078 std::vector<std::pair<G4double,G4double>*>::iterator pos2;
00079 for(pos2=vX.begin(); pos2<vX.end(); pos2++)
00080 { delete [] *pos2; }
00081 vX.clear();
00082 }
00083
00084
00085 G4QuasiElRatios* G4QuasiElRatios::GetPointer()
00086 {
00087 static G4QuasiElRatios theRatios;
00088 return &theRatios;
00089 }
00090
00091
00092 std::pair<G4double,G4double> G4QuasiElRatios::GetRatios(G4double pIU, G4int pPDG,
00093 G4int tgZ, G4int tgN)
00094 {
00095 G4double R=0.;
00096 G4double QF2In=1.;
00097 G4int tgA=tgZ+tgN;
00098 if(tgA<2) return std::make_pair(QF2In,R);
00099 std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN);
00100
00101 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.;
00102 else if(ElTot.second>0.)
00103 {
00104 R=ElTot.first/ElTot.second;
00105 QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN);
00106 }
00107 return std::make_pair(QF2In,R);
00108 }
00109
00110
00111 G4double G4QuasiElRatios::GetQF2IN_Ratio(G4double m_s, G4int A)
00112 {
00113 static const G4int nps=150;
00114 static const G4int mps=nps+1;
00115 static const G4double sma=150.;
00116 static const G4double ds=sma/nps;
00117 static const G4int nls=100;
00118 static const G4int mls=nls+1;
00119 static const G4double lsi=5.;
00120 static const G4double lsa=9.;
00121 static const G4double mi=std::exp(lsi);
00122 static const G4double min_s=std::exp(lsa);
00123 static const G4double dl=(lsa-lsi)/nls;
00124 static const G4double edl=std::exp(dl);
00125 static const G4double toler=.01;
00126 static G4double lastS=0.;
00127 static G4double lastR=0.;
00128
00129 static std::vector<G4int> vA;
00130 static std::vector<G4double> vH;
00131 static std::vector<G4int> vN;
00132 static std::vector<G4double> vM;
00133 static std::vector<G4int> vK;
00134
00135 static G4int lastA=0;
00136 static G4double lastH=0.;
00137 static G4int lastN=0;
00138 static G4double lastM=0.;
00139 static G4int lastK=0;
00140 static G4double* lastT=0;
00141 static G4double* lastL=0;
00142
00143 if(m_s<toler || A<2) return 1.;
00144 if(m_s>min_s) return 0.;
00145 if(A>238)
00146 {
00147 G4cout<<"-Warning-G4QuasiElRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl;
00148 return 0.;
00149 }
00150 G4int nDB=vA.size();
00151 if(nDB && lastA==A && m_s==lastS) return lastR;
00152 G4bool found=false;
00153 G4int i=-1;
00154 if(nDB) for (i=0; i<nDB; i++) if(A==vA[i])
00155 {
00156 found=true;
00157 break;
00158 }
00159 if(!nDB || !found)
00160 {
00161 lastA = A;
00162 lastT = new G4double[mps];
00163 lastN = static_cast<int>(m_s/ds)+1;
00164 if(lastN>nps)
00165 {
00166 lastN=nps;
00167 lastH=sma;
00168 }
00169 else lastH = lastN*ds;
00170 G4double sv=0;
00171 lastT[0]=1.;
00172 for(G4int j=1; j<=lastN; j++)
00173 {
00174 sv+=ds;
00175 lastT[j]=CalcQF2IN_Ratio(sv,A);
00176 }
00177 lastL=new G4double[mls];
00178 if(m_s>sma)
00179 {
00180 G4double ls=std::log(m_s);
00181 lastK = static_cast<int>((ls-lsi)/dl)+1;
00182 if(lastK>nls)
00183 {
00184 lastK=nls;
00185 lastM=lsa-lsi;
00186 }
00187 else lastM = lastK*dl;
00188 sv=mi;
00189 for(G4int j=0; j<=lastK; j++)
00190 {
00191 lastL[j]=CalcQF2IN_Ratio(sv,A);
00192 if(j!=lastK) sv*=edl;
00193 }
00194 }
00195 else
00196 {
00197 lastK = 0;
00198 lastM = 0.;
00199 }
00200 i++;
00201 vA.push_back(lastA);
00202 vH.push_back(lastH);
00203 vN.push_back(lastN);
00204 vM.push_back(lastM);
00205 vK.push_back(lastK);
00206 vT.push_back(lastT);
00207 vL.push_back(lastL);
00208 }
00209 else
00210 {
00211 lastA=vA[i];
00212 lastH=vH[i];
00213 lastN=vN[i];
00214 lastM=vM[i];
00215 lastK=vK[i];
00216 lastT=vT[i];
00217 lastL=vL[i];
00218 if(m_s>lastH)
00219 {
00220 G4int nextN=lastN+1;
00221 if(lastN<nps)
00222 {
00223 G4double sv=lastH;
00224
00225 lastN = static_cast<int>(m_s/ds)+1;
00226 if(lastN>nps)
00227 {
00228 lastN=nps;
00229 lastH=sma;
00230 }
00231 else lastH = lastN*ds;
00232
00233 for(G4int j=nextN; j<=lastN; j++)
00234 {
00235 sv+=ds;
00236 lastT[j]=CalcQF2IN_Ratio(sv,A);
00237 }
00238 }
00239 if(lastN>=nextN)
00240 {
00241 vH[i]=lastH;
00242 vN[i]=lastN;
00243 }
00244 G4int nextK=lastK+1;
00245 if(!lastK) nextK=0;
00246 if(m_s>sma && lastK<nls)
00247 {
00248 G4double sv=std::exp(lastM+lsi);
00249 G4double ls=std::log(m_s);
00250 lastK = static_cast<int>((ls-lsi)/dl)+1;
00251 if(lastK>nls)
00252 {
00253 lastK=nls;
00254 lastM=lsa-lsi;
00255 }
00256 else lastM = lastK*dl;
00257 for(G4int j=nextK; j<=lastK; j++)
00258 {
00259 sv*=edl;
00260 lastL[j]=CalcQF2IN_Ratio(sv,A);
00261 }
00262 }
00263 if(lastK>=nextK)
00264 {
00265 vM[i]=lastM;
00266 vK[i]=lastK;
00267 }
00268 }
00269 }
00270
00271 if(m_s<sma)
00272 {
00273 G4int n=static_cast<int>(m_s/ds);
00274 G4double d=m_s-n*ds;
00275 G4double v=lastT[n];
00276 lastR=v+d*(lastT[n+1]-v)/ds;
00277 }
00278 else
00279 {
00280 G4double ls=std::log(m_s)-lsi;
00281 G4int n=static_cast<int>(ls/dl);
00282 G4double d=ls-n*dl;
00283 G4double v=lastL[n];
00284 lastR=v+d*(lastL[n+1]-v)/dl;
00285 }
00286 if(lastR<0.) lastR=0.;
00287 if(lastR>1.) lastR=1.;
00288 return lastR;
00289 }
00290
00291
00292 G4double G4QuasiElRatios::CalcQF2IN_Ratio(G4double m_s, G4int A)
00293 {
00294 static const G4double C=1.246;
00295 G4double s2=m_s*m_s;
00296 G4double s4=s2*s2;
00297 G4double ss=std::sqrt(std::sqrt(m_s));
00298 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
00299 G4double E=.2644+.016/(1.+std::exp((29.54-m_s)/2.49));
00300 G4double F=ss*.1526*std::exp(-s2*ss*.0000859);
00301 return C*std::exp(-E*std::pow(G4double(A-1.),F))/std::pow(G4double(A),P);
00302 }
00303
00304
00305 std::pair<G4double,G4double> G4QuasiElRatios::CalcElTot(G4double p, G4int I)
00306 {
00307
00308 static const G4double lmi=3.5;
00309 static const G4double pbe=.0557;
00310 static const G4double pbt=.3;
00311 static const G4double pmi=.1;
00312 static const G4double pma=1000.;
00313 G4double El=0.;
00314 G4double To=0.;
00315 if(p<=0.)
00316 {
00317 G4cout<<"-Warning-G4QuasiElRatios::CalcElTot: p="<<p<<" is zero or negative"<<G4endl;
00318 return std::make_pair(El,To);
00319 }
00320 if (!I)
00321 {
00322 if(p<pmi)
00323 {
00324 G4double p2=p*p;
00325 El=1./(.00012+p2*.2);
00326 To=El;
00327 }
00328 else if(p>pma)
00329 {
00330 G4double lp=std::log(p)-lmi;
00331 G4double lp2=lp*lp;
00332 El=pbe*lp2+6.72;
00333 To=pbt*lp2+38.2;
00334 }
00335 else
00336 {
00337 G4double p2=p*p;
00338 G4double LE=1./(.00012+p2*.2);
00339 G4double lp=std::log(p)-lmi;
00340 G4double lp2=lp*lp;
00341 G4double rp2=1./p2;
00342 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
00343 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
00344 }
00345 }
00346 else if(I==1)
00347 {
00348 if(p<pmi)
00349 {
00350 G4double p2=p*p;
00351 El=1./(.00012+p2*(.051+.1*p2));
00352 To=El;
00353 }
00354 else if(p>pma)
00355 {
00356 G4double lp=std::log(p)-lmi;
00357 G4double lp2=lp*lp;
00358 El=pbe*lp2+6.72;
00359 To=pbt*lp2+38.2;
00360 }
00361 else
00362 {
00363 G4double p2=p*p;
00364 G4double LE=1./(.00012+p2*(.051+.1*p2));
00365 G4double lp=std::log(p)-lmi;
00366 G4double lp2=lp*lp;
00367 G4double rp2=1./p2;
00368 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
00369 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
00370 }
00371 }
00372 else if(I==2)
00373 {
00374 G4double lp=std::log(p);
00375 if(p<pmi)
00376 {
00377 G4double lr=lp+1.27;
00378 El=1.53/(lr*lr+.0676);
00379 To=El*3;
00380 }
00381 else if(p>pma)
00382 {
00383 G4double ld=lp-lmi;
00384 G4double ld2=ld*ld;
00385 G4double sp=std::sqrt(p);
00386 El=pbe*ld2+2.4+7./sp;
00387 To=pbt*ld2+22.3+12./sp;
00388 }
00389 else
00390 {
00391 G4double lr=lp+1.27;
00392 G4double LE=1.53/(lr*lr+.0676);
00393 G4double ld=lp-lmi;
00394 G4double ld2=ld*ld;
00395 G4double p2=p*p;
00396 G4double p4=p2*p2;
00397 G4double sp=std::sqrt(p);
00398 G4double lm=lp+.36;
00399 G4double md=lm*lm+.04;
00400 G4double lh=lp-.017;
00401 G4double hd=lh*lh+.0025;
00402 El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;
00403 To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd;
00404 }
00405 }
00406 else if(I==3)
00407 {
00408 G4double lp=std::log(p);
00409 if(p<pmi)
00410 {
00411 G4double lr=lp+1.27;
00412 G4double lr2=lr*lr;
00413 El=13./(lr2+lr2*lr2+.0676);
00414 To=El;
00415 }
00416 else if(p>pma)
00417 {
00418 G4double ld=lp-lmi;
00419 G4double ld2=ld*ld;
00420 G4double sp=std::sqrt(p);
00421 El=pbe*ld2+2.4+6./sp;
00422 To=pbt*ld2+22.3+5./sp;
00423 }
00424 else
00425 {
00426 G4double lr=lp+1.27;
00427 G4double lr2=lr*lr;
00428 G4double LE=13./(lr2+lr2*lr2+.0676);
00429 G4double ld=lp-lmi;
00430 G4double ld2=ld*ld;
00431 G4double p2=p*p;
00432 G4double p4=p2*p2;
00433 G4double sp=std::sqrt(p);
00434 G4double lm=lp-.32;
00435 G4double md=lm*lm+.0576;
00436 El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md;
00437 To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md;
00438 }
00439 }
00440 else if(I==4)
00441 {
00442
00443 if(p<pmi)
00444 {
00445 G4double psp=p*std::sqrt(p);
00446 El=5.2/psp;
00447 To=14./psp;
00448 }
00449 else if(p>pma)
00450 {
00451 G4double ld=std::log(p)-lmi;
00452 G4double ld2=ld*ld;
00453 El=pbe*ld2+2.23;
00454 To=pbt*ld2+19.5;
00455 }
00456 else
00457 {
00458 G4double ld=std::log(p)-lmi;
00459 G4double ld2=ld*ld;
00460 G4double sp=std::sqrt(p);
00461 G4double psp=p*sp;
00462 G4double p2=p*p;
00463 G4double p4=p2*p2;
00464 G4double lm=p-.39;
00465 G4double md=lm*lm+.000156;
00466 G4double lh=p-1.;
00467 G4double hd=lh*lh+.0156;
00468 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
00469 To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
00470 }
00471 }
00472 else if(I==5)
00473 {
00474 if(p<pmi)
00475 {
00476 G4double lr=p-.38;
00477 G4double lm=p-1.;
00478 G4double md=lm*lm+.372;
00479 El=.7/(lr*lr+.0676)+2./md;
00480 To=El+.6/md;
00481 }
00482 else if(p>pma)
00483 {
00484 G4double ld=std::log(p)-lmi;
00485 G4double ld2=ld*ld;
00486 El=pbe*ld2+2.23;
00487 To=pbt*ld2+19.5;
00488 }
00489 else
00490 {
00491 G4double ld=std::log(p)-lmi;
00492 G4double ld2=ld*ld;
00493 G4double lr=p-.38;
00494 G4double LE=.7/(lr*lr+.0676);
00495 G4double sp=std::sqrt(p);
00496 G4double p2=p*p;
00497 G4double p4=p2*p2;
00498 G4double lm=p-1.;
00499 G4double md=lm*lm+.372;
00500 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
00501 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
00502 }
00503 }
00504 else if(I==6)
00505 {
00506 if(p<pmi)
00507 {
00508 G4double p2=p*p;
00509 El=1./(.002+p2*(.12+p2));
00510 To=El;
00511 }
00512 else if(p>pma)
00513 {
00514 G4double lp=std::log(p)-lmi;
00515 G4double lp2=lp*lp;
00516 G4double sp=std::sqrt(p);
00517 El=(pbe*lp2+6.72)/(1.+2./sp);
00518 To=(pbt*lp2+38.2+900./sp)/(1.+27./sp);
00519 }
00520 else
00521 {
00522 G4double p2=p*p;
00523 G4double LE=1./(.002+p2*(.12+p2));
00524 G4double lp=std::log(p)-lmi;
00525 G4double lp2=lp*lp;
00526 G4double p4=p2*p2;
00527 G4double sp=std::sqrt(p);
00528 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
00529 To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4);
00530 }
00531 }
00532 else if(I==7)
00533 {
00534 if(p>pma)
00535 {
00536 G4double lp=std::log(p)-lmi;
00537 G4double lp2=lp*lp;
00538 El=pbe*lp2+6.72;
00539 To=pbt*lp2+38.2;
00540 }
00541 else
00542 {
00543 G4double ye=std::pow(p,1.25);
00544 G4double yt=std::pow(p,.35);
00545 G4double lp=std::log(p)-lmi;
00546 G4double lp2=lp*lp;
00547 El=80./(ye+1.)+pbe*lp2+6.72;
00548 To=(80./yt+.3)/yt+pbt*lp2+38.2;
00549 }
00550 }
00551 else
00552 {
00553 G4cout<<"*Error*G4QuasiElRatios::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl;
00554 G4Exception("G4QuasiElRatios::CalcElTot:","23",FatalException,"QEcrash");
00555 }
00556 if(El>To) El=To;
00557 return std::make_pair(El,To);
00558 }
00559
00560
00561 std::pair<G4double,G4double> G4QuasiElRatios::GetElTotXS(G4double p, G4int PDG, G4bool F)
00562 {
00563 G4int ind=0;
00564 G4bool kfl=true;
00565 G4bool kf=false;
00566 if(PDG==130||PDG==310)
00567 {
00568 kf=true;
00569 if(G4UniformRand()>.5) kfl=false;
00570 }
00571 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
00572 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
00573 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
00574 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
00575 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
00576 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
00577 else if ( PDG > 3000 && PDG < 3335) ind=6;
00578 else if ( PDG > -3335 && PDG < -2000) ind=7;
00579 else {
00580 G4cout<<"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG
00581 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
00582 G4Exception("G4QuasiElRatio::CalcElTotXS:","22",FatalException,"QEcrash");
00583 }
00584 return CalcElTot(p,ind);
00585 }
00586
00587
00588 std::pair<G4double,G4double> G4QuasiElRatios::FetchElTot(G4double p, G4int PDG, G4bool F)
00589 {
00590 static const G4int nlp=300;
00591 static const G4int mlp=nlp+1;
00592 static const G4double lpi=-5.;
00593 static const G4double lpa=10.;
00594 static const G4double mi=std::exp(lpi);
00595 static const G4double ma=std::exp(lpa);
00596 static const G4double dl=(lpa-lpi)/nlp;
00597 static const G4double edl=std::exp(dl);
00598
00599 static G4double lastP=0.;
00600 static G4int lastH=0;
00601 static G4bool lastF=true;
00602 static std::pair<G4double,G4double> lastR(0.,0.);
00603
00604 static std::vector<G4int> vI;
00605 static std::vector<G4double> vM;
00606 static std::vector<G4int> vK;
00607
00608 static G4int lastI=0;
00609 static G4double lastM=0.;
00610 static G4int lastK=0;
00611 static std::pair<G4double,G4double>* lastX=0;
00612
00613 G4int nDB=vI.size();
00614 if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;
00615
00616 lastH=PDG;
00617 lastF=F;
00618 G4int ind=-1;
00619
00620
00621 G4bool kfl=true;
00622 G4bool kf=false;
00623 if(PDG==130||PDG==310)
00624 {
00625 kf=true;
00626 if(G4UniformRand()>.5) kfl=false;
00627 }
00628 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
00629 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
00630 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
00631 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
00632 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
00633 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
00634 else if ( PDG > 3000 && PDG < 3335) ind=6;
00635 else if ( PDG > -3335 && PDG < -2000) ind=7;
00636 else {
00637 G4cout<<"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG
00638 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
00639 G4Exception("G4QuasiELRatio::FetchElTot:","22",FatalException,"QECrash");
00640 }
00641 if(nDB && lastI==ind && p>0. && p==lastP) return lastR;
00642
00643 if(p<=mi || p>=ma) return CalcElTot(p,ind);
00644 G4bool found=false;
00645 G4int i=-1;
00646 if(nDB) for (i=0; i<nDB; i++) if(ind==vI[i])
00647 {
00648 found=true;
00649 break;
00650 }
00651 G4double lp=std::log(p);
00652 if(!nDB || !found)
00653 {
00654 lastX = new std::pair<G4double,G4double>[mlp];
00655 lastI = ind;
00656 lastK = static_cast<int>((lp-lpi)/dl)+1;
00657 if(lastK>nlp)
00658 {
00659 lastK=nlp;
00660 lastM=lpa-lpi;
00661 }
00662 else lastM = lastK*dl;
00663 G4double pv=mi;
00664 for(G4int j=0; j<=lastK; j++)
00665 {
00666 lastX[j]=CalcElTot(pv,ind);
00667 if(j!=lastK) pv*=edl;
00668 }
00669 i++;
00670 vI.push_back(lastI);
00671 vM.push_back(lastM);
00672 vK.push_back(lastK);
00673 vX.push_back(lastX);
00674 }
00675 else
00676 {
00677 lastI=vI[i];
00678 lastM=vM[i];
00679 lastK=vK[i];
00680 lastX=vX[i];
00681 G4int nextK=lastK+1;
00682 G4double lpM=lastM+lpi;
00683 if(lp>lpM && lastK<nlp)
00684 {
00685 lastK = static_cast<int>((lp-lpi)/dl)+1;
00686 if(lastK>nlp)
00687 {
00688 lastK=nlp;
00689 lastM=lpa-lpi;
00690 }
00691 else lastM = lastK*dl;
00692 G4double pv=std::exp(lpM);
00693 for(G4int j=nextK; j<=lastK; j++)
00694 {
00695 pv*=edl;
00696 lastX[j]=CalcElTot(pv,ind);
00697 }
00698 }
00699 if(lastK>=nextK)
00700 {
00701 vM[i]=lastM;
00702 vK[i]=lastK;
00703 }
00704 }
00705
00706 G4double dlp=lp-lpi;
00707 G4int n=static_cast<int>(dlp/dl);
00708 G4double d=dlp-n*dl;
00709 G4double e=lastX[n].first;
00710 lastR.first=e+d*(lastX[n+1].first-e)/dl;
00711 if(lastR.first<0.) lastR.first = 0.;
00712 G4double t=lastX[n].second;
00713 lastR.second=t+d*(lastX[n+1].second-t)/dl;
00714 if(lastR.second<0.) lastR.second= 0.;
00715 if(lastR.first>lastR.second) lastR.first = lastR.second;
00716 return lastR;
00717 }
00718
00719
00720 std::pair<G4double,G4double> G4QuasiElRatios::GetElTot(G4double pIU, G4int hPDG,
00721 G4int Z, G4int N)
00722 {
00723 G4double pGeV=pIU/gigaelectronvolt;
00724 if(Z<1 && N<1)
00725 {
00726 G4cout<<"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
00727 return std::make_pair(0.,0.);
00728 }
00729 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
00730 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
00731 G4double A=(Z+N)/millibarn;
00732 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
00733 }
00734
00735
00736 std::pair<G4double,G4double> G4QuasiElRatios::GetChExFactor(G4double pIU, G4int hPDG,
00737 G4int Z, G4int N)
00738 {
00739 G4double pGeV=pIU/gigaelectronvolt;
00740 G4double resP=0.;
00741 G4double resN=0.;
00742 if(Z<1 && N<1)
00743 {
00744 G4cout<<"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
00745 return std::make_pair(resP,resN);
00746 }
00747 G4double A=Z+N;
00748 G4double pf=0.;
00749 G4double nf=0.;
00750 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N);
00751 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
00752 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
00753 {
00754 G4double dA=A+A;
00755 pf=Z/(dA+N+N);
00756 nf=N/(dA+Z+Z);
00757 }
00758 G4double mult=1.;
00759 if(pGeV>.5)
00760 {
00761 mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
00762 if(mult>1.) mult=1.;
00763 }
00764 if(pf)
00765 {
00766 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
00767 resP=pf*(hp.second/hp.first-1.)*mult;
00768 }
00769 if(nf)
00770 {
00771 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
00772 resN=nf*(hn.second/hn.first-1.)*mult;
00773 }
00774 return std::make_pair(resP,resN);
00775 }
00776
00777
00778
00779 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiElRatios::Scatter(G4int NPDG,
00780 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
00781 {
00782 static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
00783 static const G4double mProt= G4Proton::Proton()->GetPDGMass();
00784 static const G4double mDeut= G4Deuteron::Deuteron()->GetPDGMass();
00785 static const G4double mTrit= G4Triton::Triton()->GetPDGMass();
00786 static const G4double mHel3= G4He3::He3()->GetPDGMass();
00787 static const G4double mAlph= G4Alpha::Alpha()->GetPDGMass();
00788
00789 G4LorentzVector pr4M=p4M/megaelectronvolt;
00790 N4M/=megaelectronvolt;
00791 G4LorentzVector tot4M=N4M+p4M;
00792 G4double mT=mNeut;
00793 G4int Z=0;
00794 G4int N=1;
00795 if(NPDG==2212||NPDG==90001000)
00796 {
00797 mT=mProt;
00798 Z=1;
00799 N=0;
00800 }
00801 else if(NPDG==90001001)
00802 {
00803 mT=mDeut;
00804 Z=1;
00805 N=1;
00806 }
00807 else if(NPDG==90002001)
00808 {
00809 mT=mHel3;
00810 Z=2;
00811 N=1;
00812 }
00813 else if(NPDG==90001002)
00814 {
00815 mT=mTrit;
00816 Z=1;
00817 N=2;
00818 }
00819 else if(NPDG==90002002)
00820 {
00821 mT=mAlph;
00822 Z=2;
00823 N=2;
00824 }
00825 else if(NPDG!=2112&&NPDG!=90000001)
00826 {
00827 G4cout<<"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
00828 G4Exception("G4QuasiElRatios::Scatter:","21",FatalException,"QEcomplain");
00829
00830 }
00831 G4double mT2=mT*mT;
00832 G4double mP2=pr4M.m2();
00833 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
00834 G4double E2=E*E;
00835 if(E<0. || E2<mP2)
00836 {
00837 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00838 }
00839 G4double P=std::sqrt(E2-mP2);
00840
00841 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QE::Scatter: pPDG="<<pPDG<<G4endl;
00842 G4int PDG=2212;
00843 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
00844 if(!Z && N==1)
00845 {
00846 Z=1;
00847 N=0;
00848 if (PDG==2212) PDG=2112;
00849 else if(PDG==2112) PDG=2212;
00850 }
00851 G4double xSec=0.;
00852 if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG);
00853 else xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG);
00854
00855 if(xSec <= 0.)
00856 {
00857 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00858 }
00859 G4double mint=0.;
00860 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);
00861 else mint=NCSmanager->GetExchangeT(Z,N,PDG);
00862 G4double maxt=0.;
00863 if(PDG==2212) maxt=PCSmanager->GetHMaxT();
00864 else maxt=NCSmanager->GetHMaxT();
00865 G4double cost=1.-(mint+mint)/maxt;
00866 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00867 {
00868 if (cost>1.) cost=1.;
00869 else if(cost<-1.) cost=-1.;
00870 else
00871 {
00872 G4double tm=0.;
00873 if(PDG==2212) tm=PCSmanager->GetHMaxT();
00874 else tm=NCSmanager->GetHMaxT();
00875 G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
00876 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00877 }
00878 }
00879 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);
00880 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00881 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
00882 {
00883 G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
00884
00885 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00886 }
00887 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt);
00888 }
00889
00890
00891
00892
00893 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiElRatios::ChExer(G4int NPDG,
00894 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
00895 {
00896 static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
00897 static const G4double mProt= G4Proton::Proton()->GetPDGMass();
00898 G4LorentzVector pr4M=p4M/megaelectronvolt;
00899 N4M/=megaelectronvolt;
00900 G4LorentzVector tot4M=N4M+p4M;
00901 G4int Z=0;
00902 G4int N=1;
00903 G4int sPDG=0;
00904 G4double mS=0.;
00905 G4double mT=mProt;
00906 if(NPDG==2212)
00907 {
00908 mT=mNeut;
00909 Z=1;
00910 N=0;
00911 if(pPDG==-211) sPDG=111;
00912 else if(pPDG==-321)
00913 {
00914 sPDG=310;
00915 if(G4UniformRand()>.5) sPDG=130;
00916 }
00917 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321;
00918 else if(pPDG==3112) sPDG=3212;
00919 else if(pPDG==3212) sPDG=3222;
00920 else if(pPDG==3312) sPDG=3322;
00921 }
00922 else if(NPDG==2112)
00923 {
00924 if(pPDG==211) sPDG=111;
00925 else if(pPDG==321)
00926 {
00927 sPDG=310;
00928 if(G4UniformRand()>.5) sPDG=130;
00929 }
00930 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321;
00931 else if(pPDG==3222) sPDG=3212;
00932 else if(pPDG==3212) sPDG=3112;
00933 else if(pPDG==3322) sPDG=3312;
00934 }
00935 else
00936 {
00937 G4cout<<"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
00938 G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
00939
00940 }
00941 if(sPDG) mS=mNeut;
00942 else
00943 {
00944 G4cout<<"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
00945 G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
00946
00947 }
00948 G4double mT2=mT*mT;
00949 G4double mS2=mS*mS;
00950 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
00951 G4double E2=E*E;
00952 if(E<0. || E2<mS2)
00953 {
00954 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00955 }
00956 G4double P=std::sqrt(E2-mS2);
00957
00958 G4int PDG=2212;
00959 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
00960 if(!Z && N==1)
00961 {
00962 Z=1;
00963 N=0;
00964 if (PDG==2212) PDG=2112;
00965 else if(PDG==2112) PDG=2212;
00966 }
00967 G4double xSec=0.;
00968 if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG);
00969 else xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG);
00970
00971 if(xSec <= 0.)
00972 {
00973 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00974 }
00975 G4double mint=0.;
00976 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);
00977 else mint=NCSmanager->GetExchangeT(Z,N,PDG);
00978 G4double maxt=0.;
00979 if(PDG==2212) maxt=PCSmanager->GetHMaxT();
00980 else maxt=NCSmanager->GetHMaxT();
00981 G4double cost=1.-mint/maxt;
00982 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00983 {
00984 if (cost>1.) cost=1.;
00985 else if(cost<-1.) cost=-1.;
00986 else
00987 {
00988 G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
00989 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00990 }
00991 }
00992 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);
00993 pr4M=G4LorentzVector(0.,0.,0.,mS);
00994 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00995 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
00996 {
00997 G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
00998
00999 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
01000 }
01001 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt);
01002 }
01003
01004
01005 G4double G4QuasiElRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG)
01006 {
01007 p/=MeV;
01008 G4double A=Z+N;
01009 if(A<1.5) return 0.;
01010 G4double C=0.;
01011 if (pPDG==2212) C=N/(A+Z);
01012 else if(pPDG==2112) C=Z/(A+N);
01013 else G4cout<<"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
01014 C*=C;
01015
01016 G4double sp=std::sqrt(p);
01017 G4double p2=p*p;
01018 G4double p4=p2*p2;
01019 G4double dl1=std::log(p)-5.;
01020 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
01021 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p;
01022 G4double R=U/T;
01023 return C*R*R;
01024 }
01025
01026
01027 G4bool G4QuasiElRatios::RelDecayIn2(G4LorentzVector& theMomentum, G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
01028 G4LorentzVector& dir, G4double maxCost, G4double minCost)
01029 {
01030 G4double fM2 = f4Mom.m2();
01031 G4double fM = std::sqrt(fM2);
01032 G4double sM2 = s4Mom.m2();
01033 G4double sM = std::sqrt(sM2);
01034 G4double iM2 = theMomentum.m2();
01035 G4double iM = std::sqrt(iM2);
01036 G4double vP = theMomentum.rho();
01037 G4double dE = theMomentum.e();
01038 if(dE<vP)
01039 {
01040 G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
01041 G4double accuracy=.000001*vP;
01042 G4double emodif=std::fabs(dE-vP);
01043
01044
01045 G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
01046 theMomentum.setE(vP+emodif+.01*accuracy);
01047
01048 }
01049 G4ThreeVector ltb = theMomentum.boostVector();
01050 G4ThreeVector ltf = -ltb;
01051 G4LorentzVector cdir = dir;
01052 cdir.boost(ltf);
01053 G4ThreeVector vdir = cdir.vect();
01054 G4ThreeVector vx(0.,0.,1.);
01055 G4ThreeVector vy(0.,1.,0.);
01056 G4ThreeVector vz(1.,0.,0.);
01057 if(vdir.mag2() > 0.)
01058 {
01059 vx = vdir.unit();
01060 G4ThreeVector vv= vx.orthogonal();
01061 vy = vv.unit();
01062 vz = vx.cross(vy);
01063 }
01064 if(maxCost> 1.) maxCost= 1.;
01065 if(minCost<-1.) minCost=-1.;
01066 if(maxCost<-1.) maxCost=-1.;
01067 if(minCost> 1.) minCost= 1.;
01068 if(minCost> maxCost) minCost=maxCost;
01069 if(std::fabs(iM-fM-sM)<.00000001)
01070 {
01071 G4double fR=fM/iM;
01072 G4double sR=sM/iM;
01073 f4Mom=fR*theMomentum;
01074 s4Mom=sR*theMomentum;
01075 return true;
01076 }
01077 else if (iM+.001<fM+sM || iM==0.)
01078 {
01079 G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
01080 return false;
01081 }
01082 G4double d2 = iM2-fM2-sM2;
01083 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;
01084 if(p2<0.)
01085 {
01086 p2=0.;
01087 }
01088 G4double p = std::sqrt(p2);
01089 G4double ct = maxCost;
01090 if(maxCost>minCost)
01091 {
01092 G4double dcost=maxCost-minCost;
01093 ct = minCost+dcost*G4UniformRand();
01094 }
01095 G4double phi= twopi*G4UniformRand();
01096 G4double ps=0.;
01097 if(std::fabs(ct)<1.) ps = p * std::sqrt(1.-ct*ct);
01098 else
01099 {
01100 if(ct>1.) ct=1.;
01101 if(ct<-1.) ct=-1.;
01102 }
01103 G4ThreeVector pVect=(ps*std::sin(phi))*vz+(ps*std::cos(phi))*vy+p*ct*vx;
01104
01105 f4Mom.setVect(pVect);
01106 f4Mom.setE(std::sqrt(fM2+p2));
01107 s4Mom.setVect((-1)*pVect);
01108 s4Mom.setE(std::sqrt(sM2+p2));
01109
01110 if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
01111 <<f4Mom.e()-f4Mom.rho()<<G4endl;
01112 f4Mom.boost(ltb);
01113 if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
01114 <<s4Mom.e()-s4Mom.rho()<<G4endl;
01115 s4Mom.boost(ltb);
01116 return true;
01117 }
01118
01119
01120
01121
01122
01123