70 std::vector<G4double*>::iterator
pos;
71 for(pos=vT.begin(); pos<vT.end(); pos++)
74 for(pos=vL.begin(); pos<vL.end(); pos++)
78 std::vector<std::pair<G4double,G4double>*>::iterator pos2;
79 for(pos2=vX.begin(); pos2<vX.end(); pos2++)
98 if(tgA<2)
return std::make_pair(QF2In,R);
99 std::pair<G4double,G4double> ElTot=
GetElTot(pIU, pPDG, tgZ, tgN);
101 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.;
102 else if(ElTot.second>0.)
104 R=ElTot.first/ElTot.second;
107 return std::make_pair(QF2In,R);
113 static const G4int nps=150;
114 static const G4int mps=nps+1;
117 static const G4int nls=100;
118 static const G4int mls=nls+1;
121 static const G4double mi=std::exp(lsi);
122 static const G4double min_s=std::exp(lsa);
123 static const G4double dl=(lsa-lsi)/nls;
124 static const G4double edl=std::exp(dl);
129 static G4ThreadLocal std::vector<G4int> *vA_G4MT_TLS_ = 0 ;
if (!vA_G4MT_TLS_) vA_G4MT_TLS_ =
new std::vector<G4int> ; std::vector<G4int> &vA = *vA_G4MT_TLS_;
130 static G4ThreadLocal std::vector<G4double> *vH_G4MT_TLS_ = 0 ;
if (!vH_G4MT_TLS_) vH_G4MT_TLS_ =
new std::vector<G4double> ; std::vector<G4double> &vH = *vH_G4MT_TLS_;
131 static G4ThreadLocal std::vector<G4int> *vN_G4MT_TLS_ = 0 ;
if (!vN_G4MT_TLS_) vN_G4MT_TLS_ =
new std::vector<G4int> ; std::vector<G4int> &vN = *vN_G4MT_TLS_;
132 static G4ThreadLocal std::vector<G4double> *vM_G4MT_TLS_ = 0 ;
if (!vM_G4MT_TLS_) vM_G4MT_TLS_ =
new std::vector<G4double> ; std::vector<G4double> &vM = *vM_G4MT_TLS_;
133 static G4ThreadLocal std::vector<G4int> *vK_G4MT_TLS_ = 0 ;
if (!vK_G4MT_TLS_) vK_G4MT_TLS_ =
new std::vector<G4int> ; std::vector<G4int> &vK = *vK_G4MT_TLS_;
143 if(m_s<toler || A<2)
return 1.;
144 if(m_s>min_s)
return 0.;
147 G4cout<<
"-Warning-G4QuasiElRatio::GetQF2IN_Ratio:A="<<A<<
">238, return zero"<<
G4endl;
151 if(nDB && lastA==A && m_s==lastS)
return lastR;
154 if(nDB)
for (i=0; i<nDB; i++)
if(A==vA[i])
163 lastN =
static_cast<int>(m_s/ds)+1;
169 else lastH = lastN*ds;
172 for(
G4int j=1; j<=lastN; j++)
181 lastK =
static_cast<int>((ls-lsi)/dl)+1;
187 else lastM = lastK*dl;
189 for(
G4int j=0; j<=lastK; j++)
192 if(j!=lastK) sv*=edl;
225 lastN =
static_cast<int>(m_s/ds)+1;
231 else lastH = lastN*ds;
233 for(
G4int j=nextN; j<=lastN; j++)
246 if(m_s>sma && lastK<nls)
250 lastK =
static_cast<int>((ls-lsi)/dl)+1;
256 else lastM = lastK*dl;
257 for(
G4int j=nextK; j<=lastK; j++)
273 G4int n=
static_cast<int>(m_s/ds);
276 lastR=v+d*(lastT[n+1]-v)/ds;
281 G4int n=
static_cast<int>(ls/dl);
284 lastR=v+d*(lastL[n+1]-v)/dl;
286 if(lastR<0.) lastR=0.;
287 if(lastR>1.) lastR=1.;
297 G4double ss=std::sqrt(std::sqrt(m_s));
298 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
299 G4double E=.2644+.016/(1.+std::exp((29.54-m_s)/2.49));
300 G4double F=ss*.1526*std::exp(-s2*ss*.0000859);
317 G4cout<<
"-Warning-G4QuasiElRatios::CalcElTot: p="<<p<<
" is zero or negative"<<
G4endl;
318 return std::make_pair(El,To);
325 El=1./(.00012+p2*.2);
342 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
343 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
351 El=1./(.00012+p2*(.051+.1*p2));
368 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
369 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
378 El=1.53/(lr*lr+.0676);
386 El=pbe*ld2+2.4+7./
sp;
387 To=pbt*ld2+22.3+12./
sp;
402 El=LE+(pbe*ld2+2.4+7./
sp)/(1.+.7/p4)+.6/md+.05/hd;
403 To=LE*3+(pbt*ld2+22.3+12./
sp)/(1.+.4/p4)+1./md+.06/hd;
413 El=13./(lr2+lr2*lr2+.0676);
421 El=pbe*ld2+2.4+6./
sp;
422 To=pbt*ld2+22.3+5./
sp;
436 El=LE+(pbe*ld2+2.4+6./
sp)/(1.+3./p4)+.7/md;
437 To=LE+(pbt*ld2+22.3+5./
sp)/(1.+1./p4)+.8/md;
468 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
469 To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
479 El=.7/(lr*lr+.0676)+2./md;
500 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
501 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
509 El=1./(.002+p2*(.12+p2));
517 El=(pbe*lp2+6.72)/(1.+2./
sp);
518 To=(pbt*lp2+38.2+900./
sp)/(1.+27./sp);
528 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
529 To=LE+(pbt*lp2+38.2+900./
sp)/(1.+27./sp+3./p4);
547 El=80./(ye+1.)+pbe*lp2+6.72;
548 To=(80./yt+.3)/yt+pbt*lp2+38.2;
553 G4cout<<
"*Error*G4QuasiElRatios::CalcElTot:ind="<<I<<
" is not defined (0-7)"<<
G4endl;
557 return std::make_pair(El,To);
566 if(PDG==130||PDG==310)
571 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
572 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
573 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
574 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
575 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
576 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
577 else if ( PDG > 3000 && PDG < 3335) ind=6;
578 else if ( PDG > -3335 && PDG < -2000) ind=7;
580 G4cout<<
"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG
581 <<
", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<
G4endl;
590 static const G4int nlp=300;
591 static const G4int mlp=nlp+1;
594 static const G4double mi=std::exp(lpi);
595 static const G4double ma=std::exp(lpa);
596 static const G4double dl=(lpa-lpi)/nlp;
597 static const G4double edl=std::exp(dl);
602 static G4ThreadLocal std::pair<G4double,G4double>* G4MT_lastR=0;
if (!G4MT_lastR) G4MT_lastR =
new std::pair<G4double,G4double>(0.,0.); std::pair<G4double,G4double> &lastR=*G4MT_lastR;
604 static G4ThreadLocal std::vector<G4int> *vI_G4MT_TLS_ = 0 ;
if (!vI_G4MT_TLS_) vI_G4MT_TLS_ =
new std::vector<G4int> ; std::vector<G4int> &vI = *vI_G4MT_TLS_;
605 static G4ThreadLocal std::vector<G4double> *vM_G4MT_TLS_ = 0 ;
if (!vM_G4MT_TLS_) vM_G4MT_TLS_ =
new std::vector<G4double> ; std::vector<G4double> &vM = *vM_G4MT_TLS_;
606 static G4ThreadLocal std::vector<G4int> *vK_G4MT_TLS_ = 0 ;
if (!vK_G4MT_TLS_) vK_G4MT_TLS_ =
new std::vector<G4int> ; std::vector<G4int> &vK = *vK_G4MT_TLS_;
611 static G4ThreadLocal std::pair<G4double,G4double>* *lastX_G4MT_TLS_ = 0 ;
if (!lastX_G4MT_TLS_) {lastX_G4MT_TLS_ =
new std::pair<G4double,G4double>* ; *lastX_G4MT_TLS_=0 ; } std::pair<G4double,G4double>* &lastX = *lastX_G4MT_TLS_;
614 if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP)
return lastR;
623 if(PDG==130||PDG==310)
628 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
629 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
630 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
631 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
632 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
633 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
634 else if ( PDG > 3000 && PDG < 3335) ind=6;
635 else if ( PDG > -3335 && PDG < -2000) ind=7;
637 G4cout<<
"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG
638 <<
", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<
G4endl;
641 if(nDB && lastI==ind && p>0. && p==lastP)
return lastR;
643 if(p<=mi || p>=ma)
return CalcElTot(p,ind);
646 if(nDB)
for (i=0; i<nDB; i++)
if(ind==vI[i])
654 lastX =
new std::pair<G4double,G4double>[mlp];
656 lastK =
static_cast<int>((lp-lpi)/dl)+1;
662 else lastM = lastK*dl;
664 for(
G4int j=0; j<=lastK; j++)
667 if(j!=lastK) pv*=edl;
683 if(lp>lpM && lastK<nlp)
685 lastK =
static_cast<int>((lp-lpi)/dl)+1;
691 else lastM = lastK*dl;
693 for(
G4int j=nextK; j<=lastK; j++)
707 G4int n=
static_cast<int>(dlp/dl);
710 lastR.first=e+d*(lastX[n+1].first-e)/dl;
711 if(lastR.first<0.) lastR.first = 0.;
713 lastR.second=t+d*(lastX[n+1].second-t)/dl;
714 if(lastR.second<0.) lastR.second= 0.;
715 if(lastR.first>lastR.second) lastR.first = lastR.second;
726 G4cout<<
"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<
",N="<<N<<
", return zero"<<
G4endl;
727 return std::make_pair(0.,0.);
729 std::pair<G4double,G4double> hp=
FetchElTot(pGeV, hPDG,
true);
730 std::pair<G4double,G4double> hn=
FetchElTot(pGeV, hPDG,
false);
732 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
744 G4cout<<
"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<
",N="<<N<<
", return zero"<<
G4endl;
745 return std::make_pair(resP,resN);
750 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N);
751 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
752 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
761 mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
766 std::pair<G4double,G4double> hp=
FetchElTot(pGeV, hPDG,
true);
767 resP=pf*(hp.second/hp.first-1.)*mult;
771 std::pair<G4double,G4double> hn=
FetchElTot(pGeV, hPDG,
false);
772 resN=nf*(hn.second/hn.first-1.)*mult;
774 return std::make_pair(resP,resN);
795 if(NPDG==2212||NPDG==90001000)
801 else if(NPDG==90001001)
807 else if(NPDG==90002001)
813 else if(NPDG==90001002)
819 else if(NPDG==90002002)
825 else if(NPDG!=2112&&NPDG!=90000001)
827 G4cout<<
"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<
" is not 2212 or 2112"<<
G4endl;
833 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
841 if(pPDG>3400 || pPDG<-3400)
G4cout<<
"-Warning-G4QE::Scatter: pPDG="<<pPDG<<
G4endl;
843 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
848 if (PDG==2212) PDG=2112;
849 else if(PDG==2112) PDG=2212;
866 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
868 if (cost>1.) cost=1.;
869 else if(cost<-1.) cost=-1.;
875 G4cerr<<
"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<
",-t="<<mint<<
",tm="<<tm<<
G4endl;
881 if(!
RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
883 G4cerr<<
"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<
",mT="<<mT<<
",mP="<<std::sqrt(mP2)<<
G4endl;
911 if(pPDG==-211) sPDG=111;
917 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321;
918 else if(pPDG==3112) sPDG=3212;
919 else if(pPDG==3212) sPDG=3222;
920 else if(pPDG==3312) sPDG=3322;
924 if(pPDG==211) sPDG=111;
930 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321;
931 else if(pPDG==3222) sPDG=3212;
932 else if(pPDG==3212) sPDG=3112;
933 else if(pPDG==3322) sPDG=3312;
937 G4cout<<
"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<
" is not 2212 or 2112"<<
G4endl;
944 G4cout<<
"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<
", NPDG="<<NPDG<<
G4endl;
950 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
959 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
964 if (PDG==2212) PDG=2112;
965 else if(PDG==2112) PDG=2212;
982 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
984 if (cost>1.) cost=1.;
985 else if(cost<-1.) cost=-1.;
988 G4cerr<<
"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<
",t="<<mint<<
",tm="<<maxt<<
G4endl;
995 if(!
RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
997 G4cerr<<
"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<
",mT="<<mT<<
",mP="<<mS<<
G4endl;
1009 if(A<1.5)
return 0.;
1011 if (pPDG==2212) C=N/(A+Z);
1012 else if(pPDG==2112) C=Z/(A+N);
1013 else G4cout<<
"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<
G4endl;
1020 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
1021 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p;
1040 G4cerr<<
"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<
", E-p="<<dE-vP<<
G4endl;
1045 G4cerr<<
"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<
G4endl;
1046 theMomentum.setE(vP+emodif+.01*accuracy);
1057 if(vdir.mag2() > 0.)
1064 if(maxCost> 1.) maxCost= 1.;
1065 if(minCost<-1.) minCost=-1.;
1066 if(maxCost<-1.) maxCost=-1.;
1067 if(minCost> 1.) minCost= 1.;
1068 if(minCost> maxCost) minCost=maxCost;
1069 if(std::fabs(iM-fM-sM)<.00000001)
1073 f4Mom=fR*theMomentum;
1074 s4Mom=sR*theMomentum;
1077 else if (iM+.001<fM+sM || iM==0.)
1079 G4cerr<<
"***G4QH::RelDecIn2: fM="<<fM<<
"+sM="<<sM<<
">iM="<<iM<<
",d="<<iM-fM-sM<<
G4endl;
1083 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;
1097 if(std::fabs(ct)<1.) ps = p * std::sqrt(1.-ct*ct);
1103 G4ThreeVector pVect=(ps*std::sin(phi))*vz+(ps*std::cos(phi))*vy+p*ct*vx;
1105 f4Mom.setVect(pVect);
1106 f4Mom.setE(std::sqrt(fM2+p2));
1107 s4Mom.setVect((-1)*pVect);
1108 s4Mom.setE(std::sqrt(sM2+p2));
1110 if(f4Mom.e()+.001<f4Mom.rho())
G4cerr<<
"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<
",e-p="
1111 <<f4Mom.e()-f4Mom.rho()<<
G4endl;
1113 if(s4Mom.e()+.001<s4Mom.rho())
G4cerr<<
"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<
",e-p="
1114 <<s4Mom.e()-s4Mom.rho()<<
G4endl;
G4VCrossSectionDataSet * GetCrossSectionDataSet(const G4String &name, G4bool warning=true)
std::pair< G4double, G4double > GetRatios(G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
static G4QuasiElRatios * GetPointer()
static const double megaelectronvolt
static G4ThreadLocal std::vector< std::pair< G4double, G4double > * > * vX_G4MT_TLS_
CLHEP::Hep3Vector G4ThreeVector
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
G4ChipsProtonElasticXS * PCSmanager
std::pair< G4double, G4double > FetchElTot(G4double pGeV, G4int PDG, G4bool F)
std::pair< G4double, G4double > GetElTot(G4double pIU, G4int hPDG, G4int Z, G4int N)
std::pair< G4LorentzVector, G4LorentzVector > Scatter(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
std::pair< G4double, G4double > GetChExFactor(G4double pIU, G4int pPDG, G4int Z, G4int N)
G4GLOB_DLL std::ostream G4cout
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
static const char * Default_Name()
static G4CrossSectionDataSetRegistry * Instance()
static G4Triton * Triton()
static G4Proton * Proton()
G4double GetQF2IN_Ratio(G4double TotCS_mb, G4int A)
G4double CalcQF2IN_Ratio(G4double TCSmb, G4int A)
static G4Neutron * Neutron()
static const G4double A[nN]
static G4Deuteron * Deuteron()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
std::pair< G4LorentzVector, G4LorentzVector > ChExer(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
static G4ThreadLocal std::vector< G4double * > * vL_G4MT_TLS_
G4double ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG)
G4bool RelDecayIn2(G4LorentzVector &theMomentum, G4LorentzVector &f4Mom, G4LorentzVector &s4Mom, G4LorentzVector &dir, G4double maxCost=1., G4double minCost=-1.)
G4double GetPDGMass() const
std::pair< G4double, G4double > CalcElTot(G4double pGeV, G4int Index)
static const double gigaelectronvolt
static const double millibarn
static const char * Default_Name()
std::pair< G4double, G4double > GetElTotXS(G4double Mom, G4int PDG, G4bool F)
static G4ThreadLocal std::vector< G4double * > * vT_G4MT_TLS_
static const G4double pos
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
G4GLOB_DLL std::ostream G4cerr
G4ChipsNeutronElasticXS * NCSmanager
CLHEP::HepLorentzVector G4LorentzVector