58 const G4int mps=nps+1;
62 const G4int mls=nls+1;
77 const G4int mlp=nlp+1;
89 vT =
new std::vector<G4double*>;
90 vL =
new std::vector<G4double*>;
91 vX =
new std::vector<std::pair<G4double,G4double>*>;
108 lastXtot =
new std::pair<G4double,G4double>;
117 std::vector<G4double*>::iterator
pos;
118 for(pos=
vT->begin(); pos<
vT->end(); pos++)
123 for(pos=
vL->begin(); pos<
vL->end(); pos++)
128 std::vector<std::pair<G4double,G4double>*>::iterator pos2;
129 for(pos2=
vX->begin(); pos2<
vX->end(); pos2++)
143 if(tgA<2)
return std::make_pair(QF2In,R);
144 std::pair<G4double,G4double> ElTot=
GetElTot(pIU, pPDG, tgZ, tgN);
146 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.;
147 else if(ElTot.second>0.)
149 R=ElTot.first/ElTot.second;
152 return std::make_pair(QF2In,R);
159 if(m_s<toler || A<2)
return 1.;
160 if(m_s>min_s)
return 0.;
171 if(nDB)
for (i=0; i<nDB; i++)
if(A==
vARatio[i])
198 lastKRatio =
static_cast<int>((ls-lsi)/dls)+1;
267 lastKRatio =
static_cast<int>((ls-lsi)/dls)+1;
290 G4int n=
static_cast<int>(m_s/ds);
298 G4int n=
static_cast<int>(ls/dls);
313 G4double ss=std::sqrt(std::sqrt(m_s));
314 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
329 G4cout<<
"-Warning-G4QuasiElRatios::CalcElTot: p="<<p<<
" is zero or negative"<<
G4endl;
330 return std::make_pair(El,To);
337 El=1./(.00012+p2*.2);
354 El=LE+(pbe*lp2+6.72+32.6/
p)/(1.+rp2/p);
355 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
363 El=1./(.00012+p2*(.051+.1*p2));
380 El=LE+(pbe*lp2+6.72+30./
p)/(1.+.49*rp2/p);
381 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
390 El=1.53/(lr*lr+.0676);
398 El=pbe*ld2+2.4+7./
sp;
399 To=pbt*ld2+22.3+12./
sp;
414 El=LE+(pbe*ld2+2.4+7./
sp)/(1.+.7/p4)+.6/md+.05/hd;
415 To=LE*3+(pbt*ld2+22.3+12./
sp)/(1.+.4/p4)+1./md+.06/hd;
425 El=13./(lr2+lr2*lr2+.0676);
433 El=pbe*ld2+2.4+6./
sp;
434 To=pbt*ld2+22.3+5./
sp;
448 El=LE+(pbe*ld2+2.4+6./
sp)/(1.+3./p4)+.7/md;
449 To=LE+(pbt*ld2+22.3+5./
sp)/(1.+1./p4)+.8/md;
480 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
481 To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
491 El=.7/(lr*lr+.0676)+2./md;
512 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
513 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
521 El=1./(.002+p2*(.12+p2));
529 El=(pbe*lp2+6.72)/(1.+2./
sp);
530 To=(pbt*lp2+38.2+900./
sp)/(1.+27./sp);
540 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
541 To=LE+(pbt*lp2+38.2+900./
sp)/(1.+27./sp+3./p4);
559 El=80./(ye+1.)+pbe*lp2+6.72;
560 To=(80./yt+.3)/yt+pbt*lp2+38.2;
565 G4cout<<
"*Error*G4QuasiElRatios::CalcElTot:ind="<<I<<
" is not defined (0-7)"<<
G4endl;
569 return std::make_pair(El,To);
578 if(PDG==130||PDG==310)
583 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
584 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
585 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
586 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
587 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
588 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
589 else if ( PDG > 3000 && PDG < 3335) ind=6;
590 else if ( PDG > -3335 && PDG < -2000) ind=7;
592 G4cout<<
"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG
593 <<
", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<
G4endl;
613 if(PDG==130||PDG==310)
618 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
619 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
620 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
621 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
622 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
623 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
624 else if ( PDG > 3000 && PDG < 3335) ind=6;
625 else if ( PDG > -3335 && PDG < -2000) ind=7;
627 G4cout<<
"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG
628 <<
", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<
G4endl;
636 if(nDB)
for (i=0; i<nDB; i++)
if(ind==
vItot[i])
644 lastXtot =
new std::pair<G4double,G4double>[mlp];
646 lastKtot =
static_cast<int>((lp-lpi)/dlp)+1;
675 lastKtot =
static_cast<int>((lp-lpi)/dlp)+1;
697 G4int n=
static_cast<int>(dlpp/dlp);
716 G4cout<<
"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<
",N="<<N<<
", return zero"<<
G4endl;
717 return std::make_pair(0.,0.);
719 std::pair<G4double,G4double> hp=
FetchElTot(pGeV, hPDG,
true);
720 std::pair<G4double,G4double> hn=
FetchElTot(pGeV, hPDG,
false);
722 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
734 G4cout<<
"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<
",N="<<N<<
", return zero"<<
G4endl;
735 return std::make_pair(resP,resN);
740 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+
N);
741 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+
Z);
742 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
751 mult=1./(1.+
G4Log(pGeV+pGeV))/pGeV;
756 std::pair<G4double,G4double> hp=
FetchElTot(pGeV, hPDG,
true);
757 resP=pf*(hp.second/hp.first-1.)*mult;
761 std::pair<G4double,G4double> hn=
FetchElTot(pGeV, hPDG,
false);
762 resN=nf*(hn.second/hn.first-1.)*mult;
764 return std::make_pair(resP,resN);
785 if(NPDG==2212||NPDG==90001000)
791 else if(NPDG==90001001)
797 else if(NPDG==90002001)
803 else if(NPDG==90001002)
809 else if(NPDG==90002002)
815 else if(NPDG!=2112&&NPDG!=90000001)
817 G4cout<<
"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<
" is not 2212 or 2112"<<
G4endl;
831 if(pPDG>3400 || pPDG<-3400)
G4cout<<
"-Warning-G4QE::Scatter: pPDG="<<pPDG<<
G4endl;
833 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
838 if (PDG==2212) PDG=2112;
839 else if(PDG==2112) PDG=2212;
856 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
858 if (cost>1.) cost=1.;
859 else if(cost<-1.) cost=-1.;
865 G4cerr<<
"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<
",-t="<<mint<<
",tm="<<tm<<
G4endl;
871 if(!
RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
873 G4cerr<<
"G4QFR::Scat:t="<<tot4M<<tot4M.
m()<<
",mT="<<mT<<
",mP="<<std::sqrt(mP2)<<
G4endl;
901 if(pPDG==-211) sPDG=111;
907 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321;
908 else if(pPDG==3112) sPDG=3212;
909 else if(pPDG==3212) sPDG=3222;
910 else if(pPDG==3312) sPDG=3322;
914 if(pPDG==211) sPDG=111;
920 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321;
921 else if(pPDG==3222) sPDG=3212;
922 else if(pPDG==3212) sPDG=3112;
923 else if(pPDG==3322) sPDG=3312;
927 G4cout<<
"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<
" is not 2212 or 2112"<<
G4endl;
934 G4cout<<
"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<
", NPDG="<<NPDG<<
G4endl;
949 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
954 if (PDG==2212) PDG=2112;
955 else if(PDG==2112) PDG=2212;
972 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
974 if (cost>1.) cost=1.;
975 else if(cost<-1.) cost=-1.;
978 G4cerr<<
"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<
",t="<<mint<<
",tm="<<maxt<<
G4endl;
985 if(!
RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
987 G4cerr<<
"G4QFR::ChEx:t="<<tot4M<<tot4M.
m()<<
",mT="<<mT<<
",mP="<<mS<<
G4endl;
1000 if(A<1.5)
return 0.;
1002 if (pPDG==2212) Cex=N/(A+
Z);
1003 else if(pPDG==2112) Cex=Z/(A+
N);
1004 else G4cout<<
"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<
G4endl;
1011 G4double T=(6.75+.14*dl1*dl1+13./
p)/(1.+.14/p4)+.6/(p4+.00013);
1012 G4double U=(6.25+8.33e-5/p4/
p)*(p*sp+.34)/p2/
p;
1031 G4cerr<<
"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<
", E-p="<<dE-vP<<
G4endl;
1036 G4cerr<<
"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<
G4endl;
1037 theMomentum.
setE(vP+emodif+.01*accuracy);
1048 if(vdir.
mag2() > 0.)
1055 if(maxCost> 1.) maxCost= 1.;
1056 if(minCost<-1.) minCost=-1.;
1057 if(maxCost<-1.) maxCost=-1.;
1058 if(minCost> 1.) minCost= 1.;
1059 if(minCost> maxCost) minCost=maxCost;
1060 if(std::fabs(iM-fM-sM)<.00000001)
1064 f4Mom=fR*theMomentum;
1065 s4Mom=sR*theMomentum;
1068 else if (iM+.001<fM+sM || iM==0.)
1070 G4cerr<<
"***G4QH::RelDecIn2: fM="<<fM<<
"+sM="<<sM<<
">iM="<<iM<<
",d="<<iM-fM-sM<<
G4endl;
1074 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;
1088 if(std::fabs(ct)<1.) pss = p * std::sqrt(1.-ct*ct);
1094 G4ThreeVector pVect=(pss*std::sin(phi))*vz+(pss*std::cos(phi))*vy+p*ct*vx;
1097 f4Mom.
setE(std::sqrt(fM2+p2));
1099 s4Mom.
setE(std::sqrt(sM2+p2));
1101 if(f4Mom.
e()+.001<f4Mom.
rho())
G4cerr<<
"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<
",e-p="
1104 if(s4Mom.
e()+.001<s4Mom.
rho())
G4cerr<<
"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<
",e-p="