66 theNucleons(),currentNucleon(-1),
67 rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
69 probVect[0]=mediRatio;
70 for(
G4int i=1; i<256; i++) {probVect[i] = 0.;}
72 G4cout<<
"G4QNucleus::Constructor:(1) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
73 <<
", R="<<mediRatio<<
G4endl;
78 G4QHadron(90000000+s_value*1000000+z*1000+n),
Z(z),
N(n),S(s_value), dZ(0),dN(0),dS(0), maxClust(0),
79 theNucleons(), currentNucleon(-1), rho0(1.), radius(1.),
80 Tb(), TbActive(false), RhoActive(false)
82 probVect[0]=mediRatio;
83 for(
G4int i=1; i<256; i++) {probVect[i] = 0.;}
85 G4cout<<
"G4QNucleus::Construction By Z="<<z<<
",N="<<n<<
",S="<<s_value<<
G4endl;
87 SetZNSQC(z,n,s_value);
90 G4cout<<
"G4QNucleus::ConstructionByZNS: nQPDG="<<nQPDG<<
G4endl;
94 G4cout<<
"G4QNucleus::ConstructionByZNS: mass="<<mass<<
G4endl;
98 G4cout<<
"G4QNucleus::ConstructionByZNS: nQPDG set"<<
G4endl;
104 G4cout<<
"G4QNucleus::Constructor:(2) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
105 <<
", R="<<mediRatio<<
G4endl;
110 G4QHadron(nucPDG), maxClust(0), theNucleons(),
111 currentNucleon(-1), rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
113 if(nucPDG==22) nucPDG=90000000;
118 G4cout<<
"G4QNucleus::Constructor:(3) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
119 <<
", R="<<mediRatio<<
", 4M="<<p<<
G4endl;
124 G4QHadron(nucPDG, p), maxClust(0), theNucleons(),
125 currentNucleon(-1), rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
130 G4cout<<
"G4QNucleus::Constructor:(4) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
131 <<
", R="<<mediRatio<<
", 4M="<<p<<
G4endl;
136 G4QHadron(90000000+s_value*1000000+z*1000+n,p),
Z(z),
N(n),S(s_value), dZ(0),dN(0),dS(0), maxClust(0),
137 theNucleons(), currentNucleon(-1), rho0(1.),radius(1.),
138 Tb(), TbActive(false), RhoActive(false)
140 probVect[0]=mediRatio;
141 for(
G4int i=1; i<256; i++) {probVect[i] = 0.;}
148 SetZNSQC(z,n,s_value);
150 G4cout<<
"G4QNucleus::Constructor:(5) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
151 <<
", R="<<mediRatio<<
G4endl;
156 G4QHadron(nucQC), dZ(0), dN(0), dS(0), maxClust(0), theNucleons(), currentNucleon(-1),
157 rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
161 G4cout<<
"G4QNucleus::Construction By QC="<<nucQC<<
G4endl;
163 probVect[0]=mediRatio;
164 for(
G4int i=1; i<256; i++) {probVect[i] = 0.;}
174 G4cout<<
"G4QNucleus::ConstructionByQC: N="<<N<<
",Z="<<Z<<
",S="<<S<<
G4endl;
176 G4int nucPDG=90000000+S*1000000+Z*1000+N;
179 G4cout<<
"G4QNucleus::ConstructionByQC: nQPDG="<<nQPDG<<
G4endl;
184 if(nucQC.
GetTot()) mass=mPi0;
188 G4cout<<
"G4QNucleus::ConstructionByQC: mass="<<mass<<
G4endl;
192 G4cout<<
"G4QNucleus::ConstructionByQC: nQPDG set"<<
G4endl;
198 G4cout<<
"G4QNucleus::Constructor:(6) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
199 <<
", R="<<mediRatio<<
G4endl;
204 G4QHadron(nucQC,p), dZ(0), dN(0), dS(0), maxClust(0), theNucleons(), currentNucleon(-1),
205 rho0(1.), radius(1.), Tb(), TbActive(false), RhoActive(false)
208 G4cout<<
"G4QNucleus::(LV)Construction By QC="<<nucQC<<
G4endl;
210 probVect[0]=mediRatio;
211 for(
G4int i=1; i<256; i++) {probVect[i] = 0.;}
222 G4cout<<
"G4QNucleus::(LV)ConstructionByQC: N="<<N<<
",Z="<<Z<<
",S="<<S<<
G4endl;
228 G4cout<<
"G4QNucleus::Constructor:(7) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
229 <<
", R="<<mediRatio<<
G4endl;
241 maxClust = right->maxClust;
242 for(
G4int i=0; i<=maxClust; i++) probVect[i] = right->probVect[i];
243 probVect[254] = right->probVect[254];
244 probVect[255] = right->probVect[255];
246 TbActive = right->TbActive;
247 RhoActive = right->RhoActive;
253 radius = right->radius;
256 G4int nn=right->theNucleons.size();
260 theNucleons.push_back(nucleon);
264 G4cout<<
"G4QNucleus::Constructor:(8) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
265 <<
", R="<<mediRatio<<
G4endl;
278 maxClust = right.maxClust;
279 for(
G4int i=0; i<=maxClust; i++) probVect[i] = right.probVect[i];
280 probVect[254] = right.probVect[254];
281 probVect[255] = right.probVect[255];
283 TbActive = right.TbActive;
284 RhoActive = right.RhoActive;
290 radius = right.radius;
293 G4int nn=right.theNucleons.size();
297 theNucleons.push_back(nucleon);
301 G4cout<<
"G4QNucleus::Constructor:(9) N="<<freeNuc<<
", D="<<freeDib<<
", W="<<clustProb
302 <<
", R="<<mediRatio<<
G4endl;
312 TbActive = right.TbActive;
314 RhoActive = right.RhoActive;
316 radius = right.radius;
317 G4int nn = right.theNucleons.size();
321 theNucleons.push_back(nucleon);
333 maxClust = right.maxClust;
334 for(
G4int i=0; i<=maxClust; i++) probVect[i] = right.probVect[i];
335 probVect[254] = right.probVect[254];
336 probVect[255] = right.probVect[255];
343 for_each(theNucleons.begin(),theNucleons.end(),
DeleteQHadron());
366 lhs<<
"{Z="<<rhs.
GetZ()<<
",N="<<rhs.
GetN()<<
",S="<<rhs.
GetS()<<
"}";
373 static const G4int NUCPDG = 90000000;
375 G4cout<<
"G4QNucleus::InitByPDG: >Called< PDG="<<nucPDG<<
G4endl;
380 probVect[0]=mediRatio;
381 for(
G4int i=1; i<256; i++) {probVect[i] = 0.;}
387 if(nucPDG>80000000 && nucPDG<100000000)
394 G4cout<<
"G4QNucleus::InitByPDG:Z="<<Z<<
",N="<<N<<
",S="<<S<<
G4endl;
399 if(nucPDG!=NUCPDG) PDGMass=nPDG.
GetMass();
405 G4cout<<
"G4QNucleus::InitByPDG:->QPDG="<<nPDG<<
": 4M="<<p<<
G4endl;
410 G4cerr<<
"***G4QNucleus::InitByPDG:Initialized by not nuclear PDGCode="<<nucPDG<<
G4endl;
428 probVect[0]=mediRatio;
436 G4cout<<
"G4QN::UpdateCl:A="<<a<<
"(Z="<<Z<<
",N="<<N<<
",S="<<S<<
"),mR="<<mediRatio<<
G4endl;
442 G4cout<<
"***G4QNucleus::UpdateClusters:No clusters can be calculated as A="<<A<<
G4endl;
451 G4cout<<
"G4QN::UpdateCl:surf="<<surf<<
"= N="<<freeNuc<<
"+D="<<freeDib<<
",A="<<sA<<
G4endl;
454 if (din && dA<2 && a>2)
460 G4cout<<
"G4QN::UpdtC:dA="<<dA<<
",A="<<A<<
",s="<<surf<<
",S="<<sA<<
",C="<<maxClust<<
G4endl;
467 pA=0.5*freeDib*sA/surf;
477 probVect[1]= uA+dA/A;
487 sum+= probVect[2]+probVect[2];
493 G4cout<<
"G4QNucleus::UpdateClust:Only quasi-free nucleons pV[1]="<<probVect[1]<<
G4endl;
508 G4cout<<
"G4QNucl::UpdateCl:sud="<<sud<<
",v[1]=s="<<sum<<
",dA="<<dA<<
",uA="<<uA<<
G4endl;
522 G4cout<<
"G4QNucl::UpdateCl:sud="<<sud<<
",v[2]="<<prb<<
",s="<<sum<<
",pA="<<pA<<
G4endl;
532 G4cout<<
"G4QNucleus::UpdateClusters:p1="<<probVect[1]<<
", p2="<<probVect[2]<<
",sA="<<sA
533 <<
",uA="<<uA<<
",pA="<<pA<<
",wrd="<<wrd<<
",sud="<<sud<<
G4endl;
540 if(maxClust<dA) dLim=maxClust;
541 for (
int i=3; i<=dLim; i++)
546 G4cout<<
"G4QNucleus::UpdateCl:sud="<<sud<<
", v["<<i<<
"]="<<rd<<
", s="<<sum<<
G4endl;
554 G4cout<<
"G4QNucleus::UpdateCl:Cluster of "<<i<<
" baryons,pV="<<probVect[i]<<
G4endl;
559 dZ =
static_cast<int>(
static_cast<double>((dA-dS)*Z)/(Z+N) + 0.5);
563 G4cout<<
"G4QNucleus::UpdateClusters: Sum of weighted probabilities s="<<sum<<
G4endl;
615 G4QHadronVector::iterator u;
616 for(u=theNucleons.begin(); u!=theNucleons.end(); u++)
619 G4cout<<
"G4QNucleus::SubtractNucleon: LOOP 4M="<<(*u)->Get4Momentum()<<
G4endl;
628 if (NotFound)
G4Exception(
"G4QNucleus::SubtractNucleon()",
"HAD_CHPS_0000",
635 G4cout<<
"G4QNucleus::SubtractNucleon: InitialNucleus 4M="<<t4M<<
", PDG="<<tPDG<<
", nN="
636 <<theNucleons.size()<<
G4endl;
638 G4int uPDG=(*u)->GetPDGCode();
641 G4cout<<
"G4QNucleus::SubtractNucleon: subtractNucleon 4M="<<u4M<<
",PDG="<<uPDG<<
G4endl;
644 theNucleons.erase(u);
647 if (uPDG==2212) tPDG-=1000;
648 else if(uPDG==2112) tPDG--;
654 ed <<
"Impossible nucleon PDG Code: Unexpected Nucleon PDGCode ="
656 G4Exception(
"G4QNucleus::SubtractNucleon()",
"HAD_CHPS_0001",
660 G4cout<<
"G4QNucleus::SubtractNucleon: theResidualNucleus PDG="<<tPDG<<
", 4M="<<t4M
661 <<
", nN="<<theNucleons.size()<<
G4endl;
669 G4cout<<
"G4QNucleus::SubtractNucleon: rAm2="<<mR2<<
" =? 4Mm2="<<tM2<<
G4endl;
672 if(std::fabs(mR2-tM2)>.01)
G4cout<<
"*G4QNucleus::SubNucleon:rM="<<mR2<<
"#"<<tM2<<
G4endl;
677 for(u=theNucleons.begin(); u!=theNucleons.end(); u++)
682 if((*u)->GetPDGCode()==2212) m2_value=m2p;
683 G4double srE=std::sqrt(srP2+m2_value);
685 G4cout<<
"G4QNucleus::SubtractNucleon:#"<<cnt++<<
", correctedEnergy="<<tE-srE<<
G4endl;
688 (*u)->Set4Momentum(n4M);
699 G4QHadronVector::iterator u;
700 for(u=theNucleons.begin(); u!=theNucleons.end(); u++)
delete *u;
707 static const G4int NUCPDG=90000000;
708 if(cPDG>80000000&&cPDG!=NUCPDG)
711 G4int newPDG=curPDG-cPDG+NUCPDG;
724 else if(cPDG!=NUCPDG)
G4cerr<<
"***G4QN::Reduce:Subtract not nuclear PDGC="<<cPDG<<
G4endl;
731 static const G4int NUCPDG=90000000;
732 if(cPDG>80000000&&cPDG!=NUCPDG)
743 else G4cerr<<
"***G4QNucleus::Increase: PDGCode="<<cPDG<<
",4M="<<c4M<<
G4endl;
759 G4int zns=z+n+s_value;
787 if(baryn<2)
return 0;
792 G4cout<<
"G4QNucleus::SplitBaryon: B="<<baryn<<
", M="<<totM<<valQC<<
G4endl;
802 G4cout<<
"G4QNucleus::SplitBaryon: (neutron),sM="<<sM<<
",d="<<totM-sM<<
G4endl;
804 if(sM<totM+.001)
return 2112;
816 G4cout<<
"G4QNucleus::SplitBaryon: (proton),sM="<<sM<<
",d="<<totM-sM<<
G4endl;
818 if(sM<totM+.001)
return 2212;
828 G4cout<<
"G4QNucleus::SplitBaryon: (lambda),sM="<<sM<<
",d="<<totM-sM<<
G4endl;
830 if(sM<totM+.001)
return 3122;
842 G4cout<<
"G4QNucleus::SplitBaryon: (deuteron),sM="<<sM<<
",d="<<totM-sM<<
G4endl;
844 if(sM<totM+.001)
return 90001001;
853 if(NQ!=4||PQ!=4) sM+=CB;
855 G4cout<<
"G4QNucleus::SplitBaryon: (alpha),sM="<<sM<<
",d="<<totM-sM<<
G4endl;
857 if(sM<totM+.001)
return 90002002;
872 if(baryn<3)
return false;
876 G4cout<<
"G4QNucleus::Split2Baryons: B="<<baryn<<
", M="<<totM<<valQC<<
G4endl;
886 G4cout<<
"G4QNucleus::Split2Baryons: (2 neutrons), sM="<<sM<<
", d="<<totM-sM<<
G4endl;
888 if(sM<totM)
return true;
898 G4cout<<
"G4QNucleus::Split2Baryons: (2 protons), sM="<<sM<<
", d="<<totM-sM<<
G4endl;
900 if(sM<totM)
return true;
909 G4cout<<
"G4QNucleus::Split2Baryons:(proton+neutron), sM="<<sM<<
", d="<<totM-sM<<
G4endl;
911 if(sM<totM)
return true;
921 G4cout<<
"G4QNucleus::Split2Baryons:(lambda+neutron), sM="<<sM<<
", d="<<totM-sM<<
G4endl;
923 if(sM<totM)
return true;
932 G4cout<<
"G4QNucleus::Split2Baryons: (proton+lambda), sM="<<sM<<
", d="<<totM-sM<<
G4endl;
934 if(sM<totM)
return true;
943 G4cout<<
"G4QNucleus::Split2Baryons: (two lambdas), sM="<<sM<<
", d="<<totM-sM<<
G4endl;
945 if(sM<totM)
return true;
964 static const G4int gPDG = 22;
966 static const G4int nPDG = 2112;
969 static const G4int pPDG = 2212;
972 static const G4int lPDG = 3122;
977 static const G4int dPDG = 90001001;
978 static const G4int aPDG = 90002002;
998 static const G4double mN2 = mNeut*mNeut;
999 static const G4double mP2 = mProt*mProt;
1000 static const G4double mL2 = mLamb*mLamb;
1001 static const G4double mA2 = mAlph*mAlph;
1002 static const G4double mNP = mNeut+mProt;
1015 G4cout<<
"G4QNucleus::EvaporBaryon: *Called*, a="<<a<<GetThis()<<
",alph="<<evalph<<
G4endl;
1036 if(DAPBarr>SAPBarr)SAPBarr=DAPBarr;
1044 G4cout<<
"G4QN::EB:pB="<<PBarr<<
",aB="<<ABarr<<
",ppB="<<PPBarr<<
",paB="<<PABarr<<
G4endl;
1050 G4int nucPDG = -2112;
1063 if(totMass > mPi+nucM+nucM)
1070 G4cerr<<
"***G4QNucl::EvapBary: (anti) tM="<<totMass<<
"-> 2N="<<nucPDG<<
"(M="
1071 <<nucM<<
") + pi="<<piPDG<<
"(M="<<mPi<<
")"<<
G4endl;
1084 G4cerr<<
"***G4QNucleus::EvaporateBaryon: M="<<totMass
1085 <<
", M="<<totMass<<
" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<
G4endl;
1090 else if(Z==2 || N==2)
1092 G4int nucPDG = -2112;
1103 if(totMass > mPi+mPi+nucM+nucM)
1110 G4cerr<<
"***G4QNucl::EvapBary: (anti) tM="<<totMass<<
"-> 2N="<<nucPDG<<
"(M="
1111 <<nucM<<
") + 2pi="<<piPDG<<
"(M="<<mPi<<
")"<<
G4endl;
1126 G4cerr<<
"***G4QNucleus::EvaporateBaryon: M="<<totMass
1127 <<
", M="<<totMass<<
" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<
G4endl;
1138 if(!
DecayIn2(h1mom,h2mom))
return false;
1146 if(!
DecayIn2(h1mom,h2mom))
return false;
1148 else if(N==-1 && Z==-1)
1154 if(!
DecayIn2(h1mom,h2mom))
return false;
1156 else if(Z==-1 && S==-1)
1162 if(!
DecayIn2(h1mom,h2mom))
return false;
1170 if(!
DecayIn2(h1mom,h2mom))
return false;
1180 G4int nucPDG = 2112;
1193 if(totMass>mPi+nucM+nucM)
1200 G4cerr<<
"***G4QNucl::EvapBary: tM="<<totMass<<
"-> 2N="<<nucPDG<<
"(M="
1201 <<nucM<<
") + pi="<<piPDG<<
"(M="<<mPi<<
")"<<
G4endl;
1214 G4cerr<<
"***G4QNucleus::EvaporateBaryon: M="<<totMass
1215 <<
", M="<<totMass<<
" < M_2N+Pi, d="<<totMass-2*nucM-mPi<<
G4endl;
1226 if(!
DecayIn2(h1mom,h2mom))
return false;
1234 if(!
DecayIn2(h1mom,h2mom))
return false;
1241 G4cout<<
"G4QNucl::EvaporateBaryon: Photon ### d+g ###, dM="<<totMass-mNP<<
G4endl;
1255 if(!
DecayIn2(h1mom,h2mom))
return false;
1263 if(!
DecayIn2(h1mom,h2mom))
return false;
1271 if(!
DecayIn2(h1mom,h2mom))
return false;
1369 G4cout<<
"G4QNuc::EvaB:a>2, totM="<<totMass<<
" > GSMass="<<GSMass<<
",d="<<totMass-GSMass
1384 PQPDG=
G4QPDGCode(90000000+1000*(1000*S+Z-1)+N);
1388 pp2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2;
1392 pBnd=mProt-GSMass+GSResNp;
1395 G4cout<<
"G4QNuc::EvapBaryon:pm="<<eMax+sqrt(pp2m+GSResNp*GSResNp)<<
" = M="<<totMass
1396 <<
", sm="<<GSResNp+mProt+PBarr<<
",pp2="<<pp2m<<
",pB="<<pBnd<<
G4endl;
1398 pExcess=eMax-mProt+pBnd;
1403 ppQPDG=
G4QPDGCode(90000000+1000*(1000*S+Z-2)+N);
1407 G4cout<<
"G4QNucl::EvapBaryon: ppM="<<GSResPP<<
",T="<<sm-GSMass<<
",E="<<totMass-sm
1410 if(GSResPP+mProt+mProt+SPPBarr<totMass) ppFlag=
true;
1420 paQPDG =
G4QPDGCode(90000000+1000*(1000*S+Z-3)+N-2);
1423 G4double s_value=GSResPA+mAlph+mProt+SAPBarr;
1424 G4cout<<
"G4QN::EB:paM="<<GSResPA<<
",T="<<s_value-GSMass<<
",E="<<totMass-s_value<<
G4endl;
1426 if(GSResPA+mProt+SAPBarr+mAlph<totMass) paFlag=
true;
1535 aaQPDG =
G4QPDGCode(90000000+1000*(1000*S+Z-4)+N-4);
1538 G4double s_value=GSResAA+mAlph+mAlph+SAABarr;
1539 G4cout<<
"G4QNucl::EvapBaryon: a="<<GSResNP<<
",T="<<s_value-GSMass<<
",E="<<totMass-s_value
1540 <<
",A="<<SAABarr<<
G4endl;
1542 if(GSResAA+mAlph+mAlph+SAABarr<totMass) aaFlag=
true;
1546 naQPDG =
G4QPDGCode(90000000+1000*(1000*S+Z-2)+N-3);
1549 G4double s_value=GSResNA+mAlph+mNeut;
1550 G4cout<<
"G4QNucl::EvapBary: M="<<GSResNA<<
",T="<<s_value-GSMass<<
",E="<<totMass-s_value
1553 if(GSResNA+mNeut+mAlph+ABarr<totMass) naFlag=
true;
1557 laQPDG =
G4QPDGCode(90000000+1000*(1000*S+Z-1002)+N-2);
1559 if(GSResLA+mLamb+mAlph+ABarr<totMass) laFlag=
true;
1561 AQPDG =
G4QPDGCode(90000000+1000*(1000*S+Z-2)+N-2);
1565 ap2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2;
1569 aBnd=mAlph-GSMass+GSResNa;
1572 G4cout<<
"G4QNuc::EvapBar:m="<<eMax+sqrt(ap2m+GSResNa*GSResNa)<<
" = M="<<totMass
1573 <<
", sm="<<GSResNp+mProt+PBarr<<
",pp2="<<pp2m<<
",pB="<<pBnd<<
G4endl;
1575 aExcess=eMax-mAlph+aBnd;
1584 npQPDG=
G4QPDGCode(90000000+1000*(1000*S+Z-1)+N-1);
1587 G4double s_value=GSResNP+mNeut+mProt;
1588 G4cout<<
"G4QNucl::EvapBaryon: npM="<<GSResNP<<
",T="<<s_value-GSMass<<
",E="<<totMass-s_value
1591 if(GSResNP+mNeut+mProt+PBarr<totMass) npFlag=
true;
1614 plQPDG=
G4QPDGCode(90000000+1000*(1000*(S-1)+Z-1)+N);
1616 if(GSResPL+mProt+PBarr+mLamb<totMass) plFlag=
true;
1634 NQPDG=
G4QPDGCode(90000000+1000*(1000*S+Z)+N-1);
1637 G4cout<<
"G4QNucleus::EvapBaryon: M(A-N)="<<GSResNn<<
",Z="<<Z
1638 <<
",N="<<N<<
",S="<<S<<
G4endl;
1642 np2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2;
1646 nBnd=mNeut-GSMass+GSResNn;
1649 G4cout<<
"G4QNuc::EvapBaryon:nm="<<eMax+sqrt(np2m+GSResNn*GSResNn)<<
" = M="<<totMass
1650 <<
", sm="<<GSResNn+mNeut<<
",np2="<<np2m<<
",nB="<<nBnd<<
G4endl;
1652 nExcess=eMax-mNeut+nBnd;
1657 nnQPDG=
G4QPDGCode(90000000+1000*(1000*S+Z)+N-2);
1659 if(GSResNN+mNeut+mNeut<totMass) nnFlag=
true;
1679 nlQPDG=
G4QPDGCode(90000000+1000*(1000*(S-1)+Z)+N-1);
1681 if(GSResNL+mNeut+mLamb<totMass) nlFlag=
true;
1698 LQPDG=
G4QPDGCode(90000000+1000*(1000*(S-1)+Z)+N);
1702 lp2m=(tM2-mpls*mpls)*(tM2-mmin*mmin)/qtM2;
1706 lBnd=mLamb-GSMass+GSResNl;
1709 G4cout<<
"G4QNuc::EvapBaryon:lm="<<eMax+sqrt(lp2m+GSResNl*GSResNl)<<
" = M="<<totMass
1710 <<
", sm="<<GSResNl+mLamb<<
",lp2="<<lp2m<<
",lB="<<lBnd<<
G4endl;
1712 lExcess=eMax-mLamb+lBnd;
1717 llQPDG=
G4QPDGCode(90000000+1000*(1000*(S-2)+Z)+N);
1719 if(GSResLL+mLamb+mLamb<totMass) llFlag=
true;
1730 G4bool nSecF = nnFlag||npFlag||nlFlag||naFlag;
1731 G4bool pSecF = npFlag||ppFlag||plFlag||paFlag;
1732 G4bool lSecF = nlFlag||plFlag||llFlag||laFlag;
1733 G4bool aSecF = naFlag||paFlag||laFlag||aaFlag;
1738 G4bool secB = nSecF||pSecF||lSecF||aSecF;
1741 G4cout<<
"G4QNucl::EvapBary:n="<<nSecF<<
",p="<<pSecF<<
",l="<<lSecF<<
",a="<<aSecF<<
",nn="
1742 <<nnFlag<<
", np="<<npFlag<<
",pp="<<ppFlag<<
",pa="<<paFlag<<
",na="<<naFlag<<
",aa="
1750 if(!nSecF) nFlag=
false;
1751 if(!pSecF) pFlag=
false;
1752 if(!lSecF) lFlag=
false;
1753 if(!aSecF) aFlag=
false;
1755 G4cout<<
"G4QNuc::EB:nF="<<nFlag<<
",pF="<<pFlag<<
",lF="<<lFlag<<
",aF="<<aFlag<<
G4endl;
1758 if(nFlag&&nExcess>maxE) maxE=nExcess;
1759 if(pFlag&&pExcess>maxE) maxE=pExcess;
1760 if(lFlag&&lExcess>maxE) maxE=lExcess;
1761 if(lFlag&&aExcess>maxE) maxE=aExcess;
1763 if(pFlag)pMin+= PBarr;
1767 if(aFlag)aMin+= ABarr;
1769 if(nFlag&&nMin<minE) minE=nMin;
1770 if(pFlag&&pMin<minE) minE=pMin;
1771 if(lFlag&&lMin<minE) minE=lMin;
1772 if(evalph&&aFlag&&aMin<minE) minE=aMin;
1775 G4cout<<
"G4QNucleus::EvapBaryon: nE="<<nExcess<<
">"<<nMin<<
",pE="<<pExcess<<
">"<<pMin
1776 <<
",sE="<<lExcess<<
">"<<lMin<<
",E="<<aExcess<<
">"<<aMin<<
",miE="<<minE<<
"<maE="
1786 if( ( (pFlag && pExcess > pMin) ||
1787 (nFlag && nExcess > nMin) ||
1788 (lFlag && lExcess > lMin) ||
1789 (aFlag && aExcess > aMin) ) && minE<maxE )
1807 G4cout<<
"G4QNuc::EvapBary:iE="<<minE<<
",aE="<<maxE<<
",mi="<<mi<<
",mm="<<mm_value<<
",ma="
1820 G4cerr<<
"***G4QNucleus::EvapBaryon: M="<<mm_value/ma<<
",xi="<<xMi<<
",xa="<<xMa<<
G4endl;
1826 G4cout<<
"G4QNuc:EvapBaryon:mi="<<mi<<
",ma="<<ma<<
", xi="<<xMi<<
",xa="<<xMa<<
G4endl;
1831 G4cout<<
"G4QNucleus::EvaporateBaryon: Power="<<powr<<
",RevPower="<<revP<<
G4endl;
1833 G4double minR=pow(1.-xMa*xMa,powr);
1834 G4double maxR=pow(1.-xMi*xMi,powr);
1836 G4cout<<
"G4QNucleus::EvaporateBaryon: miR="<<minR<<
", maR="<<maxR<<
G4endl;
1841 while(cond&&cntr<cntm)
1850 G4cerr<<
"**G4QNucl::EvapB:R="<<R<<
",xi="<<xMi<<
" < "<<x<<
" < xa="<<xMa<<
G4endl;
1864 G4cout<<
"G4QNuc::EvapB:t="<<tk<<
",pM="<<pMin<<
",pB="<<pBnd<<
",n="<<nMin<<
",a="
1872 G4cout<<
"G4QN::EB:Proton="<<kin<<
",CB="<<PBarr<<
",B="<<pBnd<<
",M="<<pMin
1900 if(evalph&&aFlag&&tk>aMin)
1905 G4cout<<
"G4QN::EB:Alpha="<<kin<<
",CB="<<ABarr<<
",p="
1908 psum+=
CoulBarPenProb(ABarr,kin,2,4)*sqrt(kin)*evalph*Z*(Z-1)*N*(N-1)
1913 G4cout<<
"G4QNuc::EvapB:"<<r<<
",p="<<zCBPP<<
",pn="<<nCBPP<<
",pnl="<<lCBPP<<
",t="
1920 G4cout<<
"G4QNuc::EvaB:ALPHA is selected for evap, r="<<r<<
">"<<lCBPP<<
G4endl;
1924 else if(r&&r>nCBPP&&r<=lCBPP)
1927 G4cout<<
"G4QNuc::EvaB:LAMBDA is selected for evap,r="<<r<<
"<"<<lCBPP<<
G4endl;
1931 else if(r&&r>zCBPP&&r<=nCBPP)
1934 G4cout<<
"G4QNuc::EvaBar: N is selected for evapor,r="<<r<<
"<"<<nCBPP<<
G4endl;
1938 else if(r&&r<=zCBPP)
1941 G4cout<<
"G4QNuc::EvaBar: P is selected for evapor,r="<<r<<
"<"<<zCBPP<<
G4endl;
1948 G4cout<<
"G4QNuc::EvapBar:c="<<cond<<
",x="<<x<<
",cnt="<<cntr<<
",R="<<R<<
",ma="<<ma
1949 <<
",rn="<<rn<<
"<r="<<x/xMa<<
",tk="<<tk<<
",ni="<<nMin<<
",pi="<<pMin<<
G4endl;
1987 G4cout<<
"G4QNucleus::EvaporateBaryon:np2="<<p2<<
",np2m="<<np2m<<
G4endl;
2011 else G4cerr<<
"***G4QNucleus::EvaporateBaryon: PDG="<<PDG<<
G4endl;
2014 if (rEn2 > p2) rMass=sqrt(rEn2-p2);
2018 G4bool nnCond = !nnFlag || (nnFlag && GSResNN+mNeut > rMass);
2019 G4bool npCond = !npFlag || (npFlag && GSResNP+mProt+PBarr > rMass);
2020 G4bool nlCond = !nlFlag || (nlFlag && GSResNL+mLamb > rMass);
2021 G4bool naCond = !naFlag || (naFlag && GSResNA+mAlph+ABarr > rMass);
2022 G4bool pnCond = !npFlag || (npFlag && GSResNP+mNeut > rMass);
2023 if(barf) pnCond = !npFlag || (npFlag && GSResNP+mNeut+PBarr > rMass);
2024 G4bool ppCond = !ppFlag || (ppFlag && GSResPP+mProt+PPBarr > rMass);
2025 if(barf) ppCond = !ppFlag || (ppFlag && GSResPP+mProt+SPPBarr > rMass);
2026 G4bool plCond = !plFlag || (plFlag && GSResPL+mLamb > rMass);
2027 if(barf) plCond = !plFlag || (plFlag && GSResPL+mLamb+PBarr > rMass);
2028 G4bool paCond = !paFlag || (paFlag && GSResPA+mAlph+APBarr > rMass);
2029 if(barf) paCond = !paFlag || (paFlag && GSResPA+mAlph+SAPBarr > rMass);
2030 G4bool lnCond = !nlFlag || (nlFlag && GSResNL+mNeut > rMass);
2031 G4bool lpCond = !plFlag || (plFlag && GSResPL+mProt+PBarr > rMass);
2032 G4bool llCond = !llFlag || (llFlag && GSResLL+mLamb > rMass);
2033 G4bool laCond = !laFlag || (laFlag && GSResLA+mAlph+ABarr > rMass);
2034 G4bool anCond = !naFlag || (naFlag && GSResNA+mNeut > rMass);
2035 if(barf) anCond = !naFlag || (naFlag && GSResNA+mNeut+ABarr > rMass);
2036 G4bool apCond = !paFlag || (paFlag && GSResPA+mProt+PABarr > rMass);
2037 if(barf) apCond = !paFlag || (paFlag && GSResPA+mProt+SAPBarr > rMass);
2038 G4bool alCond = !laFlag || (laFlag && GSResLA+mLamb > rMass);
2039 if(barf) alCond = !laFlag || (laFlag && GSResLA+mLamb+ABarr > rMass);
2040 G4bool aaCond = !aaFlag || (aaFlag && GSResAA+mAlph+AABarr > rMass);
2041 if(barf) aaCond = !aaFlag || (aaFlag && GSResAA+mAlph+SAABarr > rMass);
2043 G4cout<<
"G4QNucl::EvaB:"<<PDG<<
", E="<<tk<<
", rM="<<rMass<<
", ";
2044 if(PDG==pPDG)
G4cout<<
"PN="<<GSResNP+mNeut<<
"("<<pnCond<<
"),PP="
2045 <<GSResPP+mProt+PPBarr<<
"("<<ppCond<<
"),PL="
2046 <<GSResPL+mLamb<<
"("<<plCond<<
"),PA="
2047 <<GSResPA+mAlph+APBarr<<
"("<<paCond;
2048 else if(PDG==nPDG)
G4cout<<
"NN="<<GSResNN+mNeut<<
"("<<nnCond<<
"),NP="
2049 <<GSResNP+mProt+PBarr<<
"("<<npCond<<
"),NL="
2050 <<GSResNL+mLamb<<
"("<<nlCond<<
"),NA="
2051 <<GSResNA+mAlph+ABarr<<
"("<<naCond;
2052 else if(PDG==lPDG)
G4cout<<
"LN="<<GSResNL+mNeut<<
"("<<lnCond<<
"),LP="
2053 <<GSResPL+mProt+PBarr<<
"("<<lpCond<<
"),LL="
2054 <<GSResLL+mLamb<<
"("<<llCond<<
"),LA="
2055 <<GSResLA+mAlph+ABarr<<
"("<<laCond;
2056 else if(PDG==aPDG)
G4cout<<
"AN="<<GSResNA+mNeut<<
"("<<anCond<<
"),AP="
2057 <<GSResPA+mProt+PABarr<<
"("<<apCond<<
"),AL="
2058 <<GSResLA+mLamb<<
"("<<alCond<<
"),AA="
2059 <<GSResAA+mAlph+AABarr<<
"("<<aaCond;
2065 if(PDG==pPDG&&(pnCond&&ppCond&&plCond&&paCond))
2068 G4cout<<
"G4QN::EB:*p*: n="<<pnCond<<
",p="<<ppCond<<
",l="<<plCond<<
",a="<<paCond
2074 if(N&&GSResNP!=GSMass&&fMass+PBarr+mNeut+GSResNP<totMass)
2076 if(barf) nLim+=(N+N)*pow(totMass-mNeut-mProt-PBarr-GSResNP,2);
2077 else nLim+=(N+N)*pow(totMass-mNeut-mProt-GSResNP,2);
2080 if(Z>1&&GSResPP!=GSMass&&fMass+mProt+SPPBarr+GSResPP<totMass)
2082 if(barf) zLim+=(Z-1)*pow(totMass-mProt-mProt-SPPBarr-GSResPP,2);
2083 else zLim+=(Z-1)*pow(totMass-mProt-mProt-GSResPP,2);
2086 if(S&&GSResPL!=GSMass&&fMass+PBarr+mLamb+GSResPL<totMass)
2088 if(barf) sLim+=(S+S)*pow(totMass-mProt-mLamb-PBarr-GSResPL,2);
2089 else sLim+=(S+S)*pow(totMass-mProt-mLamb-GSResPL,2);
2092 if(evalph&&Z>2&&N>1&&a>4&&GSResPL!=GSMass&&fMass+SAPBarr+mAlph+GSResPA<totMass)
2094 if(barf) aLim+=pow(totMass-mProt-mAlph-SAPBarr-GSResPA,2)*evalph*
2095 (Z-1)*(Z-2)*N*(N-1)*12/(a-2)/(a-3)/(a-4);
2096 else aLim+=pow(totMass-mProt-mAlph-GSResPA,2)*evalph*(Z-1)*(Z-2)*N*(N-1)
2097 *12/(a-2)/(a-3)/(a-4);
2101 G4cout<<
"G4QNuc::EvaB:p, r="<<r<<
",n="<<nLim<<
",z="<<zLim<<
",s="<<sLim<<
",a="
2105 if(!aLim) three=
false;
2116 else if(zLim<sLim&&r>zLim&&r<=sLim)
2126 else if(nLim<zLim&&r>nLim&&r<=zLim)
2148 else if(PDG==nPDG&&(nnCond&&npCond&&nlCond&&naCond))
2151 G4cout<<
"G4QN::EB:*n*: n="<<nnCond<<
",p="<<npCond<<
",l="<<nlCond<<
",a="<<naCond
2157 if(N>1&&GSResNN!=GSMass&&fMass+mNeut+GSResNN<totMass)
2158 nLim+=(N-1)*pow(totMass-mNeut-mNeut-GSResNN,2);
2160 if(Z&&GSResNP!=GSMass&&fMass+mProt+PBarr+GSResNP<totMass)
2162 if(barf) zLim+=(Z+Z)*pow(totMass-mNeut-mProt-PBarr-GSResNP,2);
2163 else zLim+=(Z+Z)*pow(totMass-mNeut-mProt-GSResNP,2);
2166 if(S&&GSResNL!=GSMass&&fMass+mLamb+GSResNL<totMass)
2167 sLim+=(S+S)*pow(totMass-mNeut-mLamb-GSResNL,2);
2169 if(evalph&&Z>1&&N>2&&a>4&&GSResNA!=GSMass&&fMass+mAlph+ABarr+GSResNA<totMass)
2171 if(barf) aLim+=pow(totMass-mNeut-mAlph-ABarr-GSResNA,2)*
2172 evalph*Z*(Z-1)*(N-1)*(N-2)*12/(a-2)/(a-3)/(a-4);
2173 else aLim+=pow(totMass-mNeut-mAlph-GSResNA,2)*
2174 evalph*Z*(Z-1)*(N-1)*(N-2)*12/(a-2)/(a-3)/(a-4);
2178 G4cout<<
"G4QN::EB:n, r="<<r<<
",n="<<nLim<<
",z="<<zLim<<
",s="<<sLim<<
",a="<<aLim
2182 if(!aLim) three=
false;
2193 else if(zLim<sLim&&r>zLim&&r<=sLim)
2203 else if(nLim<zLim&&r>nLim&&r<=zLim)
2225 else if(PDG==lPDG&&(lnCond&&lpCond&&llCond&&laCond))
2228 G4cout<<
"G4QN::EB:*l*: n="<<lnCond<<
",p="<<lpCond<<
",l="<<llCond<<
",a="<<laCond
2234 if(N&&GSResNL!=GSMass&&fMass+mNeut+GSResNL<totMass)
2235 nLim+=(N+N)*pow(totMass-mNeut-mLamb-GSResNL,2);
2237 if(Z&&GSResPL!=GSMass&&fMass+mProt+PBarr+GSResPL<totMass)
2239 if(barf) zLim+=(Z+Z)*pow(totMass-mProt-mLamb-PBarr-GSResPL,2);
2240 else zLim+=(Z+Z)*pow(totMass-mProt-mLamb-GSResPL,2);
2243 if(S>1&&GSResLL!=GSMass&&fMass+mLamb+GSResLL<totMass)
2244 sLim+=(S-1)*pow(totMass-mLamb-mLamb-GSResLL,2);
2246 if(evalph&&Z>1&&N>1&&a>4&&GSResLA!=GSMass&&fMass+mAlph+ABarr+GSResLA<totMass)
2248 if(barf) aLim+=pow(totMass-mLamb-mAlph-ABarr-GSResLA,2)*
2249 evalph*Z*(Z-1)*N*(N-1)*12/(a-2)/(a-3)/(a-4);
2250 else aLim+=pow(totMass-mLamb-mAlph-GSResLA,2)*
2251 evalph*Z*(Z-1)*N*(N-1)*12/(a-2)/(a-3)/(a-4);
2255 G4cout<<
"G4QN::EB:l, r="<<r<<
",n="<<nLim<<
",z="<<zLim<<
",s="<<sLim<<
",a="<<aLim
2259 if(!aLim) three=
false;
2270 else if(zLim<sLim&&r>zLim&&r<=sLim)
2280 else if(nLim<zLim&&r>nLim&&r<=zLim)
2302 else if(PDG==aPDG&&(anCond&&apCond&&alCond&&aaCond))
2305 G4cout<<
"G4QN::EB:*a*: n="<<anCond<<
",p="<<apCond<<
",l="<<alCond<<
",a="<<aaCond
2311 if(N>2&&GSResNA!=GSMass&&fMass+mNeut+ABarr+GSResNA<totMass)
2313 if(barf) nLim+=(N-2)*pow(totMass-mNeut-mAlph-ABarr-GSResNA,2);
2314 else nLim+=(N-2)*pow(totMass-mNeut-mAlph-GSResNA,2);
2317 if(Z>2&&GSResPA!=GSMass&&fMass+mProt+SAPBarr+GSResPA<totMass)
2319 if(barf) zLim+=(Z-2)*pow(totMass-mProt-mAlph-SAPBarr-GSResPA,2);
2320 else zLim+=(Z-2)*pow(totMass-mProt-mAlph-GSResPA,2);
2323 if(S&&GSResLA!=GSMass&&fMass+mLamb+ABarr+GSResLA<totMass)
2325 if(barf) sLim+=S*pow(totMass-mLamb-mAlph-ABarr-GSResLA,2);
2326 else sLim+=S*pow(totMass-mLamb-mAlph-GSResLA,2);
2329 if(evalph&&Z>3&&N>3&&a>7&&GSResAA!=GSMass&&fMass+mAlph+SAABarr+GSResAA<totMass)
2331 if(barf) aLim+=pow(totMass-mAlph-mAlph-SAABarr-GSResAA,2)*
2332 evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*12/(a-5)/(a-6)/(a-7);
2333 else aLim+=pow(totMass-mAlph-mAlph-GSResAA,2)*
2334 evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*12/(a-5)/(a-6)/(a-7);
2338 G4cout<<
"G4QN::EB:a, r="<<r<<
",n="<<nLim<<
",z="<<zLim<<
",s="<<sLim<<
",a="<<aLim
2342 if(!aLim) three=
false;
2353 else if(zLim<sLim&&r>zLim&&r<=sLim)
2363 else if(nLim<zLim&&r>nLim&&r<=zLim)
2388 if (rQPDG==pQPDG)rMass=mProt;
2389 else if(rQPDG==nQPDG)rMass=mNeut;
2390 else if(rQPDG==lQPDG)rMass=mLamb;
2393 G4cout<<
"G4QNucleus::EvaporateBary:evaBar="<<eMass<<bQPDG<<
",resN="<<rMass<<rQPDG
2394 <<
",secB="<<fMass<<
",three="<<three<<
G4endl;
2401 if(nFlag&&mNeut+GSResNn<totMass)
2403 G4double ken=totMass-mNeut-GSResNn;
2407 if(pFlag&&mProt+PBarr+GSResNp<totMass)
2409 G4double ken=totMass-mProt-GSResNp;
2410 if(barf) ken-=PBarr;
2414 if(lFlag&&mLamb+GSResNl<totMass)
2416 G4double ken=totMass-mLamb-GSResNl;
2420 if(evalph&&aFlag&&mAlph+GSResNa+ABarr<totMass)
2422 G4double ken=totMass-mAlph-GSResNa;
2423 if(barf) ken-=ABarr;
2424 aLim+=
CoulBarPenProb(ABarr,ken,2,4)*sqrt(ken)*evalph*Z*(Z-1)*N*(N-1)
2429 G4cout<<
"G4QNucl::EvapBar:2Decay r="<<r<<
",nLim="<<nLim<<
",zLim="<<zLim<<
",sLim="
2430 <<sLim<<
",nF="<<nFlag<<
",pF="<<pFlag<<
",lF="<<lFlag<<
",aF="<<aFlag<<
G4endl;
2439 else if(lFlag&&r>=zLim&&r<=sLim&&zLim<sLim)
2446 else if(nFlag&&r>=nLim&&r<=zLim&&nLim<zLim)
2453 else if(pFlag&&r<nLim)
2463 G4cout<<
"G4QNucleus::EvaporateBaryon: Photon #2-B#, dM="<<totMass-GSMass<<
G4endl;
2471 G4cout<<
"G4QNucl::EvaporateBaryon: b="<<eMass<<bQPDG<<
",r="<<rMass<<rQPDG<<
G4endl;
2477 G4cout<<
"G4QNucl::EvaporateBaryon:Decay in 3 particles"<<
G4endl;
2485 G4cout<<
"*G4QNucl::EvaporateBaryon:Decay M="<<totMass<<
",b="<<eMass<<bQPDG
2486 <<
",f="<<fMass<<fQPDG<<
",r="<<rMass<<rQPDG<<
G4endl;
2495 if(eMass+rMass<totMass&&cntr<cntm)
2498 G4cout<<
"G4QN::EvaB:eM="<<eMass<<
"+rM="<<rMass<<
" ="<<eMass+rMass<<
" < "<<totMass
2499 <<
",c="<<cntr<<
" < cm="<<cntm<<
", bPDG="<<bQPDG<<
", rPDG="<<rQPDG<<
G4endl;
2503 if (rQPDG==pQPDG)rMass=mProt;
2504 else if(rQPDG==nQPDG)rMass=mNeut;
2505 else if(rQPDG==lQPDG)rMass=mLamb;
2510 else if(totMass>mNeut+GSResNn)
2513 G4cout<<
"G4QNucl::EvaporateBaryon: Neutron , dM="<<totMass-GSResNn-mNeut<<
G4endl;
2520 else if(totMass>mProt+PBarr+GSResNp)
2523 G4cout<<
"G4QNucl::EvaporateBaryon: Proton , dM="<<totMass-GSResNp-mProt<<
G4endl;
2530 else if(totMass>mAlph+ABarr+GSResNa)
2533 G4cout<<
"G4QNucl::EvaporateBaryon: Alpha , dM="<<totMass-GSResNa-mAlph<<
G4endl;
2540 else if(totMass>GSMass)
2543 G4cout<<
"G4QNucl::EvaporateBaryon:Photon ### 2 ###, dM="<<totMass-GSMass<<
G4endl;
2552 G4cerr<<
"***G4QNucl::EvaporateBaryon: Cann't evaporate even gamma (1)"<<
G4endl;
2563 G4cout<<
"G4QNucleus::EvaporateBaryon: Decay in 2 baryons"<<
G4endl;
2568 if(aSecF)alp=evalph*Z*(Z-1)*N*(N-1)*10/(a-2)/(a-3)/(a-4);
2570 if(nnFlag&&totMass>mNeut+mNeut+GSResNN)
2571 nnLim+=N*(N-1)*pow(totMass-mNeut-mNeut-GSResNN,2);
2573 if(npFlag&&totMass>mNeut+mProt+PBarr+GSResNP)
2575 if(barf) nzLim+=2*N*Z*pow(totMass-mNeut-mProt-PBarr-GSResNP,2);
2576 else nzLim+=2*N*Z*pow(totMass-mNeut-mProt-GSResNP,2);
2579 if(ppFlag&&totMass>mProt+mProt+SPPBarr+GSResPP)
2581 if(barf) zzLim+=Z*(Z-1)*pow(totMass-mProt-mProt-SPPBarr-GSResPP,2);
2582 else zzLim+=Z*(Z-1)*pow(totMass-mProt-mProt-GSResPP,2);
2585 if(nlFlag&&totMass>mNeut+mLamb+GSResNL)
2586 nlLim+=2*N*S*pow(totMass-mNeut-mLamb-GSResNL,2);
2588 if(plFlag&&totMass>mProt+PBarr+mLamb+GSResPL)
2590 if(barf) zlLim+=2*Z*S*pow(totMass-mProt-mLamb-PBarr-GSResPL,2);
2591 else zlLim+=2*Z*S*pow(totMass-mProt-mLamb-GSResPL,2);
2594 if(llFlag&&totMass>mLamb+mLamb+GSResLL)
2595 llLim+=S*(S-1)*pow(totMass-mLamb-mLamb-GSResLL,2);
2597 if(naFlag&&totMass>mNeut+mAlph+ABarr+GSResNA)
2599 if(barf) naLim+=(N-3)*alp*pow(totMass-mNeut-mAlph-ABarr-GSResNA,2);
2600 else naLim+=(N-3)*alp*pow(totMass-mNeut-mAlph-GSResNA,2);
2603 if(paFlag&&totMass>mProt+SAPBarr+mAlph+GSResPA)
2605 if(barf) zaLim+=(Z-3)*alp*pow(totMass-mProt-mAlph-SAPBarr-GSResPA,2);
2606 else zaLim+=(Z-3)*alp*pow(totMass-mProt-mAlph-GSResPA,2);
2609 if(laFlag&&totMass>mLamb+mAlph+ABarr+GSResLA)
2611 if(barf) laLim+=S*alp*pow(totMass-mLamb-mAlph-ABarr-GSResLA,2);
2612 else laLim+=S*alp*pow(totMass-mLamb-mAlph-GSResLA,2);
2615 if(evalph&&aaFlag&&totMass>mAlph+mAlph+SAABarr+GSResAA)
2617 if(barf) aaLim+=alp*pow(totMass-mAlph-mAlph-SAABarr-GSResAA,2)*
2618 evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*7/(a-5)/(a-6)/(a-7);
2619 else aaLim+=alp*pow(totMass-mAlph-mAlph-GSResAA,2)*
2620 evalph*(Z-2)*(Z-3)*(N-2)*(N-3)*7/(a-5)/(a-6)/(a-7);
2631 G4cout<<
"G4QNuc::EvapBaryon: A+A, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2634 else if(aaLim&&zaLim<r&&r<laLim)
2642 G4cout<<
"G4QNuc::EvapBaryon: A+L, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2645 else if(aaLim&&naLim<r&&r<zaLim)
2653 G4cout<<
"G4QNuc::EvapBaryon: A+P, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2656 else if(aaLim&&llLim<r&&r<naLim)
2664 G4cout<<
"G4QNuc::EvapBaryon: A+N, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2667 else if(aaLim&&zlLim<r&&r<llLim)
2675 G4cout<<
"G4QNuc::EvapBaryon: L+L, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2678 else if(aaLim&&nlLim<r&&r<zlLim)
2686 G4cout<<
"G4QNuc::EvapBaryon: L+p, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2689 else if(aaLim&&zzLim<r&&r<nlLim)
2697 G4cout<<
"G4QNuc::EvapBaryon: L+n, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2701 else if(aaLim&&nzLim<r&&r<zzLim)
2709 G4cout<<
"G4QNuc::EvapBaryon: p+p, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2712 else if(aaLim&&nnLim<r&&r<nzLim)
2720 G4cout<<
"G4QNuc::EvapBaryon: n+p, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2723 else if(aaLim&&r<nnLim)
2731 G4cout<<
"G4QNuc::EvapBaryon: n+n, e="<<eMass<<
",f="<<fMass<<
",r="<<rMass<<
G4endl;
2738 G4cout<<
"G4QNucleus::EvaporateBaryon:Photon ### Decay in neutron ###"<<
G4endl;
2749 G4cout<<
"G4QNucleus::EvaporateBaryon:Photon ### Decay in proton ###"<<
G4endl;
2760 G4cout<<
"G4QNucleus::EvaporateBaryon:Photon ### Decay in alpha ###"<<
G4endl;
2771 G4cout<<
"G4QNucleus::EvaporateBaryon:Photon ### Decay in Lambda ###"<<
G4endl;
2782 G4cout<<
"G4QNuc::EvaporBaryon: Photon ### 3-Big ###,dM="<<totMass-GSMass<<
G4endl;
2798 G4cout<<
"*G4QNucl::EvaporateBaryon:Decay M="<<totMass<<
",b="<<eMass<<bQPDG
2799 <<
",f="<<fMass<<fQPDG<<
",r="<<rMass<<rQPDG<<
G4endl;
2808 G4cout<<
"G4QNuc::EvapBar:s="<<sma<<
",e="<<eMass<<
",f="<<fMass<<
",d="<<dma<<
",rM="
2816 if (rQPDG==pQPDG)rMass=mProt;
2817 else if(rQPDG==nQPDG)rMass=mNeut;
2818 else if(rQPDG==lQPDG)rMass=mLamb;
2825 G4cout<<
"***G4QNucleus::EvaporateBaryon: Emergency Decay M="<<totMass<<
",b="
2835 G4cout<<
"G4QNuc::EvapBaryon: Emergency decay is done for b="<<bQPDG<<h1->
GetQC()
2836 <<h1mom<<h1mom.
m()<<
", r="<<rQPDG<<h2->
GetQC()<<h2mom<<h2mom.
m()<<
G4endl;
2844 G4cout<<
"G4QNucleus::EvaporateBaryon: Decay in Baryon+Resid"<<
G4endl;
2848 if(nFlag&&mNeut+GSResNn<totMass)
2850 G4double ken=totMass-mNeut-GSResNn;
2854 if(pFlag&&mProt+PBarr+GSResNp<totMass)
2856 G4double ken=totMass-mProt-GSResNp;
2857 if(barf) ken-=PBarr;
2861 if(lFlag&&mLamb+GSResNl<totMass)
2863 G4double ken=totMass-mLamb-GSResNl;
2867 if(aFlag&&mAlph+GSResNa+ABarr<totMass)
2869 G4double ken=totMass-mAlph-GSResNa;
2870 if(barf) ken-=ABarr;
2871 aLim+=
CoulBarPenProb(ABarr,ken,2,4)*sqrt(ken)*evalph*Z*(Z-1)*N*(N-1)
2874 G4cout<<
"G4QNucleus::EvaporateBaryon:al="<<evalph<<
",k="<<ken<<
",P="
2880 G4cout<<
"G4QN::EB:DecIn2#2#r="<<r<<
",nL="<<nLim<<
",zL="<<zLim<<
",sL="<<sLim<<
",aL="
2881 <<aLim<<
",nF="<<nFlag<<
",pF="<<pFlag<<
",lF="<<lFlag<<
",aF="<<aFlag<<
G4endl;
2890 G4cout<<
"G4QNucleus::EvaporateBaryon: Decay in A + rA="<<GSResNa+mAlph<<
G4endl;
2893 else if(lFlag&&r>zLim&&r<sLim)
2900 G4cout<<
"G4QNucleus::EvaporateBaryon: Decay in L + rA="<<GSResNl+mLamb<<
G4endl;
2903 else if(pFlag&&r>nLim&&r<zLim)
2910 G4cout<<
"G4QNucleus::EvaporateBaryon: Decay in P + rA="<<GSResNp+mProt<<
G4endl;
2913 else if(nFlag&&r<nLim)
2920 G4cout<<
"G4QNucleus::EvaporateBaryon: Decay in N + rA="<<GSResNn+mNeut<<
G4endl;
2923 else if(mProt+GSResNp<totMass)
2926 G4cout<<
"G4QNucl::EvapBar: Emergency Proton, dM="<<totMass-GSResNp-mProt<<
G4endl;
2933 else if(mAlph+GSResNa<totMass)
2936 G4cout<<
"G4QNucl::EvapBar: Emergency Alpha, dM="<<totMass-GSResNa-mAlph<<
G4endl;
2946 G4cout<<
"G4QNuc::EvapBaryon: Photon ### 4-Big ###, dM="<<totMass-GSMass<<
G4endl;
2955 if (rQPDG==pQPDG)rMass=mProt;
2956 else if(rQPDG==nQPDG)rMass=mNeut;
2957 else if(rQPDG==lQPDG)rMass=mLamb;
2966 G4cout<<
"*G4QNucleus::EvaporateBaryon: Decay M="<<totMass<<
",b="<<bQPDG<<h1->
GetQC()
2967 <<eMass<<
",r="<<rQPDG<<h2->
GetQC()<<rMass<<
G4endl;
2972 G4cout<<
"G4QN::EvaB: **RESULT** b="<<bQPDG<<h1mom<<
", r="<<rQPDG<<h2mom<<
G4endl;
2979 G4cout<<
"G4QNucleus::EvaporateBaryon: Evaporation is done for b="<<bQPDG<<h1->
GetQC()
2980 <<h1mom<<h1mom.
m()<<
", r="<<rQPDG<<h2->
GetQC()<<h2mom<<h2mom.
m()<<
G4endl;
2987 G4cerr<<
"***G4QNucleus::EvaporateBaryon: ??? A=1"<<
G4endl;
3009 if(a<-2)
G4cerr<<
"***G4QNucleus::EvaporateBaryon: A="<<a<<
G4endl;
3022 if(r<s_value)
return 0;
3026 while(r>s_value && i<aN)
3040 static const G4int nOfMesons =45;
3041 static const G4int nOfBaryons=72;
3043 static G4int mesonPDG[nOfMesons] = {221,111,211,-211,311,321,-311,-321,331,
3045 223,113,213,-213,313,323,-313,-323,333,
3047 225,115,215,-215,315,325,-315,-325,335,
3049 227,117,217,-217,317,327,-317,-327,337,
3051 229,119,219,-219,319,329,-319,-329,339};
3053 static G4int baryonPDG[nOfBaryons]={2112,-2112,2212,-2212,3122,-3122,3112,-3112,
3055 3212,-3212,3222,-3222,3312,-3312,3322,-3322,
3057 1114,-1114,2114,-2114,2214,-2214,2224,-2224,3114,-3114,
3059 3214,-3214,3224,-3224,3314,-3314,3324,-3324,3334,-3334,
3061 2116,-2116,2216,-2216,3126,-3126,3116,-3116,
3063 3216,-3216,3226,-3226,3316,-3316,3326,-3326,
3065 1118,-1118,2118,-2118,2218,-2218,2228,-2228,3118,-3118,
3067 3218,-3218,3228,-3228,3318,-3318,3328,-3328,3338,-3338};
3072 G4int iQC = theQCandidates.size();
3073 if(iQC)
for(
G4int jq=0; jq<iQC; jq++)
delete theQCandidates[jq];
3074 theQCandidates.clear();
3075 if(maxMes>nOfMesons) maxMes=nOfMesons;
3076 if(maxMes>=0)
for (i=0; i<maxMes; i++)
3078 theQCandidates.push_back(
new G4QCandidate(mesonPDG[i]));
3080 G4cout<<
"G4QNucleus::InitCandidateVector: "<<ind++<<
", Meson # "<<i<<
" with code = "
3081 <<mesonPDG[i]<<
", QC="<<theQCandidates[i]->GetQC()<<
" is initialized"<<
G4endl;
3084 if(maxBar>nOfBaryons) maxBar=nOfBaryons;
3085 if(maxBar>=0)
for (i=0; i<maxBar; i++)
3088 G4cout<<
"G4QNucleus::InitCandidateVector: define PDG="<<baryonPDG[i]<<
G4endl;
3092 G4cout<<
"G4QNucleus::InitCandidateVector: current baryon is defined"<<
G4endl;
3094 theQCandidates.push_back(curBar);
3096 G4cout<<
"G4Nucleus::InitCandidateVector: "<<ind++<<
", Baryon # "<<i<<
" with code = "
3097 <<baryonPDG[i]<<
", QC="<<theQCandidates[i]->GetQC()<<
" is initialized"<<
G4endl;
3100 if(maxClst>=0)
for (i=0; i<maxClst; i++)
3106 theQCandidates.push_back(
new G4QCandidate(clusterPDG));
3108 G4cout<<
"G4QNucleus::InitCandidateVector:"<<ind++<<
", Cluster # "<<i<<
" with code = "
3136 G4int mCand=theQCandidates.size();
3139 if(comb<=0.) comb=1.;
3143 G4cout<<
"G4QN::PC:C#"<<mCand<<
",dZ="<<dZ<<
",dN="<<dN<<
",ZNS="<<Z<<
","<<N<<
","<<S<<
G4endl;
3198 if(piF&&gaF&&cBN!=1)
3205 if(cPDG==90001001)
G4cout<<
"G4QNuc::PrepCand: piF/gaF fragments are blocked"<<
G4endl;
3209 else if(cPDG<80000000&&(abs(cPDG)%10>4||cST>2))
3218 if(cPDG>80000000&&cPDG!=90000000)
3223 G4int nc = ac-zc-sc;
3227 if(ac<=maxClust&&pos>0.&&(pLV==zeroLV||intLV.
m()>.00001+cM))
3231 G4cout<<
"G4QNucleus::PrepareCand: ac="<<ac<<
", pV="<<pos<<
G4endl;
3234 if(ac && (piF || gaF))
3236 if (piF&&!gaF&&zc+ac) pos*=(zc+ac)/ac;
3237 else if(gaF&&!piF&&zc+dac) pos*=(zc+dac)/ac;
3240 if (ac==1&&pos>0.)dense=probVect[254]/pos;
3241 else if(ac==2&&pos>0.)dense=probVect[255]/pos;
3243 G4cout<<
"G4QNucleus::PrepC: cPDG="<<cPDG<<
",norm="<<pos<<
",zc="<<zc<<
",nc="<<nc
3244 <<
",sc="<<sc<<
",ac="<<ac<<
",ze1="<<ze1<<
",ne1="<<ne1<<
",se1="<<se1<<
G4endl;
3250 else if(nc) pos*=ne/ae;
3251 else if(sc) pos*=se/ae;
3258 G4cout<<
"G4QN::PrC:mp="<<mp<<
",pos="<<pos<<
",ae="<<ae
3259 <<
",Z="<<zc<<
",N="<<nc<<
",mac="<<mac<<
G4endl;
3266 if(ze<zc||ne<nc||se<sc) pos=0.;
3269 if (zc==2) pos*=ze*(ze-1)/aea;
3270 else if(nc==2) pos*=ne*(ne-1)/aea;
3271 else if(sc==2) pos*=se*(se-1)/aea;
3272 else if(zc==1&&nc==1) pos*=ze*ne/ae2;
3273 else if(zc==1&&sc==1) pos*=ze*se/ae2;
3274 else if(sc==1&&nc==1) pos*=se*ne/ae2;
3281 else G4cout<<
"***G4QNucl::PrepCand: z="<<zc<<
", n="<<nc<<
", s="<<sc<<
G4endl;
3287 G4cout<<
"G4QN::PrC:mp="<<mp<<
",p="<<pos<<
",A=2,(Z="<<zc<<
",N="<<nc<<
"),m="
3297 if(ac<ae1 && ac>0) comb*=(ae1-ac)/ac;
3304 G4cout<<
"G4QNuc::PrepCl:c="<<comb<<
",ac="<<ac<<
"("<<
index<<
"),m="<<mac<<
",a="
3309 if(ze1<=zc||ne1<=nc||se1<=sc) prod=0.;
3312 if(zc>0)
for(
int iz=1;
iz<=zc;
iz++) prod*=(ze1-
iz)/
iz;
3313 if(nc>0)
for(
int in=1;
in<=nc;
in++) prod*=(ne1-
in)/
in;
3314 if(sc>0)
for(
int is=1; is<=sc; is++) prod*=(se1-is)/is;
3320 if(pos)
G4cout<<
"G4QN::PreC:c="<<cPDG<<
",p="<<pos<<
",i="<<
index<<
",m="<<mac
3321 <<
",pr="<<prod<<
",c="<<cca<<
G4endl;
3330 G4cout<<
"G4QN::PrepC: ClusterPDG="<<cPDG<<
",preProb="<<pos<<
",d="<<dense<<
G4endl;
3343 G4cout<<
"G4QNucl::PrepCand:cPDG="<<cPDG<<
",pos="<<pos<<
G4endl;
3352 G4cout<<
"G4QNucl::PrepCand:covP="<<ae*sZ/ze<<
",covN="<<ae*sN/ne<<
",totP="<<sZ+sN<<
G4endl;
3365 if(rA < 0.)
G4cout<<
"-Warning-G4QNucl::CoulombBarGen: NucleusA="<<rA<<
", Z="<<rZ<<
G4endl;
3373 G4double cb=1.29*zz/(pow(ra,third)+pow(ca,third)+.1);
3380 G4cout<<
"G4QNucl::CoulBG:rA="<<cA<<
",rZ="<<cZ<<
",cA="<<cA<<
",cZ="<<cZ<<
",C="<<cb<<
G4endl;
3390 if (dA != 0.) rA-=dA;
3392 if(rA<0.)
G4cout<<
"-Warning-G4QNucl::CoulombBarrier: NucleusA="<<rA<<
", rZ="<<cZ<<
G4endl;
3395 if(delZ != 0.) rZ-=delZ;
3397 if(rA<0.)
G4cout<<
"-Warning-G4QNucl::CoulombBarrier: NucleusA="<<rA<<
", rZ="<<cZ<<
G4endl;
3407 if(cZ<=0.)
return 0.;
3413 G4double r=(pow(rA,third)+pow(cA,third))*(1.51+.00921*zz)/(1.+.009443*
zz);
3423 if(!cZ && !cA)
return Z*mProt+N*mNeut-
GetGSMass();
3425 G4int iZ=
static_cast<int>(cZ);
3426 G4int cN=
static_cast<int>(cA-cZ);
3431 G4int dz=
static_cast<int>(delZ);
3432 G4int dn=
static_cast<int>(dA-delZ);
3445 static const G4double dNeut= mNeut+mNeut;
3447 static const G4double dProt= mProt+mProt;
3452 static const G4double wellDebth=40.;
3459 if(B<1 || B>2)
return 1.;
3463 G4cout<<
"-Warning-G4QN::CBPP:SubtractedChrg="<<C<<
" >SubtractedBaryonNmbr="<<B<<
G4endl;
3487 G4cout<<
"G4QNucl::CBPenProb:GSM="<<GSM<<
",Z="<<Z<<
",N="<<N<<
",C="<<C<<
",B="<<B<<
G4endl;
3499 else if(B==1&&C==1) wD=
G4QNucleus(Z-1,N,S).GetGSMass()+mProt-GSM;
3500 else if(B==1&&C==0) wD=
G4QNucleus(Z,N-1,S).GetGSMass()+mNeut-GSM;
3507 if (!C) wD=
G4QNucleus(Z,N-2,S).GetGSMass()+dNeut-GSM;
3508 else if(C==1) wD=
G4QNucleus(Z-1,N-1,S).GetGSMass()+mDeut-GSM;
3509 else if(C==2) wD=
G4QNucleus(Z-2,N,S).GetGSMass()+dProt-GSM;
3529 G4cout<<
"G4QNucl::CBPenProb: wD="<<wD<<
",E="<<E<<
",CB="<<CB<<
G4endl;
3537 if(CBD<0.)
return 1.;
3538 if(ED<=0.)
return 0.;
3544 G4cout<<
"G4QN::CBPP:s="<<sR<<
",E="<<E<<
",w="<<wD<<
",CB="<<CB<<
",B="<<B<<
",C="<<C<<
G4endl;
3546 if(sR>=1.)
return 0.;
3562 while(x*x+y*y > 1.);
3563 std::pair<G4double, G4double> theImpactParameter;
3564 theImpactParameter.first = x*maxImpact;
3565 theImpactParameter.second = y*maxImpact;
3566 return theImpactParameter;
3573 G4cout<<
"G4QNucleus::ChooseNucleons: Nucleons search is started"<<rho0<<
G4endl;
3578 while (nucleons < theA)
3585 theNucleons.push_back(proton);
3587 else if ( (nucleons-protons) < N )
3591 theNucleons.push_back(neutron);
3593 else G4cout<<
"G4QNucleus::ChooseNucleons not efficient"<<
G4endl;
3602 static const G4double mProt2= mProt*mProt;
3603 static const G4double third= 1./3.;
3615 G4double nucDist2= nucleonDistance*nucleonDistance;
3622 while( i < theA && maxR > 0.)
3629 G4cout<<
"G4QNucl::ChoosePositions: i="<<i<<
", pos="<<aPos<<
", dens="<<density<<
G4endl;
3635 if(i==lastN) aPos=-rPos*sumPos.
unit();
3638 for(
G4int j=0; j<i && freeplace; j++)
3640 delta = places[j] - aPos;
3641 freeplace= delta.
mag2()>nucDist2;
3647 G4int nucPDG= theNucleons[i]->GetPDGCode();
3649 G4cout<<
"G4QNucl::ChoosePositions: frpl="<<freeplace<<
", nucPDG="<<nucPDG<<
G4endl;
3651 if(freeplace && nucPDG == 2212)
3654 G4double eFermi= sqrt(pFermi*pFermi+mProt2)-mProt;
3662 if( freeplace || failCNT > maxCNT )
3664 if( failCNT > maxCNT ) aPos=minPos;
3666 G4cout<<
"G4QNuc::ChoosePos:->> fill N["<<i<<
"], R="<<aPos<<
", f="<<failCNT<<
G4endl;
3678 G4cout<<
"G4QNucl::ChoosePositions: Out of the positioning LOOP"<<
G4endl;
3682 G4cout<<
"G4QNucl::ChoosePositions: The reduced summ is made"<<
G4endl;
3684 for(i=0; i<theA; i++) theNucleons[i]->
SetPosition(places[i]);
3687 G4cout <<
"G4QNucleus::ChoosePositions: The positions are defined for A="<<theA<<
G4endl;
3701 G4cout<<
"G4QNucleus::InitDensity: rA=iA=A="<<iA<<
", A^1/3="<<At<<
", A^2/3="<<At2<<
G4endl;
3708 G4cout<<
"-Warning-G4QNucl::ChoosePositions:L,iA="<<iA<<
",Radius(?)="<<radius<<
G4endl;
3711 rho0 = pow(2*
pi*radius, -1.5);
3716 G4double r0=1.16*(1.-1.16/At2)*fermi;
3720 G4cout<<
"-Warning-G4QNucl::ChoosePositions:H,iA="<<iA<<
",Radius(?)="<<radius<<
G4endl;
3724 if(!(rd<=0.1) && !(rd>-0.1))
3726 G4cout<<
"-Warning-G4QNucl::ChoosePositions:H,NAN,iA="<<iA<<
", rd="<<rd<<
G4endl;
3729 rho0=0.75/(
pi*pow(radius,3)*(1.+rd*rd*
pi2));
3742 return -exp((rPos-radius)/WoodSaxonSurf)*dens*dens*rho0/WoodSaxonSurf;
3750 return (maxRelDens>0 && maxRelDens <= 1. ) ? sqrt(-radius*log(maxRelDens) ) :
DBL_MAX;
3752 return (maxRelDens>0 && maxRelDens <= 1. ) ? (radius + WoodSaxonSurf*
3753 log((1.-maxRelDens+exp(-radius/WoodSaxonSurf))/maxRelDens) ) :
DBL_MAX;
3769 return constofpmax * pow(density*
GetA(),third);
3776 static const G4double mProt2= mProt*mProt;
3778 static const G4double third= 1./3.;
3786 G4cout<<
"G4QNucleus::ChooseFermiMomentum is called for Z="<<Z<<
", N="<<N<<
G4endl;
3788 for(i=0; i<theA; i++)
3795 if( i == am1) dir=-sumMom.
unit();
3800 if(eMax>mProt) mom=sqrt(eMax*eMax - mProt2)*rn3*
dir;
3802 else G4cerr<<
"-Warning-G4QNucleus::ChooseFermM: FailToGetProtonMomentum,p=0"<<
G4endl;
3805 else mom=ferm*rn3*
dir;
3809 G4cout<<
"G4QNucleus::ChooseFermiMomentum: for i="<<i<<
", candidate mom="<<mom<<
G4endl;
3825 for(i=0; i< theA ; i++ )
3841 G4cout<<
"G4QNucleus::ChooseFermiMomentum: FINALLY for i="<<i<<
", 4mom="<<tempV<<
G4endl;
3844 theNucleons[i]->Set4Momentum(tempV);
3847 G4cout<<
"G4QNucleus::ChooseFermiMomentum: FINALLY sum4M="<<sum<<
G4endl;
3857 for(
G4int i=0; i<theA; i++) vect[i]-=sum;
3866 G4cout<<
"-Warning-G4QNucleus::ReduceSum: *Failed* A="<<theA<<
" < 3"<<
G4endl;
3877 for(
G4int i=0; i<am1; i++) dp[i]=sum.
dot(vect[i]);
3881 for(
G4int i=0; i<am1; i++) if(dp[i]>0 && dp[i]<sum2)
3903 for(
G4int i=0; i<theA; i++) vect[i]-=sum;
3913 G4cout<<
"G4QNucleus::Init3D: is called currentNucleon="<<currentNucleon<<
G4endl;
3915 for_each(theNucleons.begin(),theNucleons.end(),
DeleteQHadron());
3916 theNucleons.clear();
3920 G4cout<<
"G4QNucleus::Init3D: Nucleons are initialized, nN="<<theNucleons.size()<<
G4endl;
3924 G4cout<<
"G4QNucl::Init3D: DensityPars for A="<<theA<<
":R="<<radius <<
",r0="<<rho0<<
G4endl;
3928 G4cout<<
"G4QNucleus::Init3D: Nucleons are positioned in the coordinate space"<<
G4endl;
3932 if(n3M.
x() || n3M.
y() || n3M.
z())
3938 G4cout<<
"G4QNucleus::Init3D: Nucleons are positioned in the momentum space"<<
G4endl;
3950 G4int theA=theNucleons.size();
3951 if(theA)
for(
G4int i=0; i<theA; i++)
3953 G4double nucr2=theNucleons[i]->GetPosition().mag2();
3954 if(nucr2 > maxradius2) maxradius2=nucr2;
3956 return sqrt(maxradius2)+nucleonDistance;
3963 G4double factor=(1.-sqrt(1.-bet2))/bet2;
3964 G4int theA=theNucleons.size();
3965 if(theA)
for (
G4int i=0; i< theA; i++)
3968 pos -= factor*(theBeta*pos)*theBeta;
3969 theNucleons[i]->SetPosition(pos);
3976 G4int theA=theNucleons.size();
3977 if(theA)
for(
G4int i=0; i<theA; i++)
3984 G4int theA=theNucleons.size();
3985 if(theA) currentNucleon=0;
3986 else G4cout<<
"-Warning-G4QNucleus::StartLoop: LOOP starts for uninited nucleons"<<
G4endl;
4023 for(
G4int i=0; i<
nt; ++i) sum+=i*Tb[i];
4026 G4cout<<
"G4QNucleus::GetTbIntegral:TI="<<sum<<
", RI="<<4*
pi*rho0*pow(radius,3)/3<<
G4endl;
4038 G4int nb =
static_cast<int>(bf);
4041 if(eb>=nt)
return 0.;
4044 return (nT-(bf-nb)*(nT-eT))/sfermi;
4053 G4cout<<
"-Warning-G4QNucleus::GetThickness: for A="<<tA<<
", => return 0"<<
G4endl;
4056 else if(tA==1)
return 0.;
4120 if (hPDG==2212) nPDG=90001000;
4121 else if(hPDG==2112) nPDG=90000001;
4122 else if(hPDG==3122||hPDG==3212) nPDG=91000000;
4123 else if(hPDG== 211) nPDG=90000999;
4124 else if(hPDG==-211) nPDG=89999001;
4125 else if(hPDG== 311) nPDG=89000001;
4126 else if(hPDG== 321) nPDG=89001000;
4127 else if(hPDG==-311) nPDG=90999999;
4128 else if(hPDG==-321) nPDG=90999000;
4129 else if(hPDG==1114) nPDG=89999002;
4130 else if(hPDG==2224) nPDG=90001999;
4131 else if(hPDG==3112) nPDG=90999000;
4132 else if(hPDG==3222) nPDG=91000999;
4133 else if(hPDG==3312) nPDG=91999000;
4134 else if(hPDG==3322) nPDG=91999999;
4135 else if(hPDG==3334) nPDG=92998999;
4136 else if(hPDG==-2212) nPDG=8999000;
4137 else if(hPDG==-2112) nPDG=8999999;
4138 else if(hPDG==-3122||hPDG==3212) nPDG=89000000;
4139 else if(hPDG==-3112) nPDG=89000999;
4140 else if(hPDG==-3222) nPDG=88999001;
4141 else if(hPDG==-3312) nPDG=88001000;
4142 else if(hPDG==-3322) nPDG=88000001;
4143 else if(hPDG==-3334) nPDG=87001001;
4151 if (nPDG==90001000) hPDG=2212;
4152 else if(nPDG==90000001) hPDG=2112;
4153 else if(nPDG==91000000) hPDG=3122;
4154 else if(nPDG==90000999) hPDG= 211;
4155 else if(nPDG==89999001) hPDG=-211;
4156 else if(nPDG==89001000) hPDG= 213;
4157 else if(nPDG==89000001) hPDG= 213;
4158 else if(nPDG==90999000) hPDG=-213;
4159 else if(nPDG==90999999) hPDG=-213;
4160 else if(nPDG==90001999) hPDG=1114;
4161 else if(nPDG==89999002) hPDG=2224;
4162 else if(nPDG==91000999) hPDG=3112;
4163 else if(nPDG==90999001) hPDG=3222;
4164 else if(nPDG==91999999) hPDG=3312;
4165 else if(nPDG==91999000) hPDG=3322;
4166 else if(nPDG==92998999) hPDG=3334;
4191 G4cout<<
"G4QNucleus::EvaporateNucleus:-Called-PDG="<<thePDG<<
",QC="<<theQC<<
G4endl;
4195 G4cout<<
"G4QNucleus::EvaporateNucleus: theBN="<<theBN<<
G4endl;
4199 if(!theBN || thePDG<80000000 || thePDG==90000000)
4204 if(thePDG==90000000)
4207 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (1) qH=0"<<
G4endl;
4210 G4cout<<
"-Warning-G4QNucleus::EvapNuc:vacuum,4Mom="<<q4M<<
G4endl;
4214 G4cout<<
"-Warning-G4QNucleus::EvapNuc:vacuum, Called for Meson PDG="<<thePDG<<
G4endl;
4215 evaHV->push_back(qH);
4222 G4int SSS=(thePDG-90000000)/1000000;
4225 theQC = qH->
GetQC();
4227 G4cout<<
"=>Hyper Change=>G4QNucleus::EvaporateNuceus: NewNucPDG="<<thePDG<<
G4endl;
4236 evaHV->push_back(qH);
4241 G4int theN=theBN-theC-theS;
4245 G4cout<<
"G4QNucleus::EvaporateNucleus(EVA):===IN==> PDG="<<thePDG<<
",4Mom="<<q4M<<
", B="
4246 <<theBN<<
", Z="<<theC<<
", N="<<theN<<
", S="<<theS<<
G4endl;
4250 G4cout<<
"-Warning-G4QNuc::EvapNuc: Evapor of anti-nuclei is not implemented yet PDG="
4252 evaHV->push_back(qH);
4255 else if(thePDG==91000000||thePDG==90001000||thePDG==90000001)
4259 if(thePDG==90001000) gsM=mProt;
4260 else if(thePDG==91000000) gsM=mLamb;
4261 if(fabs(totMass-gsM)<.001)
4266 evaHV->push_back(qH);
4269 else if(totMass<gsM)
4272 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (2) qH=0"<<
G4endl;
4278 ed <<
"Baryon is below Mass Shell: Baryon " << thePDG
4279 <<
" is belowMassShell, M=" << totMass <<
G4endl;
4280 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0000",
4287 G4cout<<
"G4QN::EvaporNucl: PDG="<<thePDG<<
",M="<<totMass<<
">"<<gsM<<
",d="<<d<<
G4endl;
4293 if(thePDG==90001000)
4300 else if(thePDG==90000001)
4318 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (3) qH=0"<<
G4endl;
4324 ed <<
"BaryonDecay In Baryon+Gam Error: h=" << thePDG <<
"(GSM="
4325 << gsM <<
")+gam>tM=" << totMass <<
G4endl;
4326 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0001",
4330 G4cout<<
"G4QNucl::EvaNuc:"<<totMass<<q4M<<
"->"<<thePDG<<h4Mom<<
"+g="<<g4Mom<<
",n="
4335 G4cout<<
"G4QNucleus::EvaporateNucleus: Hadr="<<thePDG<<h4Mom<<
G4endl;
4337 evaHV->push_back(curH);
4340 G4cout<<
"G4QNucleus::EvaporateNucleus: Gamma(pion)4M="<<g4Mom<<
G4endl;
4342 evaHV->push_back(curG);
4344 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (4) qH=0"<<
G4endl;
4349 else if(thePDG==89000000||thePDG==89999000||thePDG==89999999)
4353 if(thePDG==89999000) gsM=mProt;
4354 else if(thePDG==89000000) gsM=mLamb;
4355 if(fabs(totMass-gsM)<.001)
4360 evaHV->push_back(qH);
4363 else if(totMass<gsM)
4366 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (2a) qH=0"<<
G4endl;
4372 ed <<
"anti-Baryon is below Mass Shell: antiBaryon=" << thePDG
4373 <<
" below MassShell, M=" << totMass <<
G4endl;
4374 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0002",
4381 G4cout<<
"G4QN::EvaporNucl: PDG="<<thePDG<<
",M="<<totMass<<
">"<<gsM<<
",d="<<d<<
G4endl;
4387 if(thePDG==89999000)
4394 else if(thePDG==89999999)
4412 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (3a) qH=0"<<
G4endl;
4418 ed <<
"BaryonDecay In Baryon+Gam Error: ah=" << thePDG <<
"(GSM="
4419 << gsM <<
")+gam>tM=" << totMass <<
G4endl;
4420 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0003",
4424 G4cout<<
"G4QNucl::EvaNuc:aM="<<totMass<<q4M<<
"->"<<thePDG<<h4Mom<<
"+g="<<g4Mom<<
",n="
4429 G4cout<<
"G4QNucleus::EvaporateNucleus: antiHadr="<<thePDG<<h4Mom<<
G4endl;
4431 evaHV->push_back(curH);
4434 G4cout<<
"G4QNucleus::EvaporateNucleus: (anti) Gamma(pion)4M="<<g4Mom<<
G4endl;
4436 evaHV->push_back(curMes);
4438 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (4a) qH=0"<<
G4endl;
4443 else if((thePDG==90001999||thePDG==89999002) && totMass>1080.)
4449 if(thePDG==90001999)
4455 if(fabs(totMass-gsM-mPi)<.001)
4462 evaHV->push_back(curB);
4465 evaHV->push_back(curM);
4467 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (5) qH=0"<<
G4endl;
4472 else if(totMass<gsM+mPi)
4475 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (6) qH=0"<<
G4endl;
4481 ed <<
"Delta is below Mass Shell: Delta " << thePDG
4482 <<
" is belowMassShell, M=" << totMass <<
G4endl;
4483 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0004",
4493 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (7) qH=0"<<
G4endl;
4499 ed <<
"DeltaDecInBaryon+Pi Error: Dh=" << thePDG <<
"N+pi=" << gsM+mPi
4500 <<
">tM=" << totMass <<
G4endl;
4501 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0005",
4505 G4cout<<
"G4QNuc::EvaNuc:"<<totMass<<q4M<<
"->"<<thePDG<<h4Mom<<
"+pi="<<g4Mom<<
", nH="
4510 G4cout<<
"G4QNucleus::EvaporateNucl: Nucleon="<<thePDG<<h4Mom<<
G4endl;
4512 evaHV->push_back(curH);
4517 evaHV->push_back(curG);
4519 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (8) qH=0"<<
G4endl;
4524 else if((thePDG==89998001||thePDG==90000998) && totMass>1080.)
4530 if(thePDG==89998001)
4536 if(fabs(totMass-gsM-mPi)<.001)
4543 evaHV->push_back(curB);
4546 evaHV->push_back(curM);
4548 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (5a) qH=0"<<
G4endl;
4553 else if(totMass<gsM+mPi)
4556 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (6a) qH=0"<<
G4endl;
4562 ed <<
"anti-Delta is below Mass Shell: aDelta " << thePDG
4563 <<
" is belowMassShell, M=" << totMass <<
G4endl;
4564 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0006",
4574 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (7a) qH=0"<<
G4endl;
4580 ed <<
"AntiDeltaDecayInBaryon+Pi Error: aD=" << thePDG <<
"N+pi="
4581 << gsM+mPi <<
">tM=" << totMass <<
G4endl;
4582 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0007",
4586 G4cout<<
"G4QNuc::EvaNuc:(aD) "<<totMass<<q4M<<
"->"<<thePDG<<h4Mom<<
" + pi="<<g4Mom
4587 <<
", nH="<<evaHV->size()<<
G4endl;
4591 G4cout<<
"G4QNucleus::EvaporateNucl: Nucleon="<<thePDG<<h4Mom<<
G4endl;
4593 evaHV->push_back(curH);
4598 evaHV->push_back(curMes);
4600 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (8a) qH=0"<<
G4endl;
4607 else if((thePDG==89999003 || thePDG==90002999) && totMass>2020.)
4610 G4int nucPDG = 2112;
4613 if(thePDG==90002999)
4619 if(totMass>mPi+nucM+nucM)
4627 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (9) qH=0"<<
G4endl;
4634 ed <<
" ISO-Dibaryon DecayIn3 error: tM=" << totMass <<
"-> 2N="
4635 << nucPDG <<
"(M=" << nucM <<
") + pi=" << piPDG <<
"(M="
4637 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0008",
4641 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (10) qH=0"<<
G4endl;
4646 G4cout<<
"G4QNucleus::EvaporateNucleus: Bar1="<<nucPDG<<n14M<<
G4endl;
4648 evaHV->push_back(h1H);
4651 G4cout<<
"G4QNucleus::EvaporateNucleus: Bar2="<<nucPDG<<n24M<<
G4endl;
4653 evaHV->push_back(h2H);
4656 G4cout<<
"G4QNucleus::EvaporateNucleus: Pi="<<piPDG<<pi4M<<
G4endl;
4658 evaHV->push_back(piH);
4663 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (11) qH=0"<<
G4endl;
4670 ed <<
"ISO-DiBaryon is under MassShell: IdPDG=" << thePDG <<
", q4M="
4671 << q4M <<
", M=" << totMass <<
" < M_2N+Pi, d=" << totMass-2*nucM-mPi
4673 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0009",
4677 else if((thePDG==90000002||thePDG==90001001||thePDG==90002000)&&totMass>2020.)
4685 if (thePDG==90002000) piPDG = 211;
4686 else if(thePDG==90000002) piPDG = -211;
4702 if(totMass>mPi+n1M+n2M)
4713 ed <<
"ISO-dibaryon excit. DecayIn3 error: IsoDN, tM=" << totMass
4714 <<
"-> N1=" << n1PDG <<
"(M=" << n1M <<
") + N2=" << n2PDG
4715 <<
"(M=" << n2M <<
") + pi=" << piPDG <<
" (Mpi=" << mPi <<
")"
4717 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0010",
4721 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (12) qH=0"<<
G4endl;
4726 G4cout<<
"G4QNucleus::EvaporateNucleus: Bar1="<<n1PDG<<n14M<<
G4endl;
4728 evaHV->push_back(h1H);
4731 G4cout<<
"G4QNucleus::EvaporateNucleus: Bar2="<<n2PDG<<n24M<<
G4endl;
4733 evaHV->push_back(h2H);
4736 G4cout<<
"G4QNucleus::EvaporateNucleus: Pi="<<piPDG<<pi4M<<
G4endl;
4738 evaHV->push_back(piH);
4743 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (13) qH=0"<<
G4endl;
4750 ed <<
"IsoDiBarState is under MassShell: IdPDG=" << thePDG <<
", q4M="
4751 << q4M <<
", M=" << totMass <<
" < M1+M2+Pi, d="
4752 << totMass-n1M-n2M-mPi <<
G4endl;
4753 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0011",
4758 else if((thePDG==90000997 || thePDG==89997001) && totMass>2020.)
4761 G4int nucPDG = -2112;
4764 if(thePDG==90002999)
4770 if(totMass>mPi+nucM+nucM)
4778 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (9a) qH=0"<<
G4endl;
4785 ed <<
"Anti-ISO-DibaryonDecayIn3 error: antiM=" << totMass
4786 <<
"-> 2aN=" << nucPDG <<
"(M=" << nucM <<
") + pi=" << piPDG
4787 <<
"(M=" << mPi <<
")" <<
G4endl;
4788 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0012",
4792 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (10a) qH=0"<<
G4endl;
4797 G4cout<<
"G4QNucleus::EvaporateNucleus:(I) antiBar1="<<nucPDG<<n14M<<
G4endl;
4799 evaHV->push_back(h1H);
4802 G4cout<<
"G4QNucleus::EvaporateNucleus:(I) antiBar2="<<nucPDG<<n24M<<
G4endl;
4804 evaHV->push_back(h2H);
4807 G4cout<<
"G4QNucleus::EvaporateNucleus:(I) (anti) Pi="<<piPDG<<pi4M<<
G4endl;
4809 evaHV->push_back(piH);
4814 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (11a) qH=0"<<
G4endl;
4821 ed <<
"AntiISODiBaryon is underMassShell: aIdPDG=" << thePDG <<
", q4M="
4822 << q4M <<
", M=" << totMass <<
" < M_2N+Pi, d=" << totMass-2*nucM-mPi
4824 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0013",
4828 else if((thePDG==89999998||thePDG==89998999||thePDG==89998000)&&totMass>2020.)
4831 G4int n1PDG = -2212;
4832 G4int n2PDG = -2112;
4836 if (thePDG==89998000) piPDG = -211;
4837 else if(thePDG==89999998) piPDG = 211;
4853 if(totMass>mPi+n1M+n2M)
4864 ed <<
"AntiExcitedDibaryon DecayIn3 error: IsoDN,antiM=" << totMass
4865 <<
"-> N1=" << n1PDG <<
"(M=" << n1M <<
") + N2=" << n2PDG <<
"(M="
4866 << n2M <<
") + pi=" << piPDG <<
" (Mpi=" << mPi <<
")" <<
G4endl;
4867 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0014",
4871 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (12a) qH=0"<<
G4endl;
4876 G4cout<<
"G4QNucleus::EvaporateNucleus: antiBar1="<<n1PDG<<n14M<<
G4endl;
4878 evaHV->push_back(h1H);
4881 G4cout<<
"G4QNucleus::EvaporateNucleus: antiBar2="<<n2PDG<<n24M<<
G4endl;
4883 evaHV->push_back(h2H);
4886 G4cout<<
"G4QNucleus::EvaporateNucleus: (anti)Pi="<<piPDG<<pi4M<<
G4endl;
4888 evaHV->push_back(piH);
4893 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (13a) qH=0"<<
G4endl;
4900 ed <<
"AntiDiBarState is under MassShell: andPDG=" << thePDG <<
", q4M="
4901 << q4M <<
", M=" << totMass <<
" < M1+M2+Pi, d="
4902 << totMass-n1M-n2M-mPi <<
G4endl;
4903 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0015",
4908 else if(fabs(totMass-totGSM)<.001)
4913 if(thePDG==90004004)
4917 else if(thePDG==90004002)
4921 else if((theC==theBN||theN==theBN||theS==theBN)&&theBN>1)
4931 evaHV->push_back(qH);
4934 else if(theBN>1 && thePDG>88000000 && thePDG<89000000)
4936 G4cout<<
"---Warning---G4QNucl::EvaNuc:MustNotBeHere.PDG="<<thePDG<<
",S="<<theS<<
G4endl;
4943 G4int nucPDG = thePDG;
4944 if(bZ>=bN) nucPDG+=999000;
4951 if(bZ>bN) nucPDG+=999000;
4959 G4cout<<
"-Warning-G4QN::EvN:M="<<nucM<<
","<<nucPDG<<
",1="<<k1PDG<<
",2="<<k2PDG<<
G4endl;
4966 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (14) qH=0"<<
G4endl;
4973 ed <<
" Nucleus+2antiK DecayIn3 error: tM=" << totMass <<
"-> N="
4974 << nucPDG <<
"(M=" << nucM <<
") + k1=" << k1PDG <<
"(M=" << mK1
4975 <<
") + k2=" << k2PDG <<
"(M=" << mK2 <<
")" <<
G4endl;
4976 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0016",
4980 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (15) qH=0"<<
G4endl;
4985 G4cout<<
"G4QNucleus::EvaporateNucleus: k1="<<k1PDG<<k14M<<
G4endl;
4987 evaHV->push_back(k1H);
4990 G4cout<<
"G4QNucleus::EvaporateNucleus: k2="<<k2PDG<<k24M<<
G4endl;
4992 evaHV->push_back(k2H);
4995 G4cout<<
"G4QNucleus::EvaporateNucleus: Nuc="<<nucPDG<<n4M<<
G4endl;
4997 evaHV->push_back(nH);
5000 else if ( (thePDG > 80000000 && thePDG != 90000000) ||
5001 thePDG == 2112 || thePDG == 2212 || thePDG == 3122)
5005 if (thePDG==2112) thePDG=90000001;
5006 else if(thePDG==2212) thePDG=90001000;
5007 else if(thePDG==3122) thePDG=91000000;
5016 G4cout<<
"G4QN::EvaNuc: theBN="<<theBN<<
", bA="<<bA<<
", bZ="<<bZ<<
", bN="<<bN<<
G4endl;
5020 if(bZ==2&&bN==5)
G4cout<<
"G4QNucleus::EvaNucl: (2,5) GSM="<<GSMass<<
" > "
5022 if(bZ==1&&bN==3)
G4cout<<
"G4QNucl::EvaNucl: (1,3) GSM="<<GSMass<<
" > "
5025 G4cout<<
"G4QNucl::EvaNuc:"<<qNuc<<
",PDG="<<thePDG<<
",M="<<totMass<<
",dM="<<dM<<
G4endl;
5031 G4cout<<
"G4QNucleus::EvaporateNuc: bs="<<bsCond<<
", dbs="<<dbsCond<<
", A="<<bA<<
G4endl;
5033 if(fabs(totMass-GSMass)<.003&&!bsCond&&!dbsCond)
5038 evaHV->push_back(qH);
5041 else if ( ( bA == 1 || (!bsCond && !dbsCond) ) && totMass > GSMass+.003 )
5045 G4cout<<
"G4QN::EvaN:SplitBar, s="<<bsCond<<
",M="<<totMass<<
" > GSM="<<GSMass<<
G4endl;
5047 G4int nOfOUT = evaHV->size();
5051 G4QHadron* theLast = (*evaHV)[nOfOUT-1];
5056 G4cout<<
"G4QN::EvaNuc:*BackFus*,BN="<<lastBN<<
",nF="<<nFrag<<
",n="<<nOfOUT<<
G4endl;
5060 G4QHadron* thePrev = (*evaHV)[nOfOUT-2];
5065 G4cout<<
"G4QNucl::EvaNucl: DelTheLast, nFr="<<nFrag<<
", pPDG="<<prevPDG<<
G4endl;
5075 if(lastBN<1&&nOfOUT>1)
5077 G4QHadron* thePrev = (*evaHV)[nOfOUT-2];
5080 evaHV->push_back(theLast);
5081 evaHV->push_back(thePrev);
5095 G4cout<<
"G4QN::EvaNuc: tM="<<totMass<<
"-LM="<<lastM<<lastQC<<
"-GSM="<<GSMass<<
"="
5106 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (16) qH=0"<<
G4endl;
5110 G4cout<<
"G4QNucleus::EvaporateNucl: EVH "<<totPDG<<q4M<<
" fill AsIs"<<
G4endl;
5113 else evaHV->push_back(evH);
5123 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (17) qH=0"<<
G4endl;
5127 G4cout<<
"***G4QNucleus::EvaNucl: EVH "<<totPDG<<q4M<<
" fill AsIs"<<
G4endl;
5129 evaHV->push_back(evH);
5131 G4cout<<
"***G4QN::EvaN:DecayIn L"<<lastQC<<
"+RN"<<totQC<<
" failed"<<
G4endl;
5138 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (18) qH=0"<<
G4endl;
5144 G4cout<<
"G4QNucleus::EvaNuc:fill NH "<<totPDG<<r4Mom<<
G4endl;
5147 if(thePDG==92000000||thePDG==90002000||thePDG==90000002)
5149 else evaHV->push_back(nucH);
5169 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (19) qH=0"<<
G4endl;
5175 ed <<
" Decay in Gamma failed: h=" << thePDG <<
"(GSM=" << GSMass
5176 <<
")+g>tM=" << totMass <<
G4endl;
5177 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0017",
5181 G4cout<<
"G4QNuc::EvaNuc: "<<q4M<<
"->totResN="<<thePDG<<h4Mom<<
"+g="<<g4Mom<<
G4endl;
5185 G4cout<<
"G4QNucleus::EvaporateNucleus: Fill a Fragment="<<thePDG<<h4Mom<<
G4endl;
5187 if(thePDG==92000000||thePDG==90002000||thePDG==90000002)
DecayDibaryon(curH,evaHV);
5188 else evaHV->push_back(curH);
5191 G4cout<<
"G4QNucleus::EvaporateNucleus: Fill a Gamma="<<g4Mom<<
G4endl;
5193 evaHV->push_back(curG);
5195 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (20) qH=0"<<
G4endl;
5203 else if(totMass<GSMass+.003&&(bsCond||dbsCond))
5206 G4cout<<
"G4QN::EvN:2B="<<dbsCond<<
",B="<<bsCond<<
",M="<<totMass<<
"<"<<GSMass<<
G4endl;
5210 if(bN==4&&bZ==2&&!bS)
5217 if(bsCond==2112&&bN>0&&bA>1)
5222 if (nResPDG==90000001) nResM=mNeut;
5223 else if(nResPDG==90001000) nResM=mProt;
5224 else if(nResPDG==91000000) nResM=mLamb;
5229 if(bsCond==2212&&bZ>0&&bA>1)
5234 if (pResPDG==90000001) pResM=mNeut;
5235 else if(pResPDG==90001000) pResM=mProt;
5236 else if(pResPDG==91000000) pResM=mLamb;
5241 if(bsCond==3122&&bS>0&&bA>1)
5246 if (lResPDG==90000001) lResM=mNeut;
5247 else if(lResPDG==90001000) lResM=mProt;
5248 else if(lResPDG==91000000) lResM=mLamb;
5253 if(bsCond==90001001&&bN>0&&bZ>0&&bA>2)
5258 if (dResPDG==90000001) dResM=mNeut;
5259 else if(dResPDG==90001000) dResM=mProt;
5260 else if(dResPDG==91000000) dResM=mLamb;
5265 if(bsCond==90002002&&bN>1&&bZ>1&&bA>4)
5270 if (aResPDG==90000001) aResM=mNeut;
5271 else if(aResPDG==90001000) aResM=mProt;
5272 else if(aResPDG==91000000) aResM=mLamb;
5277 if(dbsCond&&bN>1&&bA>2)
5282 if (nnResPDG==90000001) nnResM=mNeut;
5283 else if(nnResPDG==90001000) nnResM=mProt;
5284 else if(nnResPDG==91000000) nnResM=mLamb;
5289 if(dbsCond&&bZ>1&&bA>2)
5294 if (ppResPDG==90000001) ppResM=mNeut;
5295 else if(ppResPDG==90001000) ppResM=mProt;
5296 else if(ppResPDG==91000000) ppResM=mLamb;
5301 if(dbsCond&&bN>0&&bZ>0&&bA>2)
5306 if (npResPDG==90000001) npResM=mNeut;
5307 else if(npResPDG==90001000) npResM=mProt;
5308 else if(npResPDG==91000000) npResM=mLamb;
5313 if(dbsCond&&bN>0&&bS>0&&bA>2)
5318 if (lnResPDG==90000001) lnResM=mNeut;
5319 else if(lnResPDG==90001000) lnResM=mProt;
5320 else if(lnResPDG==91000000) lnResM=mLamb;
5325 if(dbsCond&&bS>0&&bZ>0&&bA>2)
5330 if (lpResPDG==90000001) lpResM=mNeut;
5331 else if(lpResPDG==90001000) lpResM=mProt;
5332 else if(lpResPDG==91000000) lpResM=mLamb;
5337 if(dbsCond&&bS>1&&bA>2)
5342 if (llResPDG==90000001) llResM=mNeut;
5343 else if(llResPDG==90001000) llResM=mProt;
5344 else if(llResPDG==91000000) llResM=mLamb;
5348 G4cout<<
"G4QNucleus::EvaNucl: rP="<<pResPDG<<
",rN="<<nResPDG<<
",rL="<<lResPDG<<
",N="
5349 <<bN<<
",Z="<<bZ<<
",nL="<<bS<<
",totM="<<totMass<<
",n="<<totMass-nResM-mNeut
5350 <<
",p="<<totMass-pResM-mProt<<
",l="<<totMass-lResM-mLamb<<
G4endl;
5352 if ( thePDG == 90004004 ||
5353 (thePDG == 90002004 && totMass > mHel6+.003) ||
5354 (bA > 4 && bsCond && bN > 1 && bZ > 1 && totMass > aResM+mAlph) ||
5355 (bA > 1 && bsCond && ( (bN > 0 && totMass > nResM+mNeut) ||
5356 (bZ > 0 && totMass > pResM+mProt) ||
5357 (bS > 0 && totMass > lResM+mLamb) ) ) ||
5359 (( bN > 0 && bZ > 0 &&
5360 ( (bsCond && totMass > dResM+mDeut) || (dbsCond && totMass > dResM+mDeut) )
5361 ) || ( dbsCond && ( (bN > 1 && totMass > nnResM+mNeut+mNeut) ||
5362 (bZ > 1 && totMass > ppResM+mProt+mProt) ||
5363 (bS > 1 && totMass > llResM+mLamb+mLamb) ||
5364 (bN && bS && totMass > lnResM+mLamb+mNeut) ||
5365 (bZ && bS && totMass > lpResM+mLamb+mProt)
5372 G4int barPDG = 90002002;
5373 G4int resPDG = 90002002;
5379 if(gResPDG&&tMC>mHel6+.003)
5386 else if(nResPDG&&tMC>nResM+mNeut)
5393 else if(pResPDG&&totMass+.001>pResM+mProt)
5400 else if(lResPDG&&tMC>lResM+mLamb)
5407 else if(thePDG!=90004004&&bN>1&&bZ>1&&bA>4&&tMC>aResM+mAlph)
5414 else if(dResPDG&&tMC>dResM+mDeut)
5421 else if(ppResPDG&&tMC>ppResM+mProt+mProt)
5430 else if(nnResPDG&&tMC>nnResM+mNeut+mNeut)
5438 else if(npResPDG&&tMC>npResM+mProt+mNeut)
5446 else if(lnResPDG&&tMC>lnResM+mLamb+mNeut)
5454 else if(lpResPDG&&tMC>lpResM+mLamb+mProt)
5463 else if(llResPDG&&tMC>llResM+mLamb+mLamb)
5472 else if(thePDG!=90004004&&tMC>GSMass)
5479 else if(thePDG!=90004004)
5484 ed <<
"M<GSM & can't decayInPNL: PDG=" << thePDG <<
",M=" << totMass
5485 <<
"< GSM=" << GSMass <<
G4endl;
5486 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0018",
5495 evaHV->push_back(qH);
5496 G4cout<<
"---Warning---G4QNucleus::EvaNuc:rP="<<pResPDG<<
",rN="<<nResPDG<<
",rL="
5497 <<lResPDG<<
",N="<<bN<<
",Z="<<bZ<<
",L="<<bS<<
",totM="<<totMass<<
",n="
5498 <<totMass-nResM-mNeut<<
",p="<<totMass-pResM-mProt<<
",l="
5499 <<totMass-lResM-mLamb<<
G4endl;
5500 G4cout<<
"---Warning---G4QN::EvN:DecIn2Error b="<<barPDG<<
",r="<<resPDG<<
G4endl;
5506 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (21) qH=0"<<
G4endl;
5511 G4cout<<
"G4QNucleus::EvaNucleus:(1) Baryon="<<barPDG<<a4Mom<<
G4endl;
5513 evaHV->push_back(HadrB);
5516 G4cout<<
"G4QNucleus::EvaNucleus:(1) Residual="<<resPDG<<b4Mom<<
G4endl;
5520 else evaHV->push_back(HadrR);
5526 if(!qH->
DecayIn3(a4Mom,b4Mom,c4Mom))
5528 evaHV->push_back(qH);
5529 G4cout<<
"---Warning---G4QN::EvN:rNN="<<nnResPDG<<
",rNP="<<npResPDG<<
",rPP="
5530 <<ppResPDG<<
",N="<<bN<<
",Z="<<bZ<<
",L="<<bS<<
",tM="<<totMass<<
",nn="
5531 <<totMass-nnResM-mNeut-mNeut<<
",np="<<totMass-npResM-mProt-mNeut<<
",pp="
5532 <<totMass-ppResM-mProt-mProt<<
G4endl;
5533 G4cout<<
"---Warning---G4QN::EvN:DecIn2Error,b="<<barPDG<<
",r="<<resPDG<<
G4endl;
5539 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (22) qH=0"<<
G4endl;
5544 G4cout<<
"G4QNucleus::EvaporateNucleus:(2) Baryon1="<<barPDG<<a4Mom<<
G4endl;
5546 evaHV->push_back(HadrB);
5549 G4cout<<
"G4QNucleus::EvaporateNucleus:(2) Baryon2="<<thdPDG<<c4Mom<<
G4endl;
5551 evaHV->push_back(HadrC);
5554 G4cout<<
"G4QNucleus::EvaporateNucleus:(2) Residual="<<resPDG<<b4Mom<<
G4endl;
5558 else evaHV->push_back(HadrR);
5562 else if (fabs(totMass-GSMass)<.003)
5565 G4cout<<
"*|*|*|*G4QNucleus::EvaporateNuc: fill AsIs. Should never be here"<<
G4endl;
5567 evaHV->push_back(qH);
5573 G4cout<<
"***G4QNucl::EvaNuc: tM="<<totMass<<
"("<<thePDG<<
") < GSM="<<GSMass<<
", d="
5576 evaHV->push_back(qH);
5583 G4cout<<
"G4QN::EvaNuc:***EVA***tPDG="<<thePDG<<
",M="<<totMass<<
">GSM="<<GSMass<<
",d="
5601 while(evC&&evcn<evcm)
5608 G4cout<<
"***G4QNuc::EvaNuc:***EVA Failed***PDG="<<thePDG<<
",M="<<totMass<<
G4endl;
5615 evaHV->push_back(qH);
5632 G4cout<<
"G4QNucl::EvaNuc:Attempt #"<<evcn<<
" > "<<evcm<<
", rPDG="<<rPDG<<
", bPDG="
5633 <<bPDG<<
", bE="<<b4M.
e()-b4M.
m()<<
" > bCB="<<bCB<<
G4endl;
5638 G4cout<<
"G4QNucl::EvaNuc:*** EVA IS DONE *** F="<<bPDG<<b4M<<
",bB="<<bB<<
", ResNuc="
5639 <<rPDG<<r4M<<
",rB="<<rB<<
G4endl;
5642 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (23) qH=0"<<
G4endl;
5645 if(bB<2) evaHV->push_back(bHadron);
5647 else if(bB==4) evaHV->push_back(bHadron);
5657 ed <<
"Wrong evaporation act: EvaNuc:bB=" << bB
5658 <<
">2 - unexpected evaporated fragment" <<
G4endl;
5659 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0019",
5667 G4cout<<
"G4QNuc::EvaNuc:ResidDibM="<<rM<<
",GSM="<<rGSM<<
",M-GSM="<<rM-rGSM<<
G4endl;
5675 ed <<
"Evaporation below MassShell: <residual> M=" << rM <<
" < GSM="
5677 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0020",
5680 else if(fabs(rM-rGSM)<0.01&&rPDG==90001001) evaHV->push_back(rHadron);
5686 else evaHV->push_back(rHadron);
5689 G4cout<<
"G4QNucleus::EvaporateNucleus: === End of the evaporation attempt"<<
G4endl;
5695 G4cout<<
"*G4QNucleus::EvaporateNucleus: InputHadron4M="<<q4M<<
", PDG="<<thePDG<<
G4endl;
5707 if(totMass+.0001>m1+m2_value)
5710 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (24) qH=0"<<
G4endl;
5720 ed <<
"Chipol->h1+h2 DecIn2 error: tM=" << totMass <<
"-> h1M="
5721 << m1 <<
" + h2M=" << m2_value <<
G4endl;
5722 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0021",
5729 evaHV->push_back(H2);
5734 evaHV->push_back(H1);
5739 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (25) qH=0"<<
G4endl;
5745 ed <<
"Chipolino is under MassShell: M=" << totMass <<
"<" << m1
5746 <<
"+" << m2_value <<
",d=" << m1+m2_value-totMass <<
G4endl;
5747 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0022",
5754 if(fabs(totMass-totM)<0.001||abs(thePDG)-10*(abs(thePDG)/10)>2)
5759 evaHV->push_back(qH);
5761 else if ((thePDG==221||thePDG==331)&&totMass>mPi+mPi)
5764 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (26) qH=0"<<
G4endl;
5774 ed <<
"H->Pi+Pi DecayIn2 error: tM=" << totMass <<
"-> pi+ + pi-"
5776 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0023",
5781 G4cout<<
"G4QNucleus::EvaporateNucleus:(3) PiPlus="<<fq4M<<
G4endl;
5783 evaHV->push_back(H1);
5786 G4cout<<
"G4QNucleus::EvaporateNucleus:(3) PiMinus="<<qe4M<<
G4endl;
5788 evaHV->push_back(H2);
5790 else if ((thePDG==221||thePDG==331)&&totMass>mPi0+mPi0)
5793 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (27) qH=0"<<
G4endl;
5803 ed <<
"H->Pi+Pi DecayIn2 error: tM=" << totMass <<
"-> pi0 + pi0"
5805 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0024",
5810 G4cout<<
"G4QNucleus::EvaporateNucleus:(4) Pi01="<<fq4M<<
G4endl;
5812 evaHV->push_back(H1);
5815 G4cout<<
"G4QNucleus::EvaporateNucleus:(4) Pi02="<<qe4M<<
G4endl;
5817 evaHV->push_back(H2);
5819 else if (totMass>totM)
5822 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (28) qH=0"<<
G4endl;
5832 ed <<
"H*->H+gamma DecIn2 error: tM=" << totMass <<
"->h1M="
5833 << totM <<
"+gam" <<
G4endl;
5834 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0025",
5839 G4cout<<
"G4QNucleus::EvaporateNucleus:(5) tot="<<thePDG<<qe4M<<
G4endl;
5841 evaHV->push_back(H2);
5844 G4cout<<
"G4QNucleus::EvaporateNucleus:(5) GamFortot="<<fq4M<<
G4endl;
5846 evaHV->push_back(H1);
5848 else if (thePDG==111||thePDG==221||thePDG==331)
5851 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (29) qH=0"<<
G4endl;
5861 ed <<
"pi/eta->g+g DecIn2 error: tM=" << totMass
5862 <<
"->gamma + gamma" <<
G4endl;
5863 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0026",
5868 G4cout<<
"G4QNucleus::EvaporateNucleus:(6) gam1="<<qe4M<<
G4endl;
5870 evaHV->push_back(H2);
5873 G4cout<<
"G4QNucleus::EvaporateNucleus:(6) gam2="<<fq4M<<
G4endl;
5875 evaHV->push_back(H1);
5880 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (30) qH=0"<<
G4endl;
5887 ed <<
"Hadron is under MassShell: Nuc=" << thePDG << theQC
5888 <<
", q4M=" << q4M <<
", M=" << totMass <<
" < GSM=" << totM
5889 <<
", 2Pi=" << mPi+mPi <<
", 2Pi0=" << mPi0+mPi0 <<
G4endl;
5890 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0027",
5898 if(!qH)
G4cout<<
"-Warning-G4QNucleus::EvaporateNucleus: (31) qH=0"<<
G4endl;
5904 ed <<
"This is not aNucleus nor aHadron: RN=" << thePDG << theQC
5905 <<
",q4M=" << q4M <<
",qM=" << totMass <<
G4endl;
5906 G4Exception(
"G4QNucleus::EvaporateNucleus()",
"HAD_CHPS_0028",
5914 if(!qH)
G4cout<<
"G4QNucleus::EvaporateNucleus: (20) qH="<<
G4endl;
5919 G4cout<<
"G4QNucleus::EvaporateNucleus: =---=>> End. "<<
G4endl;
5942 G4cout<<
"G4QNuc:DecayIson:QC="<<qQC<<
",M="<<qM<<
",B="<<qBN<<
",S="<<qS<<
",C="<<qC<<
G4endl;
5946 G4cout<<
"--Warning(Upgrade)--G4QNuc::DecIsonuc:FillAsIs,4M="<<q4M<<
",QC="<<qQC<<
G4endl;
5947 evaHV->push_back(qH);
5963 if(-qC==qS && qS==1)
5965 if(fabs(qM-mSigM)<eps)
5967 evaHV->push_back(qH);
5970 else if(qM>mLamb+mPi)
6063 else if(qS>1 && qBN==qS)
6072 else if(!qS && qBN>1)
6081 else G4cout<<
"*?*G4QNuc::DecayIsonucleus: (1) QC="<<qQC<<
G4endl;
6085 if(qS && qS+qC==qBN)
6095 else if(qS && qC<qBN-qS)
6104 else if(qS && qBN==qS)
6108 if(fabs(qM-mSigP)<eps)
6110 evaHV->push_back(qH);
6113 else if(qM>mLamb+mPi)
6118 else if(qM>mNeut+mPi)
6166 else if(qS && qC>qBN-qS)
6169 G4int nPip = qC-qBN;
6199 else if(qBN==qC && qC>1)
6208 else if(qC<=qBN||!qBN)
G4cout<<
"*?*G4QNuc::DecayIsonucleus: (2) QC="<<qQC<<
G4endl;
6214 if(qS) ttM=qS*tMass;
6219 if(fabs(qM-sum)<eps)
6221 f4Mom=q4M*(tfM/sum);
6222 s4Mom=q4M*(tsM/sum);
6223 if(qS) t4Mom=q4M*(ttM/sum);
6225 else if(!qS && (qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom)))
6228 G4cout<<
"***G4QNuc::DecIsonuc:fPDG="<<fPDG<<
"*"<<qBN<<
"(fM="<<fMass<<
")+sPDG="<<sPDG
6229 <<
"*"<<qPN<<
"(sM="<<sMass<<
")"<<
"="<<sum<<
" > TotM="<<qM<<q4M<<qQC<<qS<<
G4endl;
6231 evaHV->push_back(qH);
6237 G4cout<<
"***G4QNuc::DecIsonuc: "<<fPDG<<
"*"<<qBN<<
"("<<fMass<<
")+"<<sPDG<<
"*"<<qPN<<
"("
6238 <<sMass<<
")+Lamb*"<<qS<<
"="<<sum<<
" > TotM="<<qM<<q4M<<qQC<<
G4endl;
6240 evaHV->push_back(qH);
6244 G4cout<<
"G4QNuc::DecayIsonucleus: *DONE* n="<<qPN<<f4Mom<<fPDG<<
", m="<<qPN<<s4Mom<<sPDG
6245 <<
", l="<<qS<<t4Mom<<
G4endl;
6251 for(
G4int ih=0; ih<qBN; ih++)
6254 evaHV->push_back(Hi);
6260 for(
G4int ip=0; ip<qPN; ip++)
6263 evaHV->push_back(Hj);
6269 for(
G4int il=0; il<qS; il++)
6272 evaHV->push_back(Hk);
6278 G4cout <<
"G4QNucleus::DecayIsonucleus: deleted at end - PDG: "
6297 static const G4double mPiN = mPi+mNeut;
6298 static const G4double mPiP = mPi+mProt;
6299 static const G4double dmPiN= mPiN+mPiN;
6300 static const G4double dmPiP= mPiP+mPiP;
6301 static const G4double nnPi = mNeut+mPiN;
6302 static const G4double ppPi = mProt+mPiP;
6303 static const G4double lnPi = mLamb+mPiN;
6304 static const G4double lpPi = mLamb+mPiP;
6305 static const G4double dNeut= mNeut+mNeut;
6306 static const G4double dProt= mProt+mProt;
6307 static const G4double dLamb= mLamb+mLamb;
6308 static const G4double dLaNe= mLamb+mNeut;
6309 static const G4double dLaPr= mLamb+mProt;
6310 static const G4double dSiPr= mSigP+mProt;
6311 static const G4double dSiNe= mSigM+mNeut;
6312 static const G4double dKsPr= mKsiZ+mProt;
6313 static const G4double dKsNe= mKsiM+mNeut;
6322 G4cout<<
"G4QNucl::DecayDibaryon: *Called* PDG="<<qPDG<<
",4Mom="<<q4M<<
", M="<<qM<<
G4endl;
6331 if (qPDG==90003998 && rM>=dmPiP)
6337 else if(qPDG==89998004 && rM>=dmPiN)
6345 else if(qPDG==90000002 && rM>=dNeut)
6352 else if(qPDG==90001001 && rM>=mDeut)
6354 if(fabs(qM-mDeut)<eps)
6356 evaHV->push_back(qH);
6359 else if(mProt+mNeut<rM)
6372 else if(qPDG==91000001 && rM>=dLaNe)
6379 else if(qPDG==91001000 && rM>=dLaPr)
6384 else if(qPDG==89999003 && rM>=nnPi)
6392 else if(qPDG==90002999 && rM>=ppPi)
6396 else if(qPDG==90999002 && rM>=lnPi)
6404 else if(qPDG==91001999 && rM>=lpPi)
6410 else if(qPDG==90999002 && rM>=dSiNe)
6417 else if(qPDG==91001999 && rM>=dSiPr)
6422 else if(qPDG==92000000 && rM>=dLamb)
6429 else if(qPDG==91999001 && rM>=dKsNe)
6436 else if(qPDG==92000999 && rM>=dKsPr)
6441 else if(qPDG!=90002000|| rM<dProt)
6453 G4cerr<<
"***G4QN::DecDiBar: badPDG="<<qPDG<<
" or smallM="<<qM<<
",2mP="<<dProt
6454 <<
",2mN="<<dNeut<<
G4endl;
6465 if(fabs(qM-sum)<eps)
6467 f4Mom=q4M*(fMass/sum);
6468 s4Mom=q4M*(sMass/sum);
6470 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
6472 G4cout<<
"---Warning---G4QN::DecDib:fPDG="<<fPDG<<
"(M="<<fMass<<
")+sPDG="<<sPDG<<
"(M="
6473 <<sMass<<
")="<<sum<<
" >? TotM="<<q4M.
m()<<q4M<<
G4endl;
6476 evaHV->push_back(qH);
6480 G4cout<<
"G4QNucleus::DecayDibaryon:(2) *DONE* f4M="<<f4Mom<<
",fPDG="<<fPDG
6481 <<
", s4M="<<s4Mom<<
",sPDG="<<sPDG<<
G4endl;
6485 evaHV->push_back(H1);
6487 evaHV->push_back(H2);
6494 if(fabs(qM-sum)<eps)
6496 f4Mom=q4M*(fMass/sum);
6497 s4Mom=q4M*(sMass/sum);
6499 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
6501 G4cout<<
"---Warning---G4QN::DecDib:fPDG="<<fPDG<<
"(M="<<fMass<<
")+sPDG="<<sPDG<<
"(M="
6502 <<sMass<<
")"<<
"="<<sum<<
">tM="<<q4M.
m()<<q4M<<
G4endl;
6505 evaHV->push_back(qH);
6509 G4cout<<
"G4QNucleus::DecayDibaryon:(3) *DONE* f4M="<<f4Mom<<
",fPDG="<<fPDG
6510 <<
", s4M="<<s4Mom<<
",sPDG="<<sPDG<<
G4endl;
6513 evaHV->push_back(H1);
6515 evaHV->push_back(H2);
6517 if(fabs(qM-sum)<eps)
6519 f4Mom=q4M*(fMass/sum);
6520 s4Mom=q4M*(sMass/sum);
6522 else if(!
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
6530 ed <<
"General DecayIn2 error: fPDG=" << fPDG <<
"(fM=" << fMass
6531 <<
") + sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum
6532 <<
" >? (DD2,Can't be here) TotM=" << q4M.
m() << q4M <<
G4endl;
6533 G4Exception(
"G4QNucleus::DecayDibaryon()",
"HAD_CHPS_0000",
6537 G4cout<<
"G4QNucl::DecayDibaryon:(4) *DONE* f4M="<<f4Mom<<
",fPDG="<<fPDG
6538 <<
", s4M="<<s4Mom<<
",sPDG="<<sPDG<<
G4endl;
6541 evaHV->push_back(H3);
6543 evaHV->push_back(H4);
6549 if(fabs(qM-sum)<eps)
6551 f4Mom=q4M*(fMass/sum);
6552 s4Mom=q4M*(sMass/sum);
6553 t4Mom=q4M*(tMass/sum);
6555 else if(qM<sum || !
G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom))
6557 G4cout<<
"---Warning---G4QN::DecDib:fPDG="<<fPDG<<
"(M="<<fMass<<
")+sPDG="<<sPDG<<
"(M="
6558 <<sMass<<
")+tPDG="<<tPDG<<
"(tM="<<tMass<<
")="<<sum<<
">TotM="<<q4M.
m()<<
G4endl;
6561 evaHV->push_back(qH);
6565 G4cout<<
"G4QNuc::DecayDibaryon:(5) *DONE* f4M="<<f4Mom<<
",fPDG="<<fPDG<<
", s4M="<<s4Mom
6566 <<
",sPDG="<<sPDG<<
", t4M="<<t4Mom<<
",tPDG="<<tPDG<<
G4endl;
6574 evaHV->push_back(H1);
6576 evaHV->push_back(H2);
6578 evaHV->push_back(H3);
6583 G4cout <<
"G4QNucleus::DecayDiBaryon: deleted at end - PDG: "
6602 static const G4double mPiN = mPi+mNeut;
6603 static const G4double mPiP = mPi+mProt;
6604 static const G4double dmPiN= mPiN+mPiN;
6605 static const G4double dmPiP= mPiP+mPiP;
6606 static const G4double nnPi = mNeut+mPiN;
6607 static const G4double ppPi = mProt+mPiP;
6608 static const G4double lnPi = mLamb+mPiN;
6609 static const G4double lpPi = mLamb+mPiP;
6610 static const G4double dNeut= mNeut+mNeut;
6611 static const G4double dProt= mProt+mProt;
6612 static const G4double dLamb= mLamb+mLamb;
6613 static const G4double dLaNe= mLamb+mNeut;
6614 static const G4double dLaPr= mLamb+mProt;
6615 static const G4double dSiPr= mSigP+mProt;
6616 static const G4double dSiNe= mSigM+mNeut;
6617 static const G4double dKsPr= mKsiZ+mProt;
6618 static const G4double dKsNe= mKsiM+mNeut;
6627 G4cout<<
"G4QNucl::DecayAntiDibar:*Called* PDG="<<qPDG<<
",4Mom="<<q4M<<
", M="<<qM<<
G4endl;
6636 if (qPDG==89996002 && rM>=dmPiP)
6642 else if(qPDG==90001996 && rM>=dmPiN)
6650 else if(qPDG==89999998 && rM>=dNeut)
6657 else if(qPDG==89998999 && rM>=mDeut)
6659 if(fabs(qM-mDeut)<eps)
6661 evaHV->push_back(qH);
6664 else if(mProt+mNeut<rM)
6675 G4cout<<
"--Warning--G4QNucl::DecayAntiDibar:ANTI-DEUTERON is created M="<<rM<<
G4endl;
6678 else if(qPDG==88999999 && rM>=dLaNe)
6685 else if(qPDG==88999999 && rM>=dLaPr)
6690 else if(qPDG==90000997 && rM>=nnPi)
6698 else if(qPDG==89997001 && rM>=ppPi)
6702 else if(qPDG==89000998 && rM>=lnPi)
6710 else if(qPDG==889998001 && rM>=lpPi)
6716 else if(qPDG==89000998 && rM>=dSiNe)
6723 else if(qPDG==88998001 && rM>=dSiPr)
6728 else if(qPDG==88000000 && rM>=dLamb)
6735 else if(qPDG==88000999 && rM>=dKsNe)
6742 else if(qPDG==87999001 && rM>=dKsPr)
6747 else if(qPDG!=89998000|| rM<dProt)
6759 G4cerr<<
"**G4QNuc::DecayAntiDiBar: badPDG="<<qPDG<<
" or smallM="<<qM<<
", 2mP="<<dProt
6760 <<
", 2mN="<<dNeut<<
G4endl;
6771 if(fabs(qM-sum)<eps)
6773 f4Mom=q4M*(fMass/sum);
6774 s4Mom=q4M*(sMass/sum);
6776 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
6778 G4cout<<
"---Warning---G4QN::DecAntiDib:fPDG="<<fPDG<<
"(M="<<fMass<<
")+sPDG="<<sPDG
6779 <<
"(M="<<sMass<<
")="<<sum<<
" >? TotM="<<q4M.
m()<<q4M<<
G4endl;
6782 evaHV->push_back(qH);
6786 G4cout<<
"G4QNucleus::DecayAntiDibaryon:(2) *DONE* f4M="<<f4Mom<<
",fPDG="<<fPDG
6787 <<
", s4M="<<s4Mom<<
",sPDG="<<sPDG<<
G4endl;
6791 evaHV->push_back(H1);
6793 evaHV->push_back(H2);
6800 if(fabs(qM-sum)<eps)
6802 f4Mom=q4M*(fMass/sum);
6803 s4Mom=q4M*(sMass/sum);
6805 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
6807 G4cout<<
"---Warning---G4QN::DecAntiDib:fPDG="<<fPDG<<
"(M="<<fMass<<
")+sPDG="<<sPDG
6808 <<
"(M="<<sMass<<
")"<<
"="<<sum<<
">tM="<<q4M.
m()<<q4M<<
G4endl;
6811 evaHV->push_back(qH);
6815 G4cout<<
"G4QNucleus::DecayAntiDibaryon:(3) *DONE* f4M="<<f4Mom<<
",fPDG="<<fPDG
6816 <<
", s4M="<<s4Mom<<
",sPDG="<<sPDG<<
G4endl;
6819 evaHV->push_back(H1);
6821 evaHV->push_back(H2);
6823 if(fabs(qM-sum)<eps)
6825 f4Mom=q4M*(fMass/sum);
6826 s4Mom=q4M*(sMass/sum);
6828 else if(!
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
6836 ed <<
" General DecayIn2 error: fPDG=" << fPDG <<
"(fM=" << fMass
6837 <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum
6838 <<
" >? (DD2,Can't be here) TotM=" << q4M.
m() << q4M <<
G4endl;
6839 G4Exception(
"G4QNucleus::DecayAntiDibaryon()",
"HAD_CHPS_0000",
6843 G4cout<<
"G4QNucl::DecayAntiDibaryon:(4) *DONE* f4M="<<f4Mom<<
",fPDG="<<fPDG
6844 <<
", s4M="<<s4Mom<<
",sPDG="<<sPDG<<
G4endl;
6847 evaHV->push_back(H3);
6849 evaHV->push_back(H4);
6855 if(fabs(qM-sum)<eps)
6857 f4Mom=q4M*(fMass/sum);
6858 s4Mom=q4M*(sMass/sum);
6859 t4Mom=q4M*(tMass/sum);
6861 else if(qM<sum || !
G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom))
6863 G4cout<<
"-Warning-G4QN::DecAntiDib:fPDG="<<fPDG<<
"(M="<<fMass<<
")+sPDG="<<sPDG<<
"(M="
6864 <<sMass<<
")+tPDG="<<tPDG<<
"(tM="<<tMass<<
")="<<sum<<
">TotM="<<q4M.
m()<<
G4endl;
6867 evaHV->push_back(qH);
6871 G4cout<<
"G4QNuc::DecayAbtiDibaryon:(5) *DONE* f4M="<<f4Mom<<
",fPDG="<<fPDG<<
", s4M="
6872 <<s4Mom<<
",sPDG="<<sPDG<<
", t4M="<<t4Mom<<
",tPDG="<<tPDG<<
G4endl;
6880 evaHV->push_back(H1);
6882 evaHV->push_back(H2);
6884 evaHV->push_back(H3);
6908 G4cout<<
"G4QNuc::DecAntS:QC="<<qQC<<
",S="<<qS<<
",B="<<qB<<
",C="<<qP<<
",4M="<<q4M<<
G4endl;
6910 G4int qN = qB-qP-qS;
6917 ed <<
"not an Anti Strange Nucleus: QC=" << qQC <<
",S=" << qS <<
",B="
6918 << qB <<
",4M=" << q4M <<
G4endl;
6919 G4Exception(
"G4QNucleus::DecayAntiStrange()",
"HAD_CHPS_0000",
6950 G4int sS=(aS-dPN)/2;
6953 if(qP>=sS && qN>=bS)
6989 G4int sS=(aS-dNP)/2;
6991 if(qN>=bS && qP>=sS)
7009 G4int qPDG=90000000+(qP-n2)*1000+(qN-n1);
7012 G4cout<<
"G4QNucleus::DecayAnStran:nK0="<<n1<<
",nK+="<<n2<<
", nucM="<<nucM<<
G4endl;
7016 if(qP>=-qS) m2_value=-qS;
7017 else if(qP>0) m1=-qS-qP;
7018 G4int sPDG=90000000+(qP-m2_value)*1000+(qN-m1);
7020 if(mucM+m1*mK+m2_value*mK0<nucM+n1*mK+n2*mK0)
7028 G4cout<<
"G4QNucleus::DecayAnStran: n1="<<n1<<
", n2="<<n2<<
", nM="<<nucM<<
G4endl;
7034 if(n2==1 && mK+nucM>qM+.0001)
7038 qPDG=90000000+qP*1000+qN-1;
7050 if(n1==1 && mK0+nucM>qM+.0001)
7054 qPDG=90000000+(qP-1)*1000+qN;
7060 G4int naPDG=90000000+(qP-1)*1000+qN;
7065 naPDG=90000000+qP*1000+qN-1;
7069 G4cout<<
"G4QNucleus::DecayAnStran:M="<<qM<<
",kM="<<k1M<<
"+nM="<<nucM<<
"="<<k1M+nucM
7070 <<
",m="<<kaM<<
"+n="<<naM<<
"="<<kaM+naM<<
G4endl;
7076 if(sum>qM+eps && n1==1)
7078 G4int naPDG=90000000+(qP-1)*1000+qN;
7084 naPDG=90000000+qP*1000+qN-1;
7101 if(fabs(qM-sum)<eps)
7103 f4Mom=q4M*(n1M/sum);
7104 s4Mom=q4M*(nucM/sum);
7106 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
7109 G4cout<<
"--Warning--G4QNuc::DASt:AsItIs, H="<<qQC<<q4M<<qM<<
" < sum="<<sum<<
"=(F)"
7110 <<nucM<<
"+(kK)"<<n1M<<
G4endl;
7112 evaHV->push_back(qH);
7116 G4cout<<
"G4QNuc::DecAntiS: nK+N "<<n1<<
"*K="<<k1PDG<<f4Mom<<
",N="<<qPDG<<s4Mom<<
G4endl;
7121 for(
G4int i1=0; i1<n1; i1++)
7124 evaHV->push_back(H1);
7130 G4cout<<
"G4QNucleus::DecAntiStr:*** After EvaporateNucleus nH="<<evaHV->size()<<
G4endl;
7141 if(fabs(qM-sum)<eps)
7143 f4Mom=q4M*(n1M/sum);
7144 s4Mom=q4M*(n2M/sum);
7145 t4Mom=q4M*(nucM/sum);
7147 else if(qM<sum || !
G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom))
7149 G4cout<<
"---Warning---G4QN::DASt:nPDG="<<qPDG<<
"(M="<<nucM<<
")+1="<<k1PDG<<
"(M="<<k1M
7150 <<
")+2="<<k2PDG<<
"(M="<<k2M<<
")="<<nucM+n1*k1M+n2*k2M<<
">tM="<<qM<<q4M<<
G4endl;
7151 evaHV->push_back(qH);
7155 G4cout<<
"G4QNuc::DecAntiS:*DONE* nPDG="<<qPDG<<
",1="<<f4Mom<<
",2="<<s4Mom<<
",n="<<t4Mom
7161 for(
G4int i1=0; i1<n1; i1++)
7164 evaHV->push_back(H1);
7167 for(
G4int i2=0; i2<n2; i2++)
7170 evaHV->push_back(H2);
7179 G4cout <<
"G4QNucleus::DecayAntiStrange: deleted at end - PDG: "
7185 G4cout<<
"G4QNucleus::DecayAntiStrange: ===> End of DecayAntiStrangness"<<
G4endl;
7201 G4cout<<
"G4QNuc::DecayMultyBaryon: *Called* PDG="<<qPDG<<
",4M="<<q4M<<qQC<<
G4endl;
7206 G4int totN=totBN-totS-totC;
7214 else if(totC==totBN)
7219 else if(totS!=totBN)
7225 ed <<
"Can not decay this PDG Code: PDG=" << qPDG <<
G4endl;
7226 G4Exception(
"G4QNucleus::DecayMultyBaryon()",
"HAD_CHPS_0000",
7236 ed <<
"Unknown PDG code of the MultiBaryon: PDG=" << qPDG <<
G4endl;
7237 G4Exception(
"G4QNucleus::DecayMultyBaryon()",
"HAD_CHPS_0001",
7241 if(totBN==1) evaHV->push_back(qH);
7247 if(fabs(qM-sum)<eps)
7254 G4cout<<
"---Warning---G4QNucl::DecayMultyBar:fPDG="<<fPDG<<
"(fM="<<fMass<<
")*2="<<sum
7255 <<
" > TotM="<<q4M.
m()<<q4M<<
G4endl;
7258 evaHV->push_back(qH);
7262 G4cout<<
"G4QNucleus::DecMulBar:*DONE* fPDG="<<fPDG<<
",f="<<f4Mom<<
",s="<<s4Mom<<
G4endl;
7266 evaHV->push_back(H1);
7268 evaHV->push_back(H2);
7276 if(fabs(qM-sum)<eps)
7284 G4cout<<
"---Warning---G4QNuc::DecayMultyBaryon: fPDG="<<fPDG<<
"(fM="<<fMass<<
")*3 = "
7285 <<3*fMass<<
" >? TotM="<<q4M.
m()<<q4M<<
G4endl;
7288 evaHV->push_back(qH);
7292 G4cout<<
"G4QNuc::DecMBar:*DONE*, fPDG="<<fPDG<<
",f="<<f4Mom<<
",s="<<s4Mom<<
",t="
7297 evaHV->push_back(H1);
7299 evaHV->push_back(H2);
7301 evaHV->push_back(H3);
7310 G4cout<<
"*G4QNul::DecMulBar:SplitMultiBar inEqParts M="<<totBN<<
"*"<<f4Mom.
m()<<
G4endl;
7311 G4cout<<
"G4QNucleus::DecMultyBaryon: *DONE* fPDG="<<fPDG<<
", f="<<f4Mom<<
G4endl;
7314 for(
G4int h=0; h<totBN; h++)
7317 evaHV->push_back(H1);
7323 G4cout <<
"G4QNucleus::DecayMultyBaryon: deleted at end - PDG: "
7342 G4cout<<
"G4QNuc::DecayAlphaDiN: *Called* PDG="<<qPDG<<
",4M="<<q4M<<
G4endl;
7346 G4int sPDG = 90002002;
7350 if(fabs(qM-mHel6)<eps)
7352 evaHV->push_back(qH);
7355 else if(mNeut+mNeut+mAlph<qM)
7366 ed <<
"Cannot decay excited He6 with this mass: M(He6=" << mHel6 <<
")="
7367 << qM <<
"<" << mNeut+mNeut+mAlph <<
G4endl;
7368 G4Exception(
"G4QNucleus::DecayAlphaDiN()",
"HAD_CHPS_0000",
7372 else if(qPDG!=90004002)
7378 ed <<
"Can not decay this PDG Code: PDG=" << qPDG <<
G4endl;
7379 G4Exception(
"G4QNucleus::DecayAlphaDiN()",
"HAD_CHPS_0001",
7386 if(fabs(qM-sum)<eps)
7388 f4Mom=q4M*(fMass/sum);
7390 t4Mom=q4M*(sMass/sum);
7392 else if(qM<sum || !
G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom))
7394 G4cout<<
"---Warning---G4QNuc::DecayAlphaDiN:fPDG="<<fPDG<<
"(M="<<fMass<<
")*2+mAlpha = "
7395 <<sum<<
" >? TotM="<<qM<<q4M<<
", d="<<sum-qM<<
G4endl;
7398 evaHV->push_back(qH);
7402 G4cout<<
"G4QNuc::DecAl2N: fPDG="<<fPDG<<
",f="<<f4Mom<<
",s="<<s4Mom<<
",t="<<t4Mom<<
G4endl;
7406 evaHV->push_back(H1);
7408 evaHV->push_back(H2);
7410 evaHV->push_back(H3);
7428 G4cout<<
"G4QNucleus::DecayAlphaBar: *Called* PDG="<<qPDG<<
",4M="<<q4M<<qQC<<
G4endl;
7434 if ( ( (!totS && !totC) || totC == totBN || totS == totBN)
7436 else if(qPDG==92001002||qPDG==92002001||qPDG==91003001||qPDG==91001003||qPDG==93001001)
7437 evaHV->push_back(qH);
7438 else if(qPDG==92000003||qPDG==92003000||qPDG==93000002||qPDG==93002000)
7449 else if(qPDG==93000002)
7456 else if(qPDG==93002000)
7468 if(fabs(qM-sum)<eps)
7470 f4Mom=q4M*(tfM/sum);
7471 s4Mom=q4M*(tsM/sum);
7473 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
7475 G4cout<<
"--Warning--G4QNuc::DecAlB:fPDG="<<fPDG<<
"(M="<<fMass<<
")*2="<<2*fMass<<
",s="
7476 <<sPDG<<
"(sM="<<sMass<<
")*3="<<3*sMass<<
"="<<sum<<
">M="<<q4M.
m()<<q4M<<
G4endl;
7479 evaHV->push_back(qH);
7483 G4cout<<
"G4QNucleus::DecAlB:*DONE*, fPDG="<<fPDG<<f4Mom<<
",sPDG="<<sPDG<<s4Mom<<
G4endl;
7488 evaHV->push_back(H1);
7489 evaHV->push_back(H1);
7492 evaHV->push_back(H2);
7493 evaHV->push_back(H2);
7494 evaHV->push_back(H2);
7496 else if(qPDG==90004001||qPDG==90001004)
7498 G4int fPDG = 90002001;
7513 if(fabs(qM-sum)<eps)
7515 f4Mom=q4M*(fMass/sum);
7516 s4Mom=q4M*(sMass/sum);
7519 else if(qM<sum || !
G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom))
7521 G4cout<<
"--Warning--G4QNuc::DecAlB:fPDG="<<fPDG<<
",M="<<fMass<<
",sPDG="<<sPDG<<
",sM="
7522 <<sMass<<
",2sM+fM="<<2*sMass+fMass<<
" > TotM="<<q4M.
m()<<q4M<<
G4endl;
7525 evaHV->push_back(qH);
7529 G4cout<<
"G4QNucl::DecAlB: *DONE*, f="<<fPDG<<f4Mom<<
", s="<<sPDG<<s4Mom<<t4Mom<<
G4endl;
7533 evaHV->push_back(H1);
7535 evaHV->push_back(H2);
7537 evaHV->push_back(H3);
7539 else if(qPDG==94000001||qPDG==94001000||qPDG==91000004||qPDG==91004000)
7550 else if(qPDG==91000004)
7557 else if(qPDG==91004000)
7567 if(fabs(qM-sum)<eps)
7569 f4Mom=q4M*((fMass+fMass)/sum);
7570 s4Mom=q4M*(sMass/sum);
7572 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
7574 G4cout<<
"--Warning--G4QNucl::DecAlphBar:fPDG="<<fPDG<<
"(2*fM="<<fMass<<
")*2="
7575 <<2*fMass<<
",sPDG="<<sPDG<<
"(sM="<<sMass<<
" > TotM="<<q4M.
m()<<q4M<<
G4endl;
7578 evaHV->push_back(qH);
7582 G4cout<<
"G4QNuc::DecAlphaB: *DONE*, fPDG="<<fPDG<<f4Mom<<
",sPDG="<<sPDG<<s4Mom<<
G4endl;
7587 evaHV->push_back(H1);
7588 evaHV->push_back(H1);
7589 evaHV->push_back(H1);
7590 evaHV->push_back(H1);
7592 evaHV->push_back(H2);
7594 else if(qPDG==90003002||qPDG==90002003||qPDG==91002002)
7596 G4int fPDG = 90002002;
7605 else if(qPDG==9100202)
7610 else if(qPDG!=90002003)
7612 evaHV->push_back(qH);
7620 G4cout<<
"***Corrected***G4QNuc::DecayAlphaBar:fPDG="<<fPDG<<
"(fM="<<fMass<<
")+ sPDG="
7621 <<sPDG<<
"(sM="<<sMass<<
")="<<fMass+sMass<<
" > TotM="<<qM<<q4M<<
G4endl;
7630 if(fabs(qM-sum)<eps)
7632 f4Mom=q4M*(fMass/sum);
7633 s4Mom=q4M*(sMass/sum);
7635 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
7637 G4cout<<
"--Warning--G4QNuc::DecAlphaBar:fPDG="<<fPDG<<
"(fM="<<fMass<<
")+sPDG="<<sPDG
7638 <<
"(sM="<<sMass<<
")="<<fMass+sMass<<
"="<<sum<<
" > TotM="<<q4M.
m()<<q4M<<
G4endl;
7641 evaHV->push_back(qH);
7645 G4cout<<
"G4QNucl::DecAlBar:*DONE*a4M="<<f4Mom<<
",s4M="<<s4Mom<<
",sPDG="<<sPDG<<
G4endl;
7649 evaHV->push_back(H1);
7651 evaHV->push_back(H2);
7653 else G4cout<<
"---Warning---G4QNucleus::DecayAlphaBar: Unknown PDG="<<qPDG<<
G4endl;
7657 G4cout <<
"G4QNucleus::DecayAlphaBar: deleted at end - PDG: "
7677 ed <<
"Not Be8 state decais in 2 alphas: qPDG=" << qPDG <<
G4endl;
7678 G4Exception(
"G4QNucleus::DecayAlphaAlpha()",
"HAD_CHPS_0000",
7684 G4cout<<
"G4QNucleus::DecayAlAl: *Called* PDG="<<qPDG<<
",M="<<qM<<q4M<<
">"<<aaGSM<<
G4endl;
7690 G4int sPDG = 90004004;
7696 if(fabs(qM-sum)<eps)
7698 f4Mom=q4M*(fMass/sum);
7699 s4Mom=q4M*(sMass/sum);
7701 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
7703 G4cout<<
"---Warning---G4QNuc::DecayAlphaAlpha:gPDG="<<fPDG<<
"(gM="<<fMass<<
")+PDG="
7704 <<sPDG<<
"(sM="<<sMass<<
")="<<sum<<
" > TotM="<<q4M.
m()<<q4M<<
G4endl;
7707 evaHV->push_back(qH);
7711 G4cout<<
"G4QNucleus::DecayAlphaAlpha: *DONE* gam4M="<<f4Mom<<
", aa4M="<<s4Mom<<
G4endl;
7714 evaHV->push_back(H1);
7718 G4int fPDG = 90002002;
7719 G4int sPDG = 90002002;
7724 if(fabs(qM-sum)<eps)
7726 f4Mom=q4M*(fMass/sum);
7729 else if(qM<sum || !
G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))
7731 G4cout<<
"---Warning---G4QNucl::DecayAlphaAlpha:fPDG="<<fPDG<<
"(fM="<<fMass<<
")*2="<<sum
7732 <<
" > TotM="<<q4M.
m()<<q4M<<
G4endl;
7735 evaHV->push_back(qH);
7739 G4cout<<
"G4QNucleus::DecayAlphaAlpha: *DONE* fal4M="<<f4Mom<<
", sal4M="<<s4Mom<<
G4endl;
7743 evaHV->push_back(H1);
7745 evaHV->push_back(H2);