62G4double G4QNucleus::nucleonDistance=.8*fermi;
63G4double G4QNucleus::WoodSaxonSurf=.545*fermi;
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();
257 for(
G4int i=0; i<nn; ++i)
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();
294 for(
G4int i=0; i<nn; ++i)
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();
318 for(
G4int i=0; i < nn; ++i)
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;
429 for (
G4int in=1; in<256; in++) probVect[in]=0.;
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);
1406 G4double sm=GSResPP+mProt+mProt+SPPBarr;
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;
3145 for (
G4int index=0; index<mCand; index++)
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;
3694 static const G4double r0sq=0.8133*fermi*fermi;
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;
3768 static const G4double constofpmax=hbarc*pow(3.*pi2,third);
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;
4001 G4double C = A*D/pi/std::log(1.+B);
4021 G4int nt = Tb.size();
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;
4034 static const G4double dfermi=fermi/10.;
4035 static const G4double sfermi=fermi*fermi;
4038 G4int nb =
static_cast<int>(bf);
4040 G4int nt = Tb.size();
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);
6234 else if(qS && (qM<sum || !
G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom)))
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);
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4QCandidate * > G4QCandidateVector
std::vector< G4QHadron * > G4QHadronVector
ostream & operator<<(ostream &lhs, G4QNucleus &rhs)
G4ThreeVector G4RandomDirection()
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
double dot(const Hep3Vector &) const
void SetPreProbability(G4double prep)
void SetPossibility(G4bool choice)
void SetDenseProbability(G4double prep)
G4int GetBaryonNumber() const
G4int GetStrangeness() const
G4int GetSPDGCode() const
G4int GetZNSPDGCode() const
G4ThreeVector Get3Momentum() const
G4LorentzVector Get4Momentum() const
void SetNFragments(const G4int &nf)
G4int GetBaryonNumber() const
G4double GetEnergy() const
G4bool DecayIn2(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom)
const G4ThreeVector & GetPosition() const
G4bool DecayIn3(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom, G4LorentzVector &t4Mom)
void SetQC(const G4QContent &newQC)
void SetPosition(const G4ThreeVector &aPosition)
void SetBindingEnergy(G4double aBindE)
G4LorentzVector theMomentum
void SetQPDG(const G4QPDGCode &QPDG)
G4int GetNFragments() const
void Set4Momentum(const G4LorentzVector &aMom)
G4QPDGCode GetQPDG() const
G4int GetStrangeness() const
G4double BindingEnergy(const G4double &cZ=0, const G4double &cA=0, G4double dZ=0., G4double dA=0.)
void SimpleSumReduction(G4ThreeVector *vectors, G4ThreeVector sum)
void InitByPDG(G4int newPDG)
G4double GetRelativeDensity(const G4ThreeVector &aPosition)
void DecayAlphaBar(G4QHadron *dB, G4QHadronVector *oHV)
std::pair< G4double, G4double > ChooseImpactXandY(G4double maxImpact)
G4double GetRelOMDensity(const G4double &r2)
G4double GetRelWSDensity(const G4double &r)
G4double GetOuterRadius()
G4double CoulombBarGen(const G4double &rZ, const G4double &rA, const G4double &cZ, const G4double &cA)
G4int RandomizeBinom(G4double p, G4int N)
G4QNucleus operator-=(const G4QNucleus &rhs)
void SubtractNucleon(G4QHadron *pNucleon)
G4int HadrToNucPDG(G4int hPDG)
void InitCandidateVector(G4QCandidateVector &theQCandidates, G4int nM=45, G4int nB=72, G4int nC=117)
void DecayAntiDibaryon(G4QHadron *dB, G4QHadronVector *oHV)
void DoLorentzBoost(const G4LorentzVector &theBoost)
void Increase(G4int PDG, G4LorentzVector LV=G4LorentzVector(0., 0., 0., 0.))
G4double GetDeriv(const G4ThreeVector &point)
void DecayIsonucleus(G4QHadron *dB, G4QHadronVector *oHV)
void DecayAntiStrange(G4QHadron *dB, G4QHadronVector *oHV)
void DecayAlphaAlpha(G4QHadron *dB, G4QHadronVector *oHV)
void DecayMultyBaryon(G4QHadron *dB, G4QHadronVector *oHV)
G4double GetRadius(const G4double maxRelativeDenisty=0.5)
G4QContent GetQCZNS() const
const G4QNucleus & operator=(const G4QNucleus &right)
void ChooseFermiMomenta()
G4double CoulBarPenProb(const G4double &CB, const G4double &E, const G4int &C, const G4int &B)
void DecayDibaryon(G4QHadron *dB, G4QHadronVector *oHV)
static void SetParameters(G4double fN=.1, G4double fD=.05, G4double cP=4., G4double mR=1., G4double nD=.8 *CLHEP::fermi)
void PrepareCandidates(G4QCandidateVector &theQCandidates, G4bool piF=false, G4bool gaF=false, G4LorentzVector LV=G4LorentzVector(0., 0., 0., 0.))
void ActivateBThickness()
G4int NucToHadrPDG(G4int nPDG)
G4bool ReduceSum(G4ThreeVector *vectors, G4ThreeVector sum)
void DecayAlphaDiN(G4QHadron *dB, G4QHadronVector *oHV)
G4double GetThickness(G4double b)
std::vector< G4double > const * GetBThickness()
G4bool EvaporateBaryon(G4QHadron *h1, G4QHadron *h2)
G4double GetFermiMomentum(G4double density)
G4int UpdateClusters(G4bool din)
void InitByQC(G4QContent newQC)
G4double FissionCoulombBarrier(const G4double &cZ, const G4double &cA, G4double dZ=0., G4double dA=0.)
G4QNucleus operator+=(const G4QNucleus &rhs)
G4QNucleus operator*=(const G4int &rhs)
G4double GetDensity(const G4ThreeVector &aPos)
void DoLorentzContraction(const G4LorentzVector &B)
void DoTranslation(const G4ThreeVector &theShift)
G4double CoulombBarrier(const G4double &cZ=1, const G4double &cA=1, G4double dZ=0., G4double dA=0.)
void EvaporateNucleus(G4QHadron *hA, G4QHadronVector *oHV)
G4double GetGSMass() const
G4QContent GetQuarkContent() const
void InitByQCode(G4int QCode)
G4double GetNuclMass(G4int Z, G4int N, G4int S)
void ConvertPDGToZNS(G4int PDG, G4int &z, G4int &n, G4int &s)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
std::ostringstream G4ExceptionDescription