59 const G4int mps=nps+1;
63 const G4int mls=nls+1;
78 const G4int mlp=nlp+1;
90 vT =
new std::vector<G4double*>;
91 vL =
new std::vector<G4double*>;
92 vX =
new std::vector<std::pair<G4double,G4double>*>;
109 lastXtot =
new std::pair<G4double,G4double>;
118 std::vector<G4double*>::iterator
pos;
119 for(pos=
vT->begin(); pos<
vT->end(); pos++)
124 for(pos=
vL->begin(); pos<
vL->end(); pos++)
129 std::vector<std::pair<G4double,G4double>*>::iterator pos2;
130 for(pos2=
vX->begin(); pos2<
vX->end(); pos2++)
144 if(tgA<2)
return std::make_pair(QF2In,R);
145 std::pair<G4double,G4double> ElTot=
GetElTot(pIU, pPDG, tgZ, tgN);
147 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.;
148 else if(ElTot.second>0.)
150 R=ElTot.first/ElTot.second;
153 return std::make_pair(QF2In,R);
160 if(m_s<toler || A<2)
return 1.;
161 if(m_s>min_s)
return 0.;
165 G4cout<<
"-Warning-G4QuasiElRatio::GetQF2IN_Ratio:A="<<A<<
">238, return zero"<<
G4endl;
172 if(nDB)
for (i=0; i<nDB; i++)
if(A==
vARatio[i])
199 lastKRatio =
static_cast<int>((ls-lsi)/dls)+1;
268 lastKRatio =
static_cast<int>((ls-lsi)/dls)+1;
291 G4int n=
static_cast<int>(m_s/ds);
299 G4int n=
static_cast<int>(ls/dls);
314 G4double ss=std::sqrt(std::sqrt(m_s));
315 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
330 G4cout<<
"-Warning-G4QuasiElRatios::CalcElTot: p="<<p<<
" is zero or negative"<<
G4endl;
331 return std::make_pair(El,To);
338 El=1./(.00012+p2*.2);
355 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
356 To=LE+(pbt*lp2+38.2+52.7*
rp2)/(1.+2.72*rp2*rp2);
364 El=1./(.00012+p2*(.051+.1*p2));
381 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
382 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*
rp2);
391 El=1.53/(lr*lr+.0676);
399 El=pbe*ld2+2.4+7./
sp;
400 To=pbt*ld2+22.3+12./
sp;
415 El=LE+(pbe*ld2+2.4+7./
sp)/(1.+.7/p4)+.6/md+.05/hd;
416 To=LE*3+(pbt*ld2+22.3+12./
sp)/(1.+.4/p4)+1./md+.06/hd;
426 El=13./(lr2+lr2*lr2+.0676);
434 El=pbe*ld2+2.4+6./
sp;
435 To=pbt*ld2+22.3+5./
sp;
449 El=LE+(pbe*ld2+2.4+6./
sp)/(1.+3./p4)+.7/md;
450 To=LE+(pbt*ld2+22.3+5./
sp)/(1.+1./p4)+.8/md;
481 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
482 To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
492 El=.7/(lr*lr+.0676)+2./md;
513 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
514 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
522 El=1./(.002+p2*(.12+p2));
530 El=(pbe*lp2+6.72)/(1.+2./
sp);
531 To=(pbt*lp2+38.2+900./
sp)/(1.+27./sp);
541 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
542 To=LE+(pbt*lp2+38.2+900./
sp)/(1.+27./sp+3./p4);
560 El=80./(ye+1.)+pbe*lp2+6.72;
561 To=(80./yt+.3)/yt+pbt*lp2+38.2;
566 G4cout<<
"*Error*G4QuasiElRatios::CalcElTot:ind="<<I<<
" is not defined (0-7)"<<
G4endl;
570 return std::make_pair(El,To);
579 if(PDG==130||PDG==310)
584 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
585 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
586 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
587 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
588 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
589 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
590 else if ( PDG > 3000 && PDG < 3335) ind=6;
591 else if ( PDG > -3335 && PDG < -2000) ind=7;
593 G4cout<<
"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG
594 <<
", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<
G4endl;
614 if(PDG==130||PDG==310)
619 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
620 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
621 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
622 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
623 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
624 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
625 else if ( PDG > 3000 && PDG < 3335) ind=6;
626 else if ( PDG > -3335 && PDG < -2000) ind=7;
628 G4cout<<
"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG
629 <<
", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<
G4endl;
634 if(p<=mip || p>=map)
return CalcElTot(p,ind);
637 if(nDB)
for (i=0; i<nDB; i++)
if(ind==
vItot[i])
645 lastXtot =
new std::pair<G4double,G4double>[mlp];
647 lastKtot =
static_cast<int>((lp-lpi)/dlp)+1;
676 lastKtot =
static_cast<int>((lp-lpi)/dlp)+1;
698 G4int n=
static_cast<int>(dlpp/dlp);
717 G4cout<<
"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<
",N="<<N<<
", return zero"<<
G4endl;
718 return std::make_pair(0.,0.);
720 std::pair<G4double,G4double> hp=
FetchElTot(pGeV, hPDG,
true);
721 std::pair<G4double,G4double> hn=
FetchElTot(pGeV, hPDG,
false);
723 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
735 G4cout<<
"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<
",N="<<N<<
", return zero"<<
G4endl;
736 return std::make_pair(resP,resN);
741 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N);
742 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
743 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
752 mult=1./(1.+
G4Log(pGeV+pGeV))/pGeV;
757 std::pair<G4double,G4double> hp=
FetchElTot(pGeV, hPDG,
true);
758 resP=pf*(hp.second/hp.first-1.)*mult;
762 std::pair<G4double,G4double> hn=
FetchElTot(pGeV, hPDG,
false);
763 resN=nf*(hn.second/hn.first-1.)*mult;
765 return std::make_pair(resP,resN);
786 if(NPDG==2212||NPDG==90001000)
792 else if(NPDG==90001001)
798 else if(NPDG==90002001)
804 else if(NPDG==90001002)
810 else if(NPDG==90002002)
816 else if(NPDG!=2112&&NPDG!=90000001)
818 G4cout<<
"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<
" is not 2212 or 2112"<<
G4endl;
824 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
832 if(pPDG>3400 || pPDG<-3400)
G4cout<<
"-Warning-G4QE::Scatter: pPDG="<<pPDG<<
G4endl;
834 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
839 if (PDG==2212) PDG=2112;
840 else if(PDG==2112) PDG=2212;
857 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
859 if (cost>1.) cost=1.;
860 else if(cost<-1.) cost=-1.;
866 G4cerr<<
"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<
",-t="<<mint<<
",tm="<<tm<<
G4endl;
872 if(!
RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
874 G4cerr<<
"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<
",mT="<<mT<<
",mP="<<std::sqrt(mP2)<<
G4endl;
902 if(pPDG==-211) sPDG=111;
908 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321;
909 else if(pPDG==3112) sPDG=3212;
910 else if(pPDG==3212) sPDG=3222;
911 else if(pPDG==3312) sPDG=3322;
915 if(pPDG==211) sPDG=111;
921 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321;
922 else if(pPDG==3222) sPDG=3212;
923 else if(pPDG==3212) sPDG=3112;
924 else if(pPDG==3322) sPDG=3312;
928 G4cout<<
"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<
" is not 2212 or 2112"<<
G4endl;
935 G4cout<<
"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<
", NPDG="<<NPDG<<
G4endl;
941 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
950 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
955 if (PDG==2212) PDG=2112;
956 else if(PDG==2112) PDG=2212;
973 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
975 if (cost>1.) cost=1.;
976 else if(cost<-1.) cost=-1.;
979 G4cerr<<
"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<
",t="<<mint<<
",tm="<<maxt<<
G4endl;
986 if(!
RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
988 G4cerr<<
"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<
",mT="<<mT<<
",mP="<<mS<<
G4endl;
1001 if(A<1.5)
return 0.;
1003 if (pPDG==2212) Cex=N/(A+Z);
1004 else if(pPDG==2112) Cex=Z/(A+N);
1005 else G4cout<<
"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<
G4endl;
1012 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
1013 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p;
1032 G4cerr<<
"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<
", E-p="<<dE-vP<<
G4endl;
1037 G4cerr<<
"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<
G4endl;
1038 theMomentum.setE(vP+emodif+.01*accuracy);
1049 if(vdir.mag2() > 0.)
1056 if(maxCost> 1.) maxCost= 1.;
1057 if(minCost<-1.) minCost=-1.;
1058 if(maxCost<-1.) maxCost=-1.;
1059 if(minCost> 1.) minCost= 1.;
1060 if(minCost> maxCost) minCost=maxCost;
1061 if(std::fabs(iM-fM-sM)<.00000001)
1065 f4Mom=fR*theMomentum;
1066 s4Mom=sR*theMomentum;
1069 else if (iM+.001<fM+sM || iM==0.)
1071 G4cerr<<
"***G4QH::RelDecIn2: fM="<<fM<<
"+sM="<<sM<<
">iM="<<iM<<
",d="<<iM-fM-sM<<
G4endl;
1075 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;
1089 if(std::fabs(ct)<1.) pss = p * std::sqrt(1.-ct*ct);
1095 G4ThreeVector pVect=(pss*std::sin(phi))*vz+(pss*std::cos(phi))*vy+p*ct*vx;
1097 f4Mom.setVect(pVect);
1098 f4Mom.setE(std::sqrt(fM2+p2));
1099 s4Mom.setVect((-1)*pVect);
1100 s4Mom.setE(std::sqrt(sM2+p2));
1102 if(f4Mom.e()+.001<f4Mom.rho())
G4cerr<<
"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<
",e-p="
1103 <<f4Mom.e()-f4Mom.rho()<<
G4endl;
1105 if(s4Mom.e()+.001<s4Mom.rho())
G4cerr<<
"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<
",e-p="
1106 <<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 G4Pow * GetInstance()
std::vector< G4double > vMRatio
G4double powA(G4double A, G4double y) const
static const double megaelectronvolt
CLHEP::Hep3Vector G4ThreeVector
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
std::vector< G4double > vHRatio
std::vector< G4double * > * vT
G4ChipsProtonElasticXS * PCSmanager
std::vector< G4int > vKtot
static const G4double rp2
std::pair< G4double, G4double > FetchElTot(G4double pGeV, G4int PDG, G4bool F)
std::vector< G4double > vMtot
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 > lastRtot
std::pair< G4double, G4double > GetChExFactor(G4double pIU, G4int pPDG, G4int Z, G4int N)
std::vector< G4int > vKRatio
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
std::vector< std::pair< G4double, G4double > * > * vX
std::vector< G4double * > * vL
static const double twopi
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()
std::vector< G4int > vItot
static G4Deuteron * Deuteron()
std::pair< G4double, G4double > * lastXtot
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)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
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::vector< G4int > vARatio
std::pair< G4double, G4double > GetElTotXS(G4double Mom, G4int PDG, G4bool F)
std::vector< G4int > vNRatio
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