75 : q4Mom(q4M), valQ(qQCont), theWorld(0), phot4M(ph4M), f2all(0), rEP(0.), rMo(0.)
78 G4cout<<
"G4Quasmon:Constructor:QC="<<qQCont<<
",Q4M="<<q4M<<
",photonE="<<ph4M.
e()<<
G4endl;
81 G4cout<<
"**>G4Q:Con:(1),T="<<Temperature<<
",S="<<SSin2Gluons<<
",E="<<EtaEtaprime<<
G4endl;
83 if(phot4M.
e()>0.) q4Mom+=phot4M;
93 status = right.status;
95 G4int nQH = right.theQHadrons.size();
96 if(nQH)
for(
G4int ih=0; ih<nQH; ih++)
99 theQHadrons.push_back(curQH);
101 theWorld = right.theWorld;
102 phot4M = right.phot4M;
103 nBarClust = right.nBarClust;
106 G4int nQC = right.theQCandidates.size();
107 if(nQC)
for(
G4int iq=0; iq<nQC; iq++)
110 theQCandidates.push_back(curQC);
120 G4cout<<
"G4Quasmon::Copy-Constructor: ***CALLED*** E="<<right->theEnvironment<<
G4endl;
122 q4Mom = right->q4Mom;
125 status = right->status;
127 G4int nQH = right->theQHadrons.size();
129 G4cout<<
"G4Quasmon::Copy-Constructor:nQH="<<nQH<<
G4endl;
131 if(nQH)
for(
G4int ih=0; ih<nQH; ih++)
134 G4cout<<
"G4Quasmon:Copy-Constructor:H#"<<ih<<
",QH="<<right->theQHadrons[ih]<<
G4endl;
137 theQHadrons.push_back(curQH);
139 theWorld = right->theWorld;
140 phot4M = right->phot4M;
141 nBarClust = right->nBarClust;
144 G4int nQC = right->theQCandidates.size();
146 G4cout<<
"G4Quasmon:Copy-Constructor: nCand="<<nQC<<
G4endl;
148 if(nQC)
for(
G4int iq=0; iq<nQC; iq++)
151 G4cout<<
"G4Quasmon:Copy-Constructor:C#"<<iq<<
",QC="<<right->theQCandidates[iq]<<
G4endl;
154 theQCandidates.push_back(curQC);
156 f2all = right->f2all;
160 G4cout<<
"G4Quasmon:Copy-Constructor: >->-> DONE <-<-<"<<
G4endl;
167 G4cout<<
"G4Quasmon::Destructor before theQCandidates delete"<<
G4endl;
171 G4cout<<
"G4Quasmon::Destructor before theQHadrons"<<
G4endl;
173 for_each(theQHadrons.begin(), theQHadrons.end(),
DeleteQHadron());
179 G4double G4Quasmon::Temperature=180.;
180 G4double G4Quasmon::SSin2Gluons=0.1;
181 G4double G4Quasmon::EtaEtaprime=0.3;
182 G4bool G4Quasmon::WeakDecays=
false;
183 G4bool G4Quasmon::ElMaDecays=
true;
194 Temperature=temperature;
209 status = right.status;
211 G4int iQH = theQHadrons.size();
212 if(iQH)
for(
G4int jh=0; jh<iQH; jh++)
delete theQHadrons[jh];
214 G4int nQH = right.theQHadrons.size();
215 if(nQH)
for(
G4int ih=0; ih<nQH; ih++)
218 theQHadrons.push_back(curQH);
220 theWorld = right.theWorld;
221 phot4M = right.phot4M;
222 nBarClust = right.nBarClust;
225 G4int iQC = theQCandidates.size();
226 if(iQC)
for(
G4int jq=0; jq<iQC; jq++)
delete theQCandidates[jq];
227 theQCandidates.clear();
228 G4int nQC = right.theQCandidates.size();
229 if(nQC)
for(
G4int iq=0; iq<nQC; iq++)
232 theQCandidates.push_back(curQC);
245 static const G4int NUCPDG = 90000000;
246 static const G4int MINPDG = 80000000;
248 static const G4double BIG = 1000000.;
249 static const G4double BIG2 = BIG*BIG;
277 static const G4double diPiM= mPi0 + mPi0;
278 static const G4double PiNM = mPi + mNeut;
280 static const G4double mPi02= mPi0* mPi0;
286 static const G4double mP2 = mProt*mProt;
296 G4cout<<
"G4Quasmon::HadrQ:*==>>START QUASMON HADRONIZATION<<==*, aP="<<addPhoton<<
",Env="
302 if(nQuasms==-1) first=
true;
303 else G4cout<<
"G4Quasmon::HadrQ: Negative #of Quasmons n="<<nQuasms<<
G4endl;
311 G4cout<<
"G4Q::HQ: PionAtRest, addP="<<addPhoton<<
", momP="<<momPhoton<<
G4endl;
317 else if(addPhoton>0.) gaF=
true;
328 G4cout<<
"G4Quasmon::HadrQ:CHIPSWorld initialized with nP="<<nP<<
",nM="<<nMesons<<
G4endl;
330 if (nP<34) nMesons = 9;
331 else if(nP<51) nMesons = 18;
332 else if(nP<65) nMesons = 27;
333 else if(nP<82) nMesons = 36;
335 if (nP<45) nBaryons = 16;
336 else if(nP<59) nBaryons = 36;
337 else if(nP<76) nBaryons = 52;
340 G4cout<<
"G4Quasmon:HadrQ: Init Candidates:"<<theEnvironment<<
",n="<<theQCandidates.size()
341 <<
",nMesons="<<nMesons<<
",nBaryons="<<nBaryons<<
",nClusters="<<nClusters<<
G4endl;
345 G4cout<<
"G4Quasmon:HadrQ:CandidatesAreInitialized,n="<<theQCandidates.size()<<
",nMesons="
346 <<nMesons<<
", nBaryons="<<nBaryons<<
", nClusters="<<nClusters<<
G4endl;
348 if(!status||q4Mom==zeroLV)
351 G4cout<<
"G4Q::HQ:NOTHING-TO-DO: Q4M="<<q4Mom<<
", QEnv="<<theEnvironment<<
G4endl;
356 ed <<
"OverheadPhoton for theZeroQuasmon: Q4M=" << q4Mom <<
",status="
357 << status <<
", phE=" << addPhoton <<
G4endl;
399 G4int oldNH=theQHadrons.size();
405 while(theEnvironment.
GetPDG()==NUCPDG || start)
409 G4int nHd=theQHadrons.size();
410 if(nHd)
for(
int ih=0; ih<nHd; ih++) ccSum+=theQHadrons[ih]->
GetCharge();
413 G4cerr<<
"*G4Q::HQ:C"<<cSum<<
",c="<<ccSum<<
",E="<<theEnvironment<<
",Q="<<valQ<<
G4endl;
414 G4cerr<<
":G4Q::HQ:oldE="<<oldEnv<<
"oldQ="<<oldCQC<<
",oN="<<oldNH<<
",N="<<nHd<<
G4endl;
415 if(nHd)
for(
int h=0; h<nHd; h++)
430 if(fabs(qM2)<.0001 || tmpEq<=tmpPq)
439 G4cout<<
"G4Q::HQ:Quasmon is gamma, Q4M="<<q4Mom<<
",E="<<theEnvironment<<
G4endl;
442 FillHadronVector(gamH);
450 G4cout<<
"G4Q::HQ:Quasmon is nothing, Q4M="<<q4Mom<<
",E="<<theEnvironment<<
G4endl;
457 else q4Mom.
setE(tmpPq*1.00001);
460 G4double qurF = quasM/(tmpEq-tmpPq);
464 G4cout<<
"G4Q::HQ: Quasm="<<q4Mom<<
",qM="<<quasM<<
",qQC="<<valQ<<
G4endl;
475 G4int dmaxActEnv=maxActEnv+maxActEnv;
498 G4cout<<
"G4Q::HQ: ePDG="<<envPDG<<
",eM="<<envM<<
",eQC="<<envQC<<
G4endl;
505 G4int totNeut=totN.GetN();
509 G4cout<<
"G4Q::HQ: tN="<<totN<<
",tGSM="<<totM<<
",tQC="<<totQC<<
G4endl;
521 resNM = resNN.GetMZNS();
522 resNPDG= resNN.GetPDG();
528 G4int totPDG=totN.GetPDG();
530 G4cout<<
"G4Q::HQ: totPDG="<<totPDG<<
",totM="<<totMass<<
",rPDG="<<resNPDG<<
G4endl;
536 G4cerr<<
"---Warning---G4Q::HQ: *Boost* tot4M="<<tot4M<<
", E-p="<<totEn-totMo<<
G4endl;
541 G4cerr<<
"G4Q::HQ: *Boost* E-p shift is corrected to "<<emodif<<
G4endl;
542 tot4M.
setE(totMo+emodif+.01*accuracy);
554 G4cout<<
"G4Q::HQ: iniPDG="<<iniPDG<<
", Z="<<iniP<<
", N="<<iniN<<
", S="<<iniS<<
G4endl;
560 if(iniBN<2||envA>0) iniRM=0.;
561 G4cout<<
"G4Q::HQ: iniRN="<<iniRN<<
", iniRM="<<iniRM<<
", iniQM="<<iniQM<<
G4endl;
568 if(envPDG==NUCPDG) bndQM=iniQM;
571 G4cout<<
"G4Q::HQ:mQ="<<quasM<<valQ<<bndQM2<<
",nQ="<<nOfQ<<
",Env="<<envM
572 <<envQC<<
",Q+E="<<quen<<
",tM="<<totPDG<<totQC<<totM<<
"<"<<totMass<<
G4endl;
578 if (bQ) s_value = 3*bQ + 2;
579 if (tQ> s_value) cQ.
DecQAQ((tQ-s_value)/2);
581 G4int rsPDG = cQ.GetSPDGCode();
582 G4cout<<
"G4Q::HQ:eN="<<envN<<
",eZ="<<envZ<<
",Q="<<rsPDG<<cQ<<
",piF="<<piF<<
",gaF="<<gaF
590 ModifyInMatterCandidates();
611 if(excE > diPiM) qCountMax=(
G4int)(excE/mPi0);
625 if(envA>0&&envA<19) pCountMax=36/envA;
629 while(kCount<kCountMax&&kCond)
635 if(excE>mPi0) miM2=mPi02;
682 G4int totN_value=totBN-totZ;
683 if(totN_value>0&&totBN>1)
686 G4double delN=tmpNN.GetMZNS()+mNeut-totM;
690 G4double delEN=envNN.GetMZNS()+mNeut-envM;
691 if(delEN>delN) delN=delEN;
694 if(delN<minK) minK=delN;
700 G4double delP=tmpPN.GetMZNS()+mProt-totM+proCB;
704 G4double delEP=envPN.GetMZNS()+mProt-envM+proCB;
705 if(delEP>delP) delP=delEP;
708 if(delP<minK) minK=delP;
710 if(totN_value>0&&totZ>0&&totBN>2)
714 G4double delD=tmpDN.GetMZNS()+mDeut-totM+proCB;
715 if(envN>0&&envZ>0&&envA>2)
718 G4double delED=envDN.GetMZNS()+mDeut-envM+proCB;
719 if(delED>delD) delD=delED;
722 if(delD<minK) minK=delD;
724 if(totN_value>1&&totZ>0&&totBN>3)
728 G4double delT=tmpTN.GetMZNS()+mTrit-totM+proCB;
729 if(envN>1&&envZ>0&&envA>3)
732 G4double delET=envTN.GetMZNS()+mTrit-envM+proCB;
733 if(delET>delT) delT=delET;
736 if(delT<minK) minK=delT;
738 if(totN_value>0&&totZ>1&&totBN>3)
742 G4double delR=tmpRN.GetMZNS()+mHel3-totM+proCB;
743 if(envN>0&&envZ>1&&envA>3)
746 G4double delER=envRN.GetMZNS()+mHel3-envM+proCB;
747 if(delER>delR) delR=delER;
750 if(delR<minK) minK=delR;
752 if(totN_value>1&&totZ>1&&totBN>4)
756 G4double delA=tmpAN.GetMZNS()+mAlph-totM+proCB;
757 if(envN>1&&envZ>1&&envA>4)
760 G4double delEA=envAN.GetMZNS()+mAlph-envM+proCB;
761 if(delEA>delA) delA=delEA;
764 if(delA*qurF<minK) minK=delA*qurF;
781 if(minK<0. || minK+minK > quasM) minK=0.;
803 G4double kpow=
static_cast<double>(nOfQ-2);
805 G4double xmi=(momPhoton-addPhoton)*quasM;
815 kMom=(1.-(1.-xmi)*pow(rn,1./kpow))*quasM/2;
833 if(cost<-1.) cost=-1.;
835 G4cout<<
"G4Q::HQ:**PHOTON out of Q**k="<<kMom<<
",ct="<<cost<<
",QM="<<quasM<<
G4endl;
872 if(!miM2) miM2=(minK+minK)*quasM;
873 if(qM2<.0001) kMom=0.;
874 else kMom = GetQPartonMomentum(maxK,miM2);
889 G4double cQM2=qM2-(kMom+kMom)*quasM;
899 if(!
G4QHadron(q4Mom).RelDecayIn2(k4Mom,cr4Mom,dir4M,cost,cost))
902 G4cerr<<
"*G4Q::HQ:PB,M="<<quasM<<
",cM="<<cQM<<
",c="<<cost<<
",F="<<gintFlag<<
G4endl;
905 if(addPhoton&&gintFlag)
909 if(qM2>0.) quasM = sqrt(qM2);
912 if(qM2<-.0001)
G4cerr<<
"--Warning-- G4Q::HQ:Phot.M2="<<qM2<<
" Cor to 0"<<
G4endl;
920 if(addPhoton&&gintFlag)
930 G4cout<<
"G4Q::HQ:(PhBack) k="<<kMom<<
",k4M="<<k4Mom<<
",ct="<<cost<<
",gF="<<gintFlag
937 G4int totCand = theQCandidates.size();
939 G4cout<<
"G4Q::HQ: ****>>K-ITERATION #"<<kCount<<
", k="<<kMom<<k4Mom<<
G4endl;
945 if(cPDG==90000001||cPDG==90001000||cPDG==91000000||cPDG<MINPDG)
949 if(cPDG==90000001 || cPDG==90001000 || cPDG==90000002 || cPDG==90001001)
950 G4cout<<
"G4Q::HQ:pos="<<poss<<
",cPDG="<<cPDG<<
",iQC="<<iniQChg<<
G4endl;
963 if(cPDG==90000001 || cPDG==90001000 || cPDG==90000002 || cPDG==90001001)
969 if(cMs) kMin=(cfM*cfM-cMs*cMs)/(cMs+cMs);
978 if(cPDG==90000001 || cPDG==90001000 || cPDG==90000002 || cPDG==90001001)
979 G4cout<<
"G4Q::HQ:i="<<
index<<
",cPDG="<<cPDG<<
",k="<<kMom<<
","<<kLS<<
">kMin="
980 <<kMin<<
",bM="<<bnM<<
",rEP="<<rEP<<
",dR="<<dR<<
",kCo="<<kCond<<
G4endl;
996 G4cout<<
"G4Q::HQ:***PHOTON OK***k="<<k4Mom<<
",Q="<<q4Mom<<
",PhE="<<addPhoton<<
G4endl;
1000 G4cout<<
"G4Q::HQ:Select="<<kMom<<
",ki="<<minK<<
",ka="<<maxK<<
",k="<<k4Mom<<kLS<<
G4endl;
1002 CalculateHadronizationProbabilities(excE,kMom,k4Mom,piF,gaF,first);
1009 G4int nCandid = theQCandidates.size();
1015 if(nCandid) maxP = theQCandidates[nCandid-1]->GetIntegProbability();
1017 G4cout<<
"G4Q::HQ:***RANDOMIZE CANDIDATEs***a#OfCand="<<nCandid<<
",maxP="<<maxP<<
G4endl;
1022 if(status == 4)
G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0001",
1030 G4cout<<
"G4Q::HQ:qCB="<<qCB<<
",Z="<<iniP<<
",B="<<iniBN<<
",eE="<<excE<<
G4endl;
1033 G4bool qenv=envPDG!=90000000&&envPDG!=90000002&&envPDG!=90002000&&iniPDG!=90002000
1036 G4cout<<
"G4Q::HQ: qCB="<<qCB<<
",pCB="<<pCB<<
",cond="<<qenv<<
",N="<<nQuasms<<
G4endl;
1039 if(!first&&excE>qCB&&envM+iniQM<totMass&&nQuasms==1&&qenv&&
G4UniformRand()<pCB)
1045 FillHadronVector(resNuc);
1048 G4cout<<
"G4Q::HQ: outQM="<<oQ4Mom.
m()<<oQ4Mom<<
",GSQM="<<iniQM<<
G4endl;
1050 if(nQuasms==1) qEnv =
G4QNucleus(env4M,envPDG);
1054 FillHadronVector(envH);
1060 G4cout<<
"G4Q::HQ: envM="<<envM<<
"=="<<eM<<
", envT="<<env4M.e()-eM<<dif<<
G4endl;
1068 G4cout<<
"G4Q::HQ:Q2E:E="<<theEnvironment<<
",valQ="<<valQ<<
",tot4M="<<tot4M<<
G4endl;
1081 else if(envPDG==NUCPDG && quasM>iniQM)
1084 G4cout<<
"***G4Q::HQ: Emergency Decay in Gamma/Pi0 + Residual GSQuasmon"<<
G4endl;
1089 if(quasM>mPi0+iniQM)
1097 if(sum>0. && fabs(quasM-sum)<eps)
1099 r4Mom=q4Mom*(gamM/sum);
1100 s4Mom=q4Mom*(iniQM/sum);
1108 ed <<
"(E=0)G/Pi0+GSQ decay error: Q=" << q4Mom << quasM
1109 <<
"->g/pi0(M=" << gamM <<
")+GSQ=" << iniPDG <<
"(M="
1110 << iniQM <<
")=" << sum <<
", d=" << sum-quasM <<
G4endl;
1111 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0002",
1115 G4cout<<
"G4Q::HQ:=== 0 ===>HadrVec, Q="<<q4Mom<<quasM<<
"->g/pi0("<<gamPDG<<
")="
1116 <<r4Mom<<gamM<<
"+GSQ="<<iniPDG<<r4Mom<<iniQM<<
G4endl;
1119 FillHadronVector(curHadr2);
1121 FillHadronVector(curHadr1);
1125 G4cout<<
"***G4Q::HQ:dE="<<excE<<
">minK="<<minK<<
",Env="<<theEnvironment<<
",k="<<kLS
1126 <<
",Q="<<valQ<<quasM<<
", nQ="<<nQuasms<<
G4endl;
1128 qEnv=theEnvironment;
1170 G4cout<<
"G4Q::HQ: fp="<<fprob<<
",QM="<<quasM<<
",QQC="<<valQ<<
",k="<<kMom<<
G4endl;
1172 while(pCount<pCountMax&&pCond)
1176 G4cout<<
"G4Q::HQ:***>New p-Attempt#"<<pCount<<
",pMax="<<pCountMax<<
",hsfl=0"<<
G4endl;
1179 while(theQCandidates[i]->GetIntegProbability() < totP) i++;
1183 ed <<
"Too big number of the candidate: Cand#"<< i <<
" >= Tot#"<< nCandid <<
G4endl;
1191 if(!sPDG&&theEnvironment.
GetPDG()!=NUCPDG&&totBN>1&&totMass>totM&&totS>=0)
1195 ed <<
"Why Fail? (1): NEED-EVAP-1:Q=" << q4Mom << valQ <<
",En="
1196 << theEnvironment <<
G4endl;
1199 qEnv=theEnvironment;
1205 ed <<
"DecayPartSelection,Evaporation: sPDG=0,E=" << theEnvironment
1206 <<
",B=" << totBN <<
",S=" << totS <<
G4endl;
1219 if(sMass+rMass>quasM)
1221 if(totBN>1&&totMass>totM&&totS>=0)
1225 ed <<
" Why Fail? (2): NEED-EVAP-2:Q=" << q4Mom << valQ <<
",E="
1226 << theEnvironment <<
G4endl;
1227 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0006",
1230 qEnv=theEnvironment;
1236 ed <<
"VirtChipo Can't EvapNucleus: QM=" << quasM <<
"<S=" << sMass
1237 <<
"+R=" << rMass <<
"=" << sMass+rMass <<
",tB=" << totBN
1238 <<
",tS=" << totS <<
",tM=" << totMass <<
">minM=" << totM <<
G4endl;
1239 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0007",
1248 ed <<
"Why Fail ? (0): NEED-EVAP-0:Q=" << q4Mom << valQ <<
",En="
1249 << theEnvironment <<
G4endl;
1252 qEnv=theEnvironment;
1259 if (rPDG>MINPDG&&rPDG!=NUCPDG)
1261 rMass=theEnvironment.
GetMZNS();
1263 valQ +=theEnvironment.
GetQC();
1268 else if(rPDG!=NUCPDG&&totBN>1&&totMass>totM&&totS>=0)
1279 if(dMpT>dMnT)dMnT=dMpT;
1281 ed <<
"Why Fail ? (3): NEED-EVAP3:s=" << sPDG <<
",Q=" << q4Mom
1282 << valQ <<
",r=" << rPDG <<
G4endl;
1285 qEnv=theEnvironment;
1288 else if(rPDG==NUCPDG)
1291 G4cout<<
"G4Quasmon::HadronizeQuasm:SafatyDecayIn PI0/GAM, rPDG="<<rPDG<<
G4endl;
1293 if(totMass-totM>mPi0)
1307 G4cout<<
"G4Q::HQ:hsflagTRUE,s="<<sPDG<<
","<<sMass<<
",r="<<rPDG<<
","<<rMass<<
G4endl;
1316 G4cout<<
"G4Q::HQ:hsfl="<<hsflag<<
", sPDG="<<sPDG<<
", i="<<i<<
G4endl;
1319 if(sPDG>MINPDG&&sPDG!=NUCPDG)
1326 G4cerr<<
"***G4Quasmon::HadronizeQ:P="<<theQCandidates[i]->GetIntegProbability()
1327 <<
",nC="<<nParCandid<<
",pP="<<sppm<<
",QM="<<quasM<<
",QC="<<valQ;
1328 for(
int ipp=0; ipp<nParCandid; ipp++)
1331 "NoParentClust forTheFragment");
1338 G4cout<<
"G4Q::HQ:p#ip="<<ip<<
",f#i="<<i<<
",tP="<<totPP<<
",sP="<<sppm<<
G4endl;
1344 pQC = pQPDG.GetQuarkContent();
1351 loli = parCluster->
GetLow();
1359 sQC = sQPDG.GetQuarkContent();
1362 else sMass = sQPDG.GetMass();
1366 G4cout<<
"G4Q::HQ:valQ="<<valQ<<
"+transQ="<<transQC<<
"("<<pPDG<<
" to "<<sPDG
1367 <<
") = curQ="<<curQ<<
",InvBinding="<<delta<<
",pM="<<pMass<<pQC<<
G4endl;
1373 sQC=theQCandidates[i]->GetQC();
1374 sMass = theQCandidates[i]->GetNBMass();
1378 G4cout<<
"G4Q::HQ: hsfl="<<hsflag<<
", valQ="<<valQ<<
"-sQ="<<sQC<<
",sM="<<sMass
1379 <<
",C="<<theQCandidates[i]->GetPDGCode()<<
",Q="<<curQ<<
",M2="<<sM2<<
G4endl;
1386 if(resTN.GetA()>0) sCB=totN.CoulombBarrier(sQC.GetCharge(),sQC.GetBaryonNumber());
1389 G4cout<<
"G4Q::HQ:rQC="<<resNQC<<
",rM="<<resTNM<<
",sM="<<sMass<<
",CB="<<sCB<<
G4endl;
1400 if(sPDG>MINPDG&&sPDG!=NUCPDG) RNQC-=pQC;
1401 if(envA-pBaryn>bEn) RNQC=curQ+bEnQC;
1407 G4cout<<
"**G4Q::HQ:*KinematicTotal*, PDG="<<RNPDG<<curQ<<
", QC="<<RNQC<<
G4endl;
1416 G4cout<<
"G4Q::HQ:M="<<bEn<<
",A="<<envA<<
",B="<<pBaryn<<
",N="<<minN<<
G4endl;
1423 G4cout<<
"*G4Q::HQ:EnvironmentA="<<envA<<
" < SecondaryFragmentA="<<pBaryn<<
G4endl;
1436 }
else if(!rqPDG||rQQ<-1) {
1438 G4cerr<<
"*G4Q::HQ:*** ResidualQuasmon *** PDG="<<rqPDG<<curQ<<
",Q="<<rQQ<<
G4endl;
1441 minSqT=10000000000.;
1442 minSqB=10000000000.;
1447 if(sPDG<MINPDG&&envPDG>MINPDG&&envPDG!=NUCPDG)
1454 G4cout<<
"G4Q::HadQ:Z="<<rqZ<<
",N="<<rqN<<
",S="<<rqS<<
",eZ="<<envZ<<
",eN="<<envN
1455 <<
",eS="<<envS<<
",ePDG="<<envPDG<<
",eM="<<envM<<
",tM="<<qpeM<<
G4endl;
1461 G4cout<<
"G4Q::HQ:rPDG="<<rqPDG<<curQ<<
",minT="<<minT<<
",minSqT="<<minSqT
1462 <<
",hsfl="<<hsflag<<
G4endl;
1468 if ( ( (sPDG < MINPDG && envPDG > MINPDG && envPDG != NUCPDG) ||
1469 (sPDG > MINPDG && sPDG != NUCPDG && envPDG > pPDG)
1470 ) && ( (rqPDG > MINPDG && rqPDG != NUCPDG) ||
1471 rqPDG==2112 || rqPDG==2212 || rqPDG==3122
1493 G4cout<<
"G4Q::HQ:***VacuumFragmentation** M="<<newT<<
",rM="<<rtM<<rtQC
1494 <<
",eM="<<envM<<
",mM="<<minT<<
G4endl;
1501 if(envA-pBaryn>bEn) reQC=bEnQC;
1507 G4cout<<
"G4Q::HQ:reQC="<<reQC<<
",rtQC="<<rtQC<<
",eA="<<envA<<
",pB="<<pBaryn
1508 <<
",bE="<<bEn<<bEnQC<<
G4endl;
1514 if (envA-pBaryn > bEn) newT=rtM-mbEn;
1517 G4cout<<
"G4Q::HQ:NuclFrM="<<newT<<
",r="<<rtM<<rtQC<<
",e="<<envM<<envQC<<
",p="
1518 <<pMass<<pQC<<
",re="<<reM<<reQC<<
",exEn="<<totMass-rtM-sMass<<
G4endl;
1521 if(minT<newT) newT=minT;
1526 G4cout<<
"G4Q::HQ:rq="<<rqPDG<<
",miT="<<minSqT<<
",miB="<<minSqB<<
",M="<<rtM<<
G4endl;
1531 ed <<
"MinResMass can't be calculated: minSqT=0(!), curQ=" << curQ <<
G4endl;
1536 if (sPDG > MINPDG && sPDG != NUCPDG) {
1538 G4cout<<
"G4Q::HQ: BoundM="<<pMass<<
",FreeM="<<sMass<<
",QM="<<quasM<<
G4endl;
1545 G4cout<<
"G4Q::HQ:Q("<<quasM<<
")->k("<<k4Mom<<
")+CRQ("<<cr4Mom.
m()<<
")"<<
G4endl;
1553 G4cout<<
"***G4Q::HQ:FailedToFind FragmM:"<<ccM2<<
"<"<<frM2<<
",M="<<pMass<<
"+k="
1554 <<k4Mom<<
"="<<sqrt(ccM2)<<cc4Mom<<
" < fM="<<sMass<<
",miK="<<minK<<
G4endl;
1560 tot4Mom=q4Mom+cl4Mom;
1561 cc4Mom=k4Mom+cl4Mom;
1566 G4cout<<
"G4Q::HQ:hsflagTRUE*NuclBINDING,ccM2="<<ccM2<<
"<frM2="<<frM2<<
G4endl;
1573 G4cout<<
"G4Q::HQ:***NUCLEAR BINDING***ccM2="<<ccM2<<
" > frM2="<<frM2<<
G4endl;
1582 G4cout<<
"G4Q::HQ:cM2="<<crMass2<<
"="<<rEP*(rEP-rMo-rMo)<<
",h="<<hili<<
",l="
1585 if(hili<loli) hili=loli;
1587 G4double fpqp2=
static_cast<double>(npqp2);
1594 if(sPDG==90001001) dM=2.25;
1595 G4cout<<
"G4Q::HQ:Is xE="<<excE<<
" > sM="<<sMass<<
"-pM="<<pMass<<
"-dM="<<dM
1596 <<
" = "<<sMass-pMass-dM<<
G4endl;
1599 G4cout<<
"G4Q::HQ: must totM="<<totMass<<
" > rTM="<<resTNM<<
"+sM="<<sMass<<
" = "
1602 if(resTNM && totMass<resTNM+sMass)
1605 G4cout<<
"***G4Quasmon::HadronizeQuasmon:***PANIC#1***TotalDE="<<excE<<
"< bE="
1606 <<sMass-pMass-dM<<
", dM="<<dM<<
", sM="<<sMass<<
", bM="<<pMass<<
G4endl;
1610 qEnv=theEnvironment;
1616 if(envA-pBaryn>bEn) tmpEQ=bEnQC;
1620 G4cout<<
"G4Q::HQ:eQC="<<envQC<<
",pQC="<<pQC<<
",rEnvM="<<tmpNM<<
",hsfl="<<hsflag
1636 G4double cta=1.-(dex/(1.-pow(loli,pw))-pMass)/kLS;
1637 if(cta>1.0001)
G4cerr<<
"Warn-G4Q::HQ: cost_max="<<cta<<
">1.CorHadrProb"<<
G4endl;
1638 G4double kap_a=ex/(1.+kLS*(1.-cta)/pMass);
1639 G4double cti=1.-(dex/(1.-pow(hili,pw))-pMass)/kLS;
1640 if(cti<-1.0001)
G4cerr<<
"Warn-G4Q::HQ: cost_i="<<cti<<
"<-1.CorHadrProb"<<
G4endl;
1641 G4double kap_i=ex/(1.+kLS*(1.-cti)/pMass);
1642 G4double q_lim=(tmpTM2-tCRN.
m2())/(tcEP+tcEP);
1643 if(cti>cta+.0001)
G4cerr<<
"**G4Q::HQ:ci="<<cti<<
">ca="<<cta<<
".CorHPro"<<
G4endl;
1644 G4cout<<
"G4Q::HQ:qi="<<kap_i<<
",ci="<<cti<<
",a="<<kap_a<<
",ca="<<cta<<
",e="<<ex
1645 <<
",q="<<q_lim<<
",S="<<tmpTM*tmpTM<<
",R2="<<tCRN.
m2()<<
","<<tcEP<<
G4endl;
1668 while(qCount<qCountMax&&qCond)
1674 G4double ctkk=1.-(dex/(1.-
z)-pMass)/kLS;
1676 if(qCount)
G4cout<<
"G4Q::HQ:qC="<<qCount<<
",ct="<<ctkk<<
",M="<<pMass<<
",z="
1677 <<z<<
",zl="<<pow(loli,pw)<<
",zh="<<pow(hili,pw)<<
",dE="
1678 <<totMass-totM<<
",bE="<<sMass-pMass<<
G4endl;
1681 G4cout<<
"G4Q::HQ:ct="<<ctkk<<
",pM="<<pMass<<
",z="<<z<<
",zl="<<pow(loli,pw)
1682 <<
",zh="<<pow(hili,pw)<<
",ex="<<ex<<
",li="<<loli<<
",hi="<<hili<<
G4endl;
1684 if(abs(ctkk)>1.00001)
1687 G4cerr<<
"***G4Q:HQ:ctkk="<<ctkk<<
",ex="<<ex<<
",z="<<z<<
",pM="<<pMass
1688 <<
",kLS="<<kLS<<
",hi="<<hili<<
",lo="<<loli<<
",n="<<npqp2<<
G4endl;
1691 if(ctkk> 1.)ctkk= 1.;
1692 if(ctkk<-1.)ctkk=-1.;
1695 G4double ctc=(cen*ctkk-kLS)/(cen-kLS*ctkk);
1700 else if(ctc<-1.) ctc=-1.;
1704 if(!
G4QHadron(cc4Mom).RelDecayIn2(kp4Mom, fr4Mom, k4Mom, ctc, ctc))
1707 ed <<
"Can't dec ColClust(Fr+kap): c4M=" << cc4Mom <<
",sM="
1708 << sMass <<
",ct=" << ctc <<
G4endl;
1709 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0012",
1719 G4double kappa=ex/(1.+kLS*(1.-ctkk)/pMass);
1720 G4cout<<
"G4Q::HQ:>>ColDec>>CF("<<ccM<<
")->F("<<sMass<<
")+q"<<kp4Mom<<
"="
1721 <<kappa<<
",hsfl="<<hsflag<<
G4endl;
1725 rQ4Mom=cr4Mom+kp4Mom;
1727 reTNM2=retN4Mom.
m2();
1740 if(fQM2>=minSqT && reTNM2>=tmpTM2 && fr4Mom.e()>=sCBE)
1751 G4cout<<
"G4Q::HQ:Attemp#"<<qCount<<
".Yes.M2="<<fQM2<<
">"<<minSqT<<
G4endl;
1763 G4cout<<
"G4Q::HQ:Attempt#"<<qCount<<
",NO.M2="<<fQM2<<
"<"<<minSqT<<
G4endl;
1772 if(reTNM2<tmpTM2 && fQM2>MEMrQM2 && fr4Mom.e()>=sCBE)
1787 else if(fr4Mom.e()<sCBE&&fr4Mom.e()>MEMsCBE&&reTNM2>=tmpTM2&&fQM2>MEMrQM2)
1834 G4cout<<
"G4Q::HadQ:RQ("<<rQ4Mom.m()<<
")=C("<<cr4Mom.
m()<<
")+q"<<kp4Mom<<
G4endl;
1847 G4cout<<
"G4Q::HQ:A="<<envA<<
",B="<<pBaryn<<
",Q="<<rQ4Mom.m()<<
","<<piF<<
G4endl;
1850 tot4Mom-=rQ4Mom+fr4Mom;
1852 G4cout<<
"G4Q::HQ:t4M="<<tot4Mom<<
",hsfl="<<hsflag<<
".Is kt="<<kt<<
">"<<minSqB
1853 <<
" or kn="<<kn<<
">"<<minSqN<<
"? m2="<<m2_value<<
", sPDG="<<sPDG<<
G4endl;
1859 if(kn<minSqN && ku<minSqT)
1876 G4cout<<
"G4Q::HQ:**hsflag=1** No, sPDG="<<sPDG<<
", kt="<<kt<<
"<"<<minSqB
1877 <<
" or kn="<<kn<<
"<"<<minSqN<<
G4endl;
1881 else G4cout<<
"G4Q::HQ:YES,t="<<kt<<
">"<<minSqB<<
",n="<<kn<<
">"<<minSqN<<
G4endl;
1887 kt = (quasM-dk)*(quasM-sM2/dk);
1889 if(kt>0.) rQM=sqrt(kt);
1892 if(!
G4QHadron(q4Mom).DecayIn2(fr4Mom, rQ4Mom))
1895 ed <<
"Can't dec Quasmon in Fr+rQuas: q4M=" << q4Mom <<
", sM="
1896 << sMass <<
", rQM=" << rQM <<
G4endl;
1900 if(envPDG>MINPDG&&envPDG!=NUCPDG)
1904 if(envA>bEn) TCRN+=bEn4M;
1916 G4cout<<
"G4Q::HQ:k="<<kMom<<
".F:"<<kt<<
">"<<minSqB<<
",N:"<<kn<<
">"<<minSqN<<
" &tM="
1917 <<tM<<
">rtM="<<rtM<<
" & hsfl="<<hsflag<<
" to avoid decay R+S="<<sPDG<<
G4endl;
1924 if(kt>minSqB+.01 && tM>rtM && !hsflag)
1950 G4cout<<
"G4Q::HQ:YES forFragment it's big enough:kn="<<kn<<
">"<<minSqN<<
G4endl;
1959 G4cout<<
"G4Q::HQ:NO,hsfl=1,kt="<<kt<<
"<"<<minSqB<<
" or M="<<tM<<
"<"<<rtM<<
G4endl;
1963 G4cout<<
"G4Q::HQ:******>>rM="<<rMass<<
",sqM2="<<sqrt(m2_value)<<
",hsfl="<<hsflag<<
G4endl;
1965 rMass=sqrt(m2_value);
1971 if(rPDG==111&&sPDG!=111&&rrn>.5) rPDG=221;
1973 G4int aPDG = abs(rPDG);
1979 ed <<
"Unidentifiable residual Hadron: rQ =" << curQ <<
", rPDG=" << rPDG
1980 <<
"(b=" << rb <<
") + sPDG=" << sPDG <<
"(sM=" << sMass <<
")" <<
G4endl;
1981 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0014",
1987 if(rPDG==221 || rPDG==331) rcMass=mPi0;
1992 if(sPDG<MINPDG&&envPDG==NUCPDG)rFl=
G4UniformRand()<bs*bs*bs/m3_value;
1994 G4cout<<
"G4Q::HQ: sPDG="<<sPDG<<
", hsflag="<<hsflag<<
", rPDG="<<rPDG<<curQ<<
",rM="
1995 <<rMass<<
",rb="<<rb<<
",F="<<rFl<<
",v="<<m2_value<<
", bs="<<bs<<
G4endl;
1997 if(!hsflag&&rFl&&rPDG&& rPDG!=10 && rb<2 && aPDG!=1114 && aPDG!=2224 && aPDG!=3334)
2003 if(rPDG && rPDG!=10)
2005 if (rPDG== 3122) rPDG= 3212;
2006 else if(rPDG==-3122) rPDG=-3212;
2007 if(rPDG>0)repPDG=rPDG+2;
2009 if(repPDG>0)redPDG=repPDG+2;
2010 else redPDG=repPDG-2;
2011 if(redPDG>0)refPDG=redPDG+2;
2012 else refPDG=redPDG-2;
2013 if(refPDG>0)regPDG=refPDG+2;
2014 else regPDG=refPDG-2;
2017 G4cout<<
"G4Q::HQ:QuasM="<<quasM<<valQ<<
")->H("<<sPDG<<
")+R("<<rPDG<<
")"<<
",rp="
2018 <<repPDG<<
",rd="<<redPDG<<
",rf="<<refPDG<<
",rg="<<regPDG<<
G4endl;
2024 sB = sqrt((resM2+repM2)/2.);
2025 G4double pB = sqrt((repM2+redM2)/2.);
2026 G4double dB = sqrt((redM2+refM2)/2.);
2028 if(!cB) fB= sqrt((refM2+
G4QPDGCode(regPDG).GetMass2())/2.);
2030 G4double rM = GetRandomMass(repPDG,dif);
2031 G4double dM = GetRandomMass(redPDG,dif);
2032 G4double fM = GetRandomMass(refPDG,dif);
2034 G4cout<<
"G4Q::HQ: rM="<<rM<<
",rMa="<<rMass<<
",sB="<<sB<<
"(bQ="<<bQ<<
"),pB="<<pB
2035 <<
",dM="<<dM<<
",dB="<<dB<<
",fM="<<fM<<
",fB="<<fB<<
G4endl;
2037 if(((rM>0 && rMass<pB && rMass>sB) || (dM>0 && rMass>pB && rMass<dB) ||
2038 (fM>0 && rMass>dB && rMass<fB)) && theEnvironment.
GetPDG()==NUCPDG)
2040 if (rMass>pB && rMass<dB && dM>0)
2045 else if(rMass>dB && rMass<fB && dM>0)
2051 G4cout<<
"G4Q::HQ:s="<<sPDG<<
",Q=>rM="<<rMass<<
"(minQ="<<rPDG<<curQ<<
")+sB="
2054 if(quasM<rM+sMass &&(sPDG==221||sPDG==331))
2062 if(fabs(quasM-sum)<eps)
2064 r4Mom=q4Mom*(rM/sum);
2065 s4Mom=q4Mom*(sMass/sum);
2073 ed <<
"H+Res Decay failed: rPD=" << repPDG <<
"(rM=" << rMass
2074 <<
")+sPD=" << sPDG <<
"(sM=" << sMass <<
"), Env="
2075 << theEnvironment <<
G4endl;
2076 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0015",
2080 G4cout<<
"G4Q::HQ:=== 1 ===> HadronVec, Q="<<q4Mom<<
" -> s4M="<<s4Mom<<
"("
2081 <<sPDG<<
"), r4M="<<r4Mom<<
"("<<repPDG<<
")"<<
G4endl;
2085 FillHadronVector(curHadr1);
2087 if (sPDG==2112) sPDG=90000001;
2088 else if(sPDG==2212) sPDG=90001000;
2089 else if(sPDG==3122) sPDG=91000000;
2092 FillHadronVector(curHadr2);
2094 qEnv=theEnvironment;
2101 fdul = rFl && rPDG!=10;
2106 if (kn > minSqN && ku > minSqT)
2118 G4cout<<
"G4Q::HQ:P-Attempt#"<<pCount<<
" *Yes* sPDG="<<sPDG<<
",ku="<<ku<<
">"
2119 <<minSqT<<
" || kn="<<kn<<
">"<<minSqN<<
G4endl;
2132 G4cout<<
"G4Q::HQ:P-Attempt#"<<pCount<<
",No. ku="<<ku<<
"<"<<minSqT<<
" or kn="<<kn
2133 <<
"<"<<minSqN<<
" or E="<<fr4Mom.e()<<
"<"<<sCBE<<
G4endl;
2145 if (kn < minSqN && ku < minSqT)
2155 if(ku < minSqT && ku > PMEMktM2)
2183 G4cout<<
"G4Q::HQ:RemTheBest rPDG="<<rPDG<<
",sPDG="<<sPDG<<
",kt="<<kt<<
G4endl;
2214 G4cout<<
"G4Q::HQ:RemTheFirst rPDG="<<rPDG<<
",sPDG="<<sPDG<<
",kt="<<kt<<
G4endl;
2250 G4cout<<
"G4Q::HQ:>rPDG="<<rPDG<<curQ<<
",sPDG="<<sPDG<<
",kt="<<kt<<
",F="<<fprob
2251 <<
",totQC="<<totQC<<
",sQC="<<sQC<<
G4endl;
2257 if(rPDG==111&&sPDG!=111&&rrr>.5) rPDG=221;
2258 if(rPDG==221&&sPDG!=221&&sPDG!=331&&rrr<.5) rPDG=111;
2265 ed <<
"Unidentifiable residual Hadron: Q=" << curQ <<
",r=" << rPDG
2266 <<
"+s=" << sPDG <<
"(sM=" << sMass <<
")" <<
G4endl;
2269 if(rPDG==221||rPDG==331) reMass=mPi0;
2274 if ( ( ( (sPDG < MINPDG && envPDG > MINPDG && envPDG != NUCPDG) ||
2275 (sPDG > MINPDG && sPDG!=NUCPDG && envPDG > pPDG)
2277 ) || iniBN > 1 || rPDG == 10
2281 G4cout <<
"G4Q::HQ:Is hsfl="<<hsflag<<
" or fdul="<<fdul<<
" or [rM="<<rMass<<
"<"<<reMass
2282 <<
" + "<<aMass<<
" or rM2="<<reTNM2<<
" < miM2="<<tmpTM2<<
" and ePDG="<<envPDG
2283 <<
">pPDG="<<pPDG<<
"] to fail?"<<
G4endl;
2287 (sPDG < MINPDG && rMass < reMass+aMass) ||
2288 (sPDG > MINPDG && envPDG > pPDG && reTNM2 < tmpTM2) ||
2294 G4cout<<
"G4Q::HQ: Yes(No), hsf="<<hsflag<<
",sPDG="<<sPDG<<
",pM="<<pMass<<
",Env="
2295 <<envPDG<<
",QM="<<quasM<<valQ<<
", fpr="<<fprob<<
G4endl;
2298 if(hsflag) rMass=rQPDG.
GetMass();
2299 if(sPDG>MINPDG&&sPDG!=NUCPDG)
2312 ed <<
"Why Fail? (5): NEEDS-EVAP-5,Q=" << q4Mom << valQ <<
",QEnv="
2313 << theEnvironment <<
G4endl;
2316 qEnv=theEnvironment;
2326 if(tmpTM>retNM) tmpT=
G4QNucleus(tmpTQ,retN4Mom);
2338 if(rCB < 0.) rCB=0.;
2339 G4int sB=sQPDG.GetBaryNum();
2341 if(sCB < 0.) sCB=0.;
2354 G4cout<<
"G4Q::HQ:tM="<<totMass<<
",totQC="<<totQC<<
",rtQC="<<tmpTQ<<
",pQC="<<pQC
2355 <<
",sB="<<sB<<
",resB="<<tmpT.GetA()<<
G4endl;
2358 if(nQuasms==1 && tmpNM+rMass+rCB+sMass+sCB < totMass &&
2369 G4cout<<
"G4Q::HQ: *YES*,tM="<<ctM<<
"="<<totMass<<
",fM="<<cfM<<
"="<<sMass<<
G4endl;
2372 if(fabs(totMass-sum)<eps)
2374 re4M=tot4M*(tmpNM/sum);
2375 rq4M=tot4M*(rMass/sum);
2376 fr4M=tot4M*(sMass/sum);
2381 ed <<
"DecayIn Frag+ResQ+ResE failed: Decay (" << totMass <<
") in Fragm("
2382 << sMass <<
")+ResQ("<< rMass <<
")+ResEnv("<< tmpNM <<
")="<< sum <<
G4endl;
2386 FillHadronVector(resQH);
2390 FillHadronVector(envH);
2397 G4cout<<
"**G4Q::HQ:(3)**KeepEnvironmentMoving**, nQ="<<nQuasms<<
G4endl;
2401 FillHadronVector(candH);
2407 else if(nQuasms==1 && tmpTM+sMass+sCB < totMass)
2412 G4cout<<
"G4Q::HQ:(2)*KeepEnvironmentMoving*,nQ="<<nQuasms<<
",Env="<<qEnv<<
G4endl;
2420 G4cout<<
"G4Q::HQ: rTotM="<<retN4Mom.
m()<<
" >? GSM="<<tmpTM<<
",d4M="<<d4M<<
",dQC="
2424 FillHadronVector(candHadr);
2428 G4cout<<
"G4Q::HQ:sM="<<sMass<<
"="<<frM<<
", fT="<<fr4Mom.e()-frM<<
",dif24M="<<dif2
2434 else if(totBN>1 &&totMass>totM &&totS>=0&&envPDG>MINPDG&&envPDG!=NUCPDG)
2443 ed <<
"Why Fail?(6): ProductMasses>totalMass: EV-6: TotEVAPORATION:s=" << sPDG
2444 <<
",T=" << kinE <<
",RM=" << retN4Mom.
m() <<
"<" << tmpTM <<
",tQC="
2445 << transQC <<
",E=" << excE <<
",sM=" << sumM <<
">tM=" << totMass <<
",nQ="
2450 G4cout<<
"G4Q::HQ:Q="<<q4Mom<<quasM<<
",E="<<theEnvironment<<
",P="<<phot4M<<
G4endl;
2452 qEnv=theEnvironment;
2455 else if(totBN==1 && nQuasms==1)
2458 G4cout<<
"G4Q::HQ:tB=1,nQ=1,Z="<<totZ<<
",S="<<totS<<totQC<<
",M="<<totMass<<
G4endl;
2462 G4int nucPDG = 2212;
2468 if(!totZ&&totMass>mLamb+mPi0)
2475 else if(abs(totZ)==1&&totMass>mLamb+mPi)
2480 if(totZ>0) piPDG = 211;
2486 ed <<
"Pi + Lambda decay error:Z=" << totZ <<
",S=" << totS
2487 << totQC <<
",tM=" << totMass <<
G4endl;
2488 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0020",
2494 if(!totZ&&totMass>mNeut+mK0)
2501 else if(totZ==2&&totMass>mProt+mK)
2511 else if(totZ==1&&totMass>=mNeut+mK)
2521 ed <<
"K + Nucleon decay error: Z=" << totZ <<
",S=" << totS
2522 << totQC <<
",tM=" << totMass <<
G4endl;
2523 G4Exception(
"G4Quasmon::HadronizeQuasmon",
"HAD_CHPS_0021",
2528 else if(totMass>PiNM&&!totS)
2535 else if(!totZ&&totMass>mNeut+mPi0)
2549 else if(totZ==1&&totMass>mProt+mPi0)
2569 ed <<
"Pi + Nucleon decay error: Z=" << totZ <<
",B=" << totBN
2570 <<
",E=" << envQC <<
",Q=" << valQ <<
G4endl;
2571 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0022",
2582 else if(totZ<0||totZ>1)
2585 ed <<
"Photon+Nucleon decay error: Z=" << totZ <<
",B=" << totBN
2586 <<
",E=" << envQC <<
",Q=" << valQ <<
G4endl;
2587 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0023",
2594 if(fabs(totMass-sum)<eps)
2596 pi4M=tot4M*(piM/sum);
2597 nuc4M=tot4M*(nucM/sum);
2602 ed <<
"Gam/Pi/K+N decay error: T=" << tot4M << totMass
2603 <<
"->gam/pi/K(" << piM <<
")+N=" << nucPDG <<
"(" << nucM
2604 <<
")=" << sum <<
G4endl;
2605 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0024",
2609 G4cout<<
"G4Q::HQ:T="<<tot4M<<totMass<<
"->GPK="<<piPDG<<pi4M<<
"+B="<<nucPDG<<nuc4M
2613 FillHadronVector(piH);
2615 FillHadronVector(nucH);
2621 else G4cout<<
"***G4Q::HQ: B="<<totBN<<
",tM="<<totMass<<
" > M="<<totM<<
",S="<<totS
2622 <<
", envPDG="<<envPDG<<
G4endl;
2627 G4cout<<
"G4Q::HQ:f="<<fprob<<
",d="<<dm<<
",rPDG="<<rPDG<<
",rM="<<rMass<<
",M="<<reMass
2635 FillHadronVector(quasH);
2637 qEnv=theEnvironment;
2640 else G4cerr<<
"---Warning---G4Q::HQ:Q=H,q="<<iniPDG<<
",s="<<sPDG<<
",d="<<dm<<
G4endl;
2644 if(rPDG!=10&&rMass>dm&&!rWi)
2652 if(fabs(sMM-ddm)<1.5*sWi-.001 && ddm>mmm)
2657 sMass=GetRandomMass(sPDG,ddm);
2658 if(fabs(sMass)<.001)
2661 G4cerr<<
"***G4Q::HQ:ChangeToM=0, "<<sPDG<<
",new="<<ddm<<
",old="<<msm<<
G4endl;
2665 if(sMass<ddm) sMass=ddm;
2667 G4cout<<
"G4Q::HQ: sPDG="<<sPDG<<
",sM="<<sMass<<
",d="<<ddm<<
",isM="<<msm<<
",W="
2687 G4cout<<
"G4Q::HQ:BEFrPDGcor,d="<<dm<<
",R="<<rnd<<
",r="<<rPDG<<
",rM="<<rMass<<
G4endl;
2691 if(rPDG==111 && sPDG!=111 && dm>548.)
2694 if(dm>958.) rPDG=331;
2697 if(rPDG==221 && dm>958. && rnd>.5 ) rPDG=331;
2698 if(rPDG==331 &&(dm<958. || rnd<.5)) rPDG=221;
2700 if(rPDG==221 && dm<548.) rPDG=111;
2703 if ( ( (rPDG == 111 && sPDG!= 111) || rPDG == 221) &&
2704 rMass > 544. && dm > 544. && rnd > .5) rPDG=113;
2706 if ( ( (rPDG == 111 && sPDG != 111) || rPDG == 221) &&
2707 rMass > 782. && dm > 782. && rnd < .5) rPDG = 223;
2709 if ( rPDG == 331 && rMass > 1020. && dm > 1020. && rnd < .5) rPDG=333;
2711 if(rPDG== 211 && dm>544. && rnd>.5) rPDG= 213;
2712 if(rPDG==-211 && dm>544. && rnd>.5) rPDG=-213;
2714 G4cout<<
"G4Q::HQ:rCor,Q="<<quasM<<
",sM="<<sMass<<
",r="<<rPDG<<
",rM="<<rMass<<
G4endl;
2716 if (rPDG < MINPDG && rPDG != 2212 && rPDG != 2112 && rPDG != 3122 && rPDG != 10)
2718 reMass=GetRandomMass(rPDG,dm);
2720 G4cout<<
"G4Q::HQ:dm="<<dm<<
", ResQM="<<reMass<<
" is changed to PDG="<<rPDG<<
G4endl;
2724 if(sPDG==221 || sPDG==331)
2726 if (sPDG==221) dm+=mEta-mPi0;
2727 else if(sPDG==331) dm+=mEtaP-mPi0;
2739 if(dm<mPi0-.00001&&rPDG==111)
2744 else reMass=GetRandomMass(rPDG,dm);
2745 if(reMass==0.)
G4cerr<<
"-W-G4Q::HQ:2,M="<<quasM<<
",r="<<rPDG<<
",d="<<dm<<
G4endl;
2758 ed <<
"Why Fail? (7): NeedsEvap7:s=" << sPDG <<
",Q=" << q4Mom
2759 << valQ <<
",r=" << rPDG <<
G4endl;
2762 qEnv=theEnvironment;
2767 else if(rPDG==NUCPDG)
2783 if(fRQW<.001) fRQW=.001;
2785 G4int sChg=sQPDG.GetCharge();
2786 G4int sBaryn=sQPDG.GetBaryNum();
2789 G4cout<<
"G4Q::HQ:h="<<sCB<<
",C="<<sChg<<
",B="<<sBaryn<<
",E="<<theEnvironment<<
G4endl;
2795 G4cout<<
"G4Q::HQ:rqCB="<<rCB<<
",rqC="<<rChg<<
",rqB="<<sBaryn<<
",rM="<<rQPDG<<
",reM="
2798 if ( totBN > 1 && totS >= 0 && envPDG > MINPDG && envPDG != NUCPDG &&
2799 (reMass+sMass > quasM || sCB+rCB+reMass+sMass+envM > totMass ||
2800 (!RQB && quasM < diPiM)
2807 ed <<
"Why Fail? (8): RQM+SM=" << reMass+sMass <<
">QM=" << quasM <<
", sCB="
2808 << sCB <<
" + rCB=" << rCB <<
" + rM=" << reMass <<
" + sMass=" << sMass
2809 <<
" + eM=" << envM <<
" = " << sCB+rCB+reMass+sMass+envM <<
">tM=" << totMass
2810 <<
"," << reMass+sMass+envM <<
G4endl;
2813 qEnv=theEnvironment;
2819 ed <<
"Residual Particle is Vacuum: rPDG=90000000, MV=" << reMass <<
G4endl;
2822 if(rPDG==2212&&sPDG==311&&reMass+sMass>quasM)
2824 if(mNeut+mK<=quasM+.001)
2833 G4cout<<
"G4Q::HQ:NCB="<<rCB<<
",NC="<<rChg<<
",sB="<<sBaryn<<
",r="<<rQPDG<<
G4endl;
2839 if(mNeut+mK<=quasM) sMass=quasM-mNeut;
2842 sChg=sQPDG.GetCharge();
2843 sBaryn=sQPDG.GetBaryNum();
2846 G4cout<<
"G4Q::HQ:KCB="<<sCB<<
",KC="<<sChg<<
",frB="<<sBaryn<<
",E="<<theEnvironment
2854 ed<<
"Can't decay Q in N and K: (NK) QM="<< quasM<<
",d="<< quasM-mNeut-mK<<
G4endl;
2859 G4cout<<
"G4Q::HQ: ****** Before reM="<<reMass<<
", rM="<<rMass<<
G4endl;
2862 if(tmpQPDG.GetWidth()<.000001) reMass=tmpQPDG.GetMass();
2863 if(!reMass) reMass=rMass;
2865 G4cout<<
"G4Q::HQ: Decay in sM="<<sMass<<
" + reM="<<reMass<<
" (rM="<<rMass<<
G4endl;
2872 G4cout<<
"G4Q::HQ:Q->RQ+QEX s="<<sPDG<<
",pM="<<pMass<<
",E="<<theEnvironment<<
G4endl;
2878 if(fabs(tmM-sum)<eps)
2880 r4Mom=q4Mom*(reMass/sum);
2881 s4Mom=q4Mom*(sMass/sum);
2887 G4cerr<<
"---Warning---G4Q::HQ:M="<<tmM<<
"=>rPDG="<<rPDG<<
"(rM="<<reMass<<
")+sPDG="
2888 <<sPDG<<
"(sM="<<sMass<<
")="<<sum<<
",resNQC="<<resNQC<<
G4endl;
2892 if(sPDG==311 && tmpQPDG.GetCharge()>0)
2894 G4QContent crQC=tmpQPDG.GetQuarkContent()-KpQC+K0QC;
2907 if(fabs(tmM-sum)<eps)
2909 r4Mom=q4Mom*(reMass/sum);
2910 s4Mom=q4Mom*(sMass/sum);
2915 ed <<
"Hadron+K+ DecayIn2: (I) KCor M=" << tmM <<
"=>rPDG=" << rPDG <<
"(rM="
2916 << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum <<
G4endl;
2917 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0029",
2924 ed <<
"Hadron+K+ DecayIn2:(O) KCor M=" << tmM <<
"=>rPDG=" << rPDG <<
"(rM="
2925 << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum <<
G4endl;
2929 else if(sPDG==321 && tmpQPDG.GetCharge()<=tmpQPDG.GetBaryNum())
2931 G4QContent crQC=tmpQPDG.GetQuarkContent()-K0QC+KpQC;
2944 if(fabs(tmM-sum)<eps)
2946 r4Mom=q4Mom*(reMass/sum);
2947 s4Mom=q4Mom*(sMass/sum);
2952 ed <<
"Hadron+K0 DecayIn2: (I) K0Cor M=" << tmM <<
"=>rPDG=" << rPDG
2953 <<
"(rM=" << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")="
2955 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0031",
2962 ed <<
"Hadron+K0 DecayIn2: (O) K0Cor M=" << tmM <<
"=>rPDG=" << rPDG <<
"(rM="
2963 << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum <<
G4endl;
2967 else if(sPDG==211 && tmpQPDG.GetCharge()<tmpQPDG.GetBaryNum())
2969 G4QContent crQC=tmpQPDG.GetQuarkContent()-Pi0QC+PiQC;
2982 if(fabs(tmM-sum)<eps)
2984 r4Mom=q4Mom*(reMass/sum);
2985 s4Mom=q4Mom*(sMass/sum);
2990 ed <<
"Hadron+Pi0 DecayIn2: (I) Pi+/Pi0Cor M=" << tmM <<
"=>rPDG=" << rPDG
2991 <<
"(rM=" << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum
2993 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0033",
3007 if(fabs(tmM-sum)<eps)
3009 r4Mom=q4Mom*(reMass/sum);
3010 s4Mom=q4Mom*(sMass/sum);
3015 ed <<
"Hadron+Gamma DecayIn2: (I) Pi+/GamCor M=" << tmM <<
"=>rPDG=" << rPDG
3016 <<
"(rM=" << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum
3018 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0034",
3025 ed <<
"Hadron+Pi0/Gam DecayIn2: (O) Pi+/Pi0Cor M=" << tmM <<
"=>rPDG=" << rPDG
3026 <<
"(rM="<< reMass <<
")+sPDG="<< sPDG <<
"(sM="<< sMass <<
")="<< sum <<
G4endl;
3030 else if(sPDG==-211 && tmpQPDG.GetCharge()>0)
3032 G4QContent crQC=tmpQPDG.GetQuarkContent()-Pi0QC+PiMQC;
3045 if(fabs(tmM-sum)<eps)
3047 r4Mom=q4Mom*(reMass/sum);
3048 s4Mom=q4Mom*(sMass/sum);
3053 ed <<
"Hadron+Pi0 DecayIn2: (I) Pi-/Pi0Cor M=" << tmM <<
"=>rPDG=" << rPDG
3054 <<
"(rM=" << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum
3056 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0036",
3070 if(fabs(tmM-sum)<eps)
3072 r4Mom=q4Mom*(reMass/sum);
3073 s4Mom=q4Mom*(sMass/sum);
3078 ed <<
"Hadron+Gamma DecayIn2: (I) Pi-/GamCor M=" << tmM <<
"=>rPDG=" << rPDG
3079 <<
"(rM=" << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum
3081 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0037",
3086 else if((sPDG==221 || sPDG==331) && tmM>mPi0+reMass)
3092 if(fabs(tmM-sum)<eps)
3094 r4Mom=q4Mom*(reMass/sum);
3095 s4Mom=q4Mom*(sMass/sum);
3100 ed <<
"Hadron+Pi0 DecayIn2: Eta/Pi0Cor M="<< tmM <<
"=>rPDG="<< rPDG <<
"(rM="
3101 << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum <<
G4endl;
3105 else if((sPDG==111 || sPDG==221 || sPDG==331) && tmM>reMass)
3111 if(fabs(tmM-reMass)<eps)
3113 r4Mom=q4Mom*(reMass/sum);
3114 s4Mom=q4Mom*(sMass/sum);
3119 ed <<
"QHadron+Gamma DecayIn2: PiCor M=" << tmM <<
"=>rPDG=" << rPDG <<
"(rM="
3120 << reMass <<
")+sPDG="<< sPDG <<
"(sM=" << sMass <<
")=" << reMass <<
G4endl;
3124 else if(iniBN>0 && iniS>0)
3127 G4QContent lanQC=tmpQPDG.GetQuarkContent()+tmpSQC+K0QC;
3133 G4cout<<
"G4Q::HQ:LsPDG="<<sPDG<<
",rPDG="<<rPDG<<
",Z="<<nucZ<<
",M="<<nucM<<
G4endl;
3135 if((
G4UniformRand()<.33333 || mPi+nreM>tmM) && mPi0+nreZ<tmM)
3145 if(fabs(tmM-sum)<eps)
3147 r4Mom=q4Mom*(reMass/sum);
3148 s4Mom=q4Mom*(sMass/sum);
3153 ed <<
"(L->n)+Pi0 DecayIn2: LamPi0 Cor M="<< tmM <<
"=>rPDG="<< rPDG <<
"(rM="
3154 << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum <<
G4endl;
3155 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0040",
3159 else if(mPi+nreM<tmM)
3163 curQ+=tmpSQC+K0QC-PiMQC;
3169 if(fabs(tmM-sum)<eps)
3171 r4Mom=q4Mom*(reMass/sum);
3172 s4Mom=q4Mom*(sMass/sum);
3177 ed <<
"(L->n)+Pi- DecayIn2: LamPiM Cor M=" << tmM <<
"=>rPDG=" << rPDG
3178 <<
"(rM=" << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum
3180 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0041",
3194 if(fabs(tmM-sum)<eps)
3196 r4Mom=q4Mom*(reMass/sum);
3197 s4Mom=q4Mom*(sMass/sum);
3202 ed <<
"(L->n)+Gamma DecayIn2: LamNGam Cor M=" << tmM <<
"=>rPDG=" << rPDG
3203 <<
"(rM=" << reMass <<
")+sPDG=" << sPDG <<
"(sM=" << sMass <<
")=" << sum
3205 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0042",
3212 ed <<
"LamTo0N with Pi DecayIn2: LamToN M=" << tmM << totQC <<
"=>rM="
3213 << nucM.GetPDG() <<
"," << nucZ.GetPDG() <<
"(" << nreM <<
"," << nreZ
3214 <<
")+PiM/PiZ=" << mPi+nreM <<
"," << mPi0+nreZ <<
G4endl;
3215 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0043",
3230 if(fabs(tmM-reMass)<eps)
3232 r4Mom=q4Mom*(reMass/sum);
3233 s4Mom=q4Mom*(sMass/sum);
3238 ed <<
"QHadron+Gamma DecayIn2: gam+TQ M=" << tmM <<
"=>rPDG=" << rPDG <<
"(rM="
3239 << reMass <<
")+sPDG="<< sPDG <<
"(sM=" << sMass <<
")=" << reMass <<
G4endl;
3240 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0044",
3244 else if(totMass>resTNM+sMass)
3249 G4cout<<
"G4Q::HQ:EMERGENCY,rEM="<<resTN<<resTNM<<
",fM="<<sMass<<
",tM="<<totMass
3250 <<
",d="<<totMass-resTNM-sMass<<
G4endl;
3253 if(fabs(totMass-sum)<eps)
3255 re4M=tot4M*(resTNM/sum);
3256 rs4M=tot4M*(sMass/sum);
3261 ed <<
"DecayIn2 Frag+ResE failed: HadrQ:Decay T=" << totMass
3262 <<
"->R=" << resTNM <<
"+S=" << sMass <<
")=" << sum <<
G4endl;
3269 FillHadronVector(fragH);
3272 resTN.Set4Momentum(re4M);
3278 FillHadronVector(envH);
3285 else if(totMass>totM)
3290 G4cout<<
"G4Q::HQ:EMERGENSY,minM="<<totM<<
" < totM="<<totMass<<
G4endl;
3292 if(fabs(totMass-totM)<eps) re4M=tot4M*(resTNM/sum);
3296 ed<<
"DecayIn2 gam+TotN failed:HadrQ:Decay,T="<<totMass<<
"->g+M="<<totM<<
G4endl;
3297 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0046",
3304 FillHadronVector(fragH);
3307 totN.Set4Momentum(re4M);
3313 FillHadronVector(envH);
3322 G4cerr<<
"***G4Q::HQ:M="<<tmM<<
"=>rPDG="<<rPDG<<
"(rM="<<reMass<<
")+sPDG="
3323 <<sPDG<<
"(sM="<<sMass<<
")="<<sum<<
",QM="<<iniQM<<
G4endl;
3324 if(fabs(tmM-sum)<1.)
3327 r4Mom=q4Mom*(reMass/sum);
3328 s4Mom=q4Mom*(sMass/sum);
3331 else G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0047",
3337 G4cout<<
"G4Q::HQ:=2.3=>QHVect s4M="<<s4Mom<<
",sPDG="<<sPDG<<
", r4M/M="<<r4Mom<<reMass
3338 <<
",fR="<<freeRQM<<
",fW="<<fRQW<<
",PDG="<<rPDG<<
",r="<<rCB<<
",s="<<sCB<<
G4endl;
3345 G4cout<<
"****G4Q::HQ:E-9: sKE="<<sKE<<
"<sCB="<<sCB<<
G4endl;
3347 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0048",
3351 qEnv=theEnvironment;
3354 else if(abs(reMass-freeRQM)<fRQW||envPDG==NUCPDG)
3357 FillHadronVector(curHadr1);
3359 FillHadronVector(curHadr2);
3361 G4cout<<
"G4Q::HQ:DecayQuasmon "<<q4Mom<<
" in 4M="<<r4Mom+s4Mom<<
" RQ="<<rPDG<<r4Mom
3362 <<
" + Fragment="<<sPDG<<s4Mom<<
", Env="<<theEnvironment<<
G4endl;
3364 if(sPDG>MINPDG) theEnvironment.
Reduce(pPDG);
3366 qEnv=theEnvironment;
3375 G4double resTotNM=resTotN.GetMZNS();
3376 if(resTotN4Mom.
m()<resTotNM)
3381 ed<<
"Why Fail?(10):NEEDS-EVAP-10,M="<<resTotN4Mom.
m()<<
"<miM="<<resTotNM<<
G4endl;
3385 qEnv=theEnvironment;
3391 FillHadronVector(curHadr2);
3395 theEnvironment.
Reduce(pPDG);
3400 G4cout<<
"OK***>G4Q::HQ:S="<<sPDG<<s4Mom<<
",Env="<<theEnvironment<<
",Q="<<q4Mom
3410 G4cout<<
"***>G4Q::HQ:After,S="<<sPDG<<s4Mom<<
",Env="<<theEnvironment<<
",Q="
3411 <<q4Mom<<valQ<<curQ<<
G4endl;
3413 qEnv=theEnvironment;
3419 else G4cout<<
"G4Q::HQ:NO-OK,h="<<hsflag<<
",d="<<fdul<<
",M="<<rMass<<
"<"<<reMass<<
",M2="
3420 <<reTNM2<<
"<I="<<tmpTM2<<
",sP="<<sPDG<<
",eP="<<envPDG<<
",pP="<<pPDG<<
G4endl;
3428 G4cout<<
"G4Q::HQ:>>"<<sPDG<<fr4Mom<<fr4Mom.m()<<
"="<<sMass<<
",T="<<frKin<<
",E="<<ePDG
3437 G4cout<<
"G4Q::HQ: sPDG="<<sPDG<<
", sM="<<sMass<<
", SQ="<<SQ<<
G4endl;
3439 if(!sPDG&&SQ<0&&nQuasms==1)
3447 G4int rKPPDG = rKPN.GetPDG();
3450 G4int rK0PDG = rK0N.GetPDG();
3452 if ( (rKPM+mK > totMass && rK0M+mK0 > totMass) ||
3458 ed <<
"Why PANIC? (2): ***PANIC#2***tM=" << totMass <<
"<KM=" << mK <<
","
3459 << mK0 <<
",rM=" << rKPM <<
"," << rK0M <<
",d=" << mK+rKPM-totMass <<
","
3460 << mK0+rK0M-totMass <<
G4endl;
3464 qEnv=theEnvironment;
3467 if(rKPM + mK > rK0M + mK0)
3485 if(fabs(ctM-sum)<eps)
3487 r4Mom=tot4M*(rMass/sum);
3488 s4Mom=tot4M*(sMass/sum);
3493 ed <<
"HadrQuasm:K+ResNuc DecayIn2 didn't succeed: tM=" << ctM
3494 << totQC <<
" => rPDG=" << rPDG <<
"(rM=" << rMass <<
") + sPDG="
3495 << sPDG <<
"(sM=" << sMass <<
")=" << sum <<
G4endl;
3499 G4cout<<
"G4Q::HQ:===2.4===>HadrVec s="<<sPDG<<s4Mom<<
",r="<<rPDG<<r4Mom<<
G4endl;
3503 FillHadronVector(curHadr1);
3505 FillHadronVector(curHadr2);
3511 if(quasM<rMass+sMass&&(sPDG==221||sPDG==331))
3518 if (iniS<0&&iniQChg+iniQChg>=iniBN)
3523 rPDG = totQN.GetPDG();
3524 rMass= totQN.GetMZNS();
3531 rPDG = totQN.GetPDG();
3532 rMass= totQN.GetMZNS();
3534 else if(iniQChg>iniBN-iniS)
3539 rPDG = totQN.GetPDG();
3540 rMass= totQN.GetMZNS();
3547 rPDG = totQN.GetPDG();
3548 rMass= totQN.GetMZNS();
3550 else if(quasM>iniQM+mPi0)
3567 G4cout<<
"G4Q::HQ:MQ="<<q4Mom.
m()<<
"->sPDG="<<sPDG<<
"(M="<<sMass<<
") + rPDG="<<rPDG
3568 <<
"(M="<<rMass<<
")"<<
",S="<<rMass+sMass<<
G4endl;
3570 if(q4Mom.
m()+.003<rMass+sMass)
3573 G4cerr<<
"G4Q::HQ:***PANIC#3***tM="<<q4Mom.
m()<<
"<rM="<<rMass<<
",sM="<<sMass
3574 <<
",d="<<rMass+sMass-q4Mom.
m()<<
G4endl;
3578 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0052",
3582 qEnv=theEnvironment;
3589 if(fabs(cqM-sum)<eps)
3591 resQ4Mom=q4Mom*(rMass/sum);
3592 s4Mom=q4Mom*(sMass/sum);
3600 ed <<
"Quasm+Hadr DecayIn2 error: MQ=" << cqM <<
"-> rPDG=" << rPDG
3601 <<
", (M=" << rMass <<
") + sPDG=" << sPDG <<
"(M=" << sMass
3602 <<
")=" << sum <<
G4endl;
3603 G4Exception(
"G4Quasmon::HadronizeQuasmon()",
"HAD_CHPS_0053",
3607 G4cout<<
"G4Q::HQ:Decay of Quasmon="<<q4Mom<<
"->s="<<sPDG<<s4Mom<<
"+R="<<resQ4Mom
3614 theQHadrons.push_back(candHadr);
3616 FillHadronVector(candHRes);
3618 qEnv=theEnvironment;
3626 if(theEnvironment.
GetPDG()>NUCPDG)
3631 G4double outT = s4Mom.e()-s4Mom.m();
3632 outProb = theEnvironment.
CoulBarPenProb(outCB,outT,outChg,outBar);
3636 G4cout<<
"G4Q::HQ: for "<<sPDG<<
", rnd="<<rnd<<
" < outP="<<outProb<<
" ?"<<
G4endl;
3640 FillHadronVector(candHadr);
3668 G4double outT = fr4Mom.e()-fr4Mom.m();
3672 theEnvironment.
Reduce(pPDG);
3695 if(sPDG>MINPDG) valQ += transQC;
3697 FillHadronVector(candHadr);
3699 G4cout<<
"G4Q::HQ:QuarkExchHadronizThroughCB Q="<<valQ<<
",trQC="<<transQC<<
G4endl;
3703 G4cout<<
"G4Q::HQ:status=1, ------>> NuclearMatter SUBCHECK ---->>"<<sumL<<
G4endl;
3713 G4cout<<
"G4Q::HQ:CBIsn'tPEN,P="<<outProb<<
",T="<<outT<<
",M="<<fr4Mom.m()<<
G4endl;
3729 if(CheckGroundState())
3733 qEnv=theEnvironment;
3739 G4int sumC = eZ+valQ.GetCharge()+ccheck;
3740 G4int curPDG=valQ.GetSPDGCode();
3741 G4cout<<
"G4Q::HQ:Z="<<eZ<<valQ<<
"***>FinalCHECK***>>4M="<<sumLor<<
",Ch="<<sumC<<
G4endl;
3742 if(!curPDG)
G4cout<<
"***G4Q::HQ: Quasmon-Tripolino QC="<<valQ<<
G4endl;
3743 G4cout<<
"G4Q::HQ:=------=> ResidualQ 4M="<<q4Mom<<
", QC="<<valQ<<
G4endl;
3747 G4int ecSum=theEnvironment.
GetZ()+valQ.GetCharge();
3748 G4int nHe=theQHadrons.size();
3749 if(nHe)
for(
int ih=0; ih<nHe; ih++) ecSum+=theQHadrons[ih]->
GetCharge();
3752 G4cerr<<
"***G4Q::HQ:C"<<cSum<<
",c="<<ecSum<<
",E="<<theEnvironment<<
",Q="<<valQ<<
G4endl;
3753 G4cerr<<
":G4Q::HQ:*END*,oE="<<oldEnv<<
"oQ="<<oldCQC<<
",oN="<<oldNH<<
",N="<<nHe<<
G4endl;
3754 if(nHe)
for(
G4int h=0; h<nHe; h++)
3762 G4cout<<
"G4Q::HQ: Q="<<q4Mom<<valQ<<
",E="<<theEnvironment<<
", status="<<status<<
G4endl;
3764 qEnv=theEnvironment;
3772 FillHadronVector(thisQuasmon);
3774 G4int nHadrs=theQHadrons.size();
3776 G4cout<<
"G4Q::DecayQuasmon:After decay (FillHadronVector byItself) nH="<<nHadrs<<
G4endl;
3778 if(nHadrs)
for (
int hadron=0; hadron<nHadrs; hadron++)
3781 theFragments->push_back(curHadr);
3784 else G4cerr<<
"*******G4Quasmon::DecayQuasmon: *** Nothing is in the output ***"<<
G4endl;
3788 return theFragments;
3792 void G4Quasmon::FillHadronVector(
G4QHadron* qH)
3831 if(thePDG==113 && fabs(t.
m()-770.)<.001)
3833 G4cerr<<
"G4Q::FillHadronVector: PDG="<<thePDG<<
",M="<<t.
m()<<
G4endl;
3835 G4Exception(
"G4Quasmon::FillHadronVector()",
"HAD_CHPS_0000",
3840 G4cout<<
"G4Q::FillHadronVector:Hadron's PDG="<<thePDG<<
",4Mom="<<t<<
",m="<<t.
m()<<
G4endl;
3842 if(thePDG>80000000 && (thePDG<90000000 || thePDG%1000>500 || thePDG%1000000>500000)
3843 && thePDG!=90002999 && thePDG!=89999003 && thePDG!=90003998 && thePDG!=89998004
3844 && thePDG!=90003999 && thePDG!=89999004 && thePDG!=90004998 && thePDG!=89998005)
3846 if (thePDG==90999999) thePDG=-311;
3847 else if(thePDG==90999000) thePDG=-321;
3848 else if(thePDG==89000001) thePDG=311;
3849 else if(thePDG==89001000) thePDG=321;
3850 else if(thePDG==90000999) thePDG=211;
3851 else if(thePDG==89999001) thePDG=-211;
3852 else if(thePDG==89999999) thePDG=-2112;
3853 else if(thePDG==89999000) thePDG=-2212;
3854 else if(thePDG==89000000) thePDG=-3122;
3855 else if(thePDG==88999002) thePDG=-3222;
3856 else if(thePDG==89000999) thePDG=-3222;
3857 else if(thePDG==88000001) thePDG=-3322;
3858 else if(thePDG==88001000) thePDG=-3312;
3859 else if(thePDG==89999002) thePDG=1114;
3860 else if(thePDG==90001999) thePDG=2224;
3861 else if(thePDG==91000000) thePDG=3122;
3862 else if(thePDG==90999001) thePDG=3112;
3863 else if(thePDG==91000999) thePDG=3222;
3864 else if(thePDG==91999000) thePDG=3312;
3865 else if(thePDG==91999999) thePDG=3322;
3866 else if(thePDG==92998999) thePDG=3112;
3880 if(!h1PDG || !h2PDG)
3882 G4cerr<<
"***FillHV:h1QC="<<h1QC<<
"(PDG="<<h1PDG<<
"),h2QC="<<h2QC<<
"(PDG="<<h2PDG<<
")"
3885 "Chipolino can't be defragmented");
3895 G4cerr<<
"***G4Q::FillHadrV:ChipQC"<<chipQC<<
":PDG1="<<h1PDG<<
",PDG2="<<h2PDG<<
G4endl;
3896 theQHadrons.push_back(qH);
3902 FillHadronVector(fHadr);
3904 FillHadronVector(sHadr);
3907 else if(thePDG>80000000&&thePDG!=90000000)
3917 G4int nZ =qNuc.GetZ();
3918 G4int nS =qNuc.GetS();
3919 G4int bA =qNuc.GetA();
3921 G4cout<<
"G4Quasm::FillHadrVect:Nucl="<<qNuc<<
",nPDG="<<thePDG<<
",GSM="<<GSMass<<
G4endl;
3923 if((nN<0 || nZ<0 || nS<0) && bA>0)
3928 G4int newS=newNpm.GetStrangeness();
3929 if(newS>0) newNpm=
G4QNucleus(totQC+PiQC+newS*K0QC);
3930 G4int PDG2=newNpm.GetPDG();
3931 G4double m2_value=newNpm.GetMZNS();
3937 PDG2 =newNp.GetPDG();
3938 m2_value =newNp.GetMZNS();
3941 if (m3_value+mK0<m2_value+mK)
3946 PDG2=newN0.GetPDG();
3949 else if(nS>0&&nZ+nN>0)
3956 PDG2 =newNp.GetPDG();
3957 m2_value =newNp.GetMZNS();
3964 PDG2 =newNp.GetPDG();
3965 m2_value =newNp.GetMZNS();
3972 PDG2 =newNpp.GetPDG();
3973 m2_value =newNpp.GetMZNS();
3975 if(fragMas>m1+m2_value)
3982 ed <<
"Mes+ResA DecayIn2 did not succeed: QM=" << t.
m() <<
"-> Mes=" << PDG1
3983 <<
"(M=" << m1 <<
") + ResA=" << PDG2 <<
"(M=" << m2_value <<
")" <<
G4endl;
3988 theQHadrons.push_back(H1);
3990 FillHadronVector(H2);
3992 else if(fabs(m1+m2_value-fragMas)<0.01)
4002 theQHadrons.push_back(H1);
4004 FillHadronVector(H2);
4009 G4cerr<<
"-Warning-G4Q::FillHVec:PDG="<<thePDG<<
"("<<t.
m()<<
","<<fragMas<<
") < Mes="
4010 <<PDG1<<
"("<<m1<<
") + ResA="<<PDG2<<
"("<<m2_value<<
"), d="<<fragMas-m1-m2_value<<
G4endl;
4013 theQHadrons.push_back(qH);
4016 else if(abs(fragMas-GSMass)<.1)
4027 nResPDG=resN.GetPDG();
4028 if (nResPDG==90000001) nResM=mNeut;
4029 else if(nResPDG==90001000) nResM=mProt;
4030 else if(nResPDG==91000000) nResM=mLamb;
4031 else nResM=resN.GetMZNS();
4039 pResPDG=resN.GetPDG();
4040 if (pResPDG==90000001) pResM=mNeut;
4041 else if(pResPDG==90001000) pResM=mProt;
4042 else if(pResPDG==91000000) pResM=mLamb;
4043 else pResM =resN.GetMZNS();
4051 lResPDG=resN.GetPDG();
4052 if (lResPDG==90000001) lResM=mNeut;
4053 else if(lResPDG==90001000) lResM=mProt;
4054 else if(lResPDG==91000000) lResM=mLamb;
4055 else lResM =resN.GetMZNS();
4058 G4cout<<
"G4Quasm::FillHadrVec:rP="<<pResPDG<<
",rN="<<nResPDG<<
",rL="<<lResPDG<<
",nN="
4059 <<nN<<
",nZ="<<nZ<<
",nL="<<nS<<
",totM="<<fragMas<<
",n="<<fragMas-nResM-mNeut
4060 <<
",p="<<fragMas-pResM-mProt<<
",l="<<fragMas-lResM-mLamb<<
G4endl;
4062 if ( thePDG == 90004004 ||
4063 (bA > 1 && ( (nN > 0 && fragMas > nResM+mNeut) ||
4064 (nZ > 0 && fragMas > pResM+mProt) ||
4065 (nS > 0 && fragMas > lResM+mLamb) ) ) )
4067 G4int barPDG = 90002002;
4068 G4int resPDG = 90002002;
4072 if (fragMas > nResM+mNeut) {
4078 else if(fragMas>pResM+mProt)
4085 else if(fragMas>lResM+mLamb)
4092 else if(thePDG!=90004004 && fragMas>GSMass)
4099 else if(thePDG!=90004004)
4102 ed <<
"Below GSM but cann't decay: PDG=" << thePDG <<
",M=" << fragMas
4103 <<
"<GSM=" << GSMass <<
G4endl;
4110 theQHadrons.push_back(qH);
4111 G4cerr<<
"---Warning---G4Q::FillHadronVector: Be8 decay did not succeed"<<
G4endl;
4121 FillHadronVector(HadrB);
4123 FillHadronVector(HadrR);
4129 G4cout<<
"G4Quasm::FillHadrVect: Leave as it is"<<
G4endl;
4131 theQHadrons.push_back(qH);
4134 else if (fragMas < GSMass)
4136 G4cerr<<
"***G4Quasmon::FillHV:M="<<fragMas<<
">GSM="<<GSMass<<
"(PDG="<<thePDG<<
"),d="
4137 <<fragMas-GSMass<<
", NZS="<<nN<<
","<<nZ<<
","<<nS<<
G4endl;
4139 G4cout<<
"****>>G4Quasm::FillHadrVect: Leave as it is Instead of Exception"<<
G4endl;
4140 theQHadrons.push_back(qH);
4142 else if (bA==1 && fragMas>GSMass)
4146 if(fragMas>mPi0+GSMass)
4155 theQHadrons.push_back(qH);
4156 G4cerr<<
"---Warning---G4Q::FillHadrVect:N*->gamma/pi0+N decay error"<<
G4endl;
4161 FillHadronVector(HadrB);
4163 theQHadrons.push_back(qH);
4169 G4cout<<
"G4Quasm::FillHadrVect:Evaporate "<<thePDG<<
",tM="<<fragMas<<
" > GS="<<GSMass
4170 <<qNuc.Get4Momentum()<<
", m="<<qNuc.Get4Momentum().m()<<
G4endl;
4174 if(!qNuc.EvaporateBaryon(bHadron,rHadron))
4176 G4cerr<<
"---Warning---G4Q::FillHV:Evaporate PDG="<<thePDG<<
",M="<<fragMas<<
G4endl;
4179 theQHadrons.push_back(qH);
4190 G4int tmpS=tmpQHadVec->size();
4191 theQHadrons.resize(tmpS+theQHadrons.size());
4192 copy( tmpQHadVec->begin(), tmpQHadVec->end(), theQHadrons.end()-tmpS);
4193 tmpQHadVec->clear();
4196 else FillHadronVector(bHadron);
4200 G4int tmpS=tmpQHadVec->size();
4201 theQHadrons.resize(tmpS+theQHadrons.size());
4202 copy( tmpQHadVec->begin(), tmpQHadVec->end(), theQHadrons.end()-tmpS);
4203 tmpQHadVec->clear();
4206 else FillHadronVector(rHadron);
4217 G4cout<<
"G4Q::FillHV: ---DECAY IS DONE--- with nH="<<tmpQHadVec->size()<<
G4endl;
4219 G4int tmpS=tmpQHadVec->size();
4220 theQHadrons.resize(tmpS+theQHadrons.size());
4221 copy( tmpQHadVec->begin(), tmpQHadVec->end(), theQHadrons.end()-tmpS);
4223 G4cout<<
"G4Q::FillHV: -->Products are added to QHV, nQHV="<<theQHadrons.size()<<
G4endl;
4225 tmpQHadVec->clear();
4228 G4cout<<
"G4Q::FillHV: TemporaryQHV of DecayProducts is deleted"<<
G4endl;
4241 G4cout<<
"G4Quas::GetQPartonMomentum:***called*** kMax="<<kMax<<
",mC="<<sqrt(mC2)<<
G4endl;
4259 if (kLim<kMax) kMax = kLim;
4260 if (kMin<0 || kMax<0 || qMass<=0. || nOfQ<2)
4263 ed <<
"Can not generate quark-parton: kMax=" << kMax <<
", kMin=" << kMin
4264 <<
", kLim=" << kLim <<
", MQ=" << qMass <<
", n=" << nOfQ <<
G4endl;
4268 G4cout<<
"G4Q::GetQPM: kLim="<<kLim<<
",kMin="<<kMin<<
",kMax="<<kMax<<
",nQ="<<nOfQ<<
G4endl;
4270 if(kMin>kMax||nOfQ==2)
return kMax;
4280 if (kMax>=kLim) vRndm = vRndm*pow((1.-xMin),n)*(1.+n*xMin);
4284 G4double vRmin = pow((1.-xMin),n)*(1.+n*xMin);
4285 G4double vRmax = pow((1.-xMax),n)*(1.+n*xMax);
4286 vRndm = vRmax + vRndm*(vRmin-vRmax);
4292 G4double vRmax = pow((1.-xMax),n)*(1.+n*xMax);
4293 vRndm = vRmax + vRndm*(1.-vRmax);
4295 if (vRndm<=0. || vRndm>1.)
4299 if(vRndm<=0.) vRndm=1.e-9;
4300 else if(vRndm>1.) vRndm=1.;
4302 if (n==1)
return kLim*sqrt(1.-vRndm);
4305 G4double x = 1.-pow(vRndm*(1+n*vRndm)/(fn+1.),1./fn);
4309 G4double df = 1./
static_cast<double>(nOfQ);
4310 G4double f = df*(
static_cast<int>(nOfQ*nOfQ*n*x/5.)+(nOfQ/2));
4321 G4cout<<
"G4Q::GetQPMom: f="<<
f<<
" is changed to 99"<<
G4endl;
4325 if(x<1.
e-27) x=1.e-27;
4326 else if(x>.999999999) x=.999999999;
4329 if (n>2) p = pow(r,n-1);
4334 G4cout<<
"G4Q::GetQPMom:--->> First x="<<x<<
", n="<<n<<
", f="<<
f<<
", d/R(first)="
4338 if(nitMax>100)nitMax=100;
4339 while( abs(
d/vRndm) > 0.001 && it <= nitMax)
4341 x = x +
f*od/(r*nx*(fn+1.));
4342 if(x<1.
e-27) x=1.e-27;
4343 else if(x>.999999999) x=.999999999;
4345 if (n>2) p = pow(r,n-1);
4350 if ((od>0&&d<0)||(od<0&&d>0))
4352 if (
f>1.0001)
f=1.+(
f-1.)/2.;
4353 if (
f<0.9999)
f=1.+(1.-
f)*2;
4355 if(x<1.
e-27) x=1.e-27;
4356 else if(x>.999999999) x=.999999999;
4358 if (n>2) p = pow(r,n-1);
4366 if (
f>1.0001&&
f<27.)
f=1.+(
f-1.)*2;
4367 if (
f<0.99999999999)
f=1.+(1.-
f)/2.;
4371 G4cout<<
"G4Q::GetQPMom: Iter#"<<it<<
": (c="<<c<<
" - R="<<vRndm<<
")/R ="<<d/vRndm
4372 <<
", x="<<x<<
", f="<<
f<<
G4endl;
4376 if(fabs(d)>fabs(od) && n>99 && x!=xMin && x!=xMax)
4386 if(it>nitMax)
G4cout<<
"G4Q::GetQPMom: a#of iterations > nitMax="<<nitMax<<
G4endl;
4391 if(kCand>=kMax)kCand=kMax-.001;
4392 if(kCand<=kMin)kCand=kMin+.001;
4432 G4double qMOverT = qMass/Temperature;
4433 G4int valc = valQ.GetTot();
4459 nOfQ = RandomPoisson((1.+sqrt(1.+qMOverT*qMOverT))/2.);
4461 if (!ev && nOfQ<2) nOfQ=2;
4462 else if ( ev && nOfQ<3) nOfQ=3;
4465 G4cout<<
"G4Q::Calc#ofQP:QM="<<q4Mom<<qMass<<
",T="<<Temperature<<
",QC="<<valQ<<
",n="<<nOfQ
4468 G4int absb = abs(valQ.GetBaryonNumber());
4470 if(absb)tabn=3*absb;
4472 if (nOfQ<tabn) nOfQ=tabn;
4473 G4int nSeaPairs = (nOfQ-valc)/2;
4474 G4int stran = abs(valQ.GetS());
4475 G4int astra = abs(valQ.GetAS());
4476 if(astra>stran) stran=astra;
4477 G4int nMaxStrangeSea=
static_cast<int>((qMass-stran*mK0)/(mK0+mK0));
4478 if (absb) nMaxStrangeSea=
static_cast<int>((qMass-absb)/672.);
4480 G4cout<<
"G4Q::Calc#ofQP:"<<valQ<<
",INtot="<<valc<<
",nOfQ="<<nOfQ<<
",SeaPairs="<<nSeaPairs
4488 if(nSeaPairs>0)valQ.IncQAQ(nSeaPairs,SSin2Gluons);
4490 else morDec=valQ.DecQAQ(-nSeaPairs);
4491 if(morDec)
G4cout<<
"G4Q::Calc#ofQP: "<<morDec<<
" pairs can be reduced more"<<
G4endl;
4493 G4int sSea=valQ.GetS();
4494 G4int asSea=valQ.GetAS();
4495 if(asSea<sSea) sSea=asSea;
4496 if(sSea>nMaxStrangeSea)
4499 G4cout<<
"G4Q::Calc#ofQP:**Reduce** S="<<sSea<<
",aS="<<asSea<<
",maxS="<<nMaxStrangeSea
4502 sSea-=nMaxStrangeSea;
4505 valQ.IncQAQ(sSea,0.);
4514 G4cout<<
"G4Quasmon::Calc#ofQP: *** RESULT IN*** nQ="<<nOfQ<<
", FinalQC="<<valQ<<
G4endl;
4520 void G4Quasmon::ModifyInMatterCandidates()
4533 for (
unsigned ind=0; ind<theQCandidates.size(); ind++)
4543 if(cPDG>80000000&&cPDG!=90000000)
4545 if(totMass<tmpTM+frM)
4548 G4cout<<
"G4Q::ModInMatCand:C="<<cPDG<<tmpT<<tmpTM<<
"+"<<frM<<
"="<<tmpTM+frM<<
">tM="
4554 G4int cP = cNuc.GetZ();
4555 G4int cN = cNuc.GetN();
4556 G4int cL = cNuc.GetS();
4559 if(cPDG==90001000)
G4cout<<
"G4Q::MIM:->>cPDG=90001000<<-,possibility="<<poss<<
G4endl;
4561 if(eP>=cP&&eN>=cN&&eL>=cL&&poss)
4566 if(cP==eP&&cN==eN&&cL==eL)clME=cQPDG.GetMass();
4569 renvM = cQPDG.GetNuclMass(eP-cP,eN-cN,eL-cL);
4572 if(cP==tP&&cN==tN&&cL==tL)clMN=cQPDG.GetMass();
4575 renvM = cQPDG.GetNuclMass(tP-cP,tN-cN,tL-cL);
4583 G4cout<<
"G4Q:ModInMatCand:C="<<cPDG<<cNuc<<clME<<
","<<clMN<<
",E="<<envPDGC<<
",M="
4593 void G4Quasmon::CalculateHadronizationProbabilities
4613 if(vap>0)vaf = var*mQk/kVal/vap;
4615 G4double accumulatedProbability = 0.;
4616 G4double secondAccumProbability = 0.;
4617 G4int qBar =valQ.GetBaryonNumber();
4618 G4int nofU = valQ.GetU()- valQ.GetAU();
4619 G4int dofU = nofU+nofU;
4620 G4int nofD = valQ.GetD()- valQ.GetAD();
4621 G4int dofD = nofD+nofD;
4622 G4int qChg = valQ.GetCharge();
4623 G4int qIso = qBar-qChg-qChg;
4626 G4int maxC = theQCandidates.size();
4627 G4double totZ = theEnvironment.
GetZ() + valQ.GetCharge();
4638 G4int absb = abs(qBar);
4640 G4cout<<
"G4Q::CalcHadronizationProbab:Q="<<mQ<<valQ<<
",v="<<vaf<<
",r="<<var<<
",mC="<<maxB
4641 <<
",vap="<<vap<<
",k="<<kVal<<
G4endl;
4644 unsigned nHC=theQCandidates.size();
4652 G4int aPDG = abs(cPDG);
4656 if ( (aPDG > 80000000 && envA > 0) || aPDG < 80000000)
4666 G4int cU=candQC.GetU()-candQC.GetAU();
4668 G4int cD=candQC.GetD()-candQC.GetAD();
4670 G4int dUD=abs(cU-cD);
4675 G4bool pPrint= (abs(cPDG)%10 <3 && cPDG <80000000) || (cPDG >80000000 && frM <5000.);
4681 if(pPrint)
G4cout<<
"G4Q::CHP:==****==>> c="<<cPDG<<
",dUD="<<dUD<<
",pos="<<pos<<
",eA="
4682 <<envA<<
",tM="<<totMass<<
" > tmpTM+frM="<<tmpTM+frM<<
G4endl;
4685 if(pos&&(cPDG<80000000||(cPDG>80000000&&cPDG!=90000000&&dUD<2)))
4692 G4int cC = candQC.GetCharge();
4697 G4int cNQ= candQC.GetTot()-2;
4702 if(pPrint)
G4cout<<
"G4Q::CHP:B="<<baryn<<
",C="<<cC<<
",CB="<<CB<<
",#q="<<cNQ<<
G4endl;
4704 if(cPDG>80000000&&cPDG!=90000000&&baryn<=envA)
4706 G4int parentCounter=0;
4716 if (dofU<=nofD) iQmax=1;
4717 else if(dofD<=nofU) iQmin=1;
4725 if(pPrint)
G4cout<<
"G4Q::CHP:***!!!***>>F="<<cPDG<<
",mF="<<frM<<
",iq:"<<iQmin<<
","
4726 <<iQmax<<
",kLS="<<kLS<<
",kQCM="<<kVal<<
",eA="<<envA<<
G4endl;
4728 if(iQmax>iQmin)
for(
int iq=iQmin; iq<iQmax; iq++)
4736 if(pPrint)
G4cout<<
"G4Q::CHP:photon cap(gaF) is enhanced for Uquark"<<
G4endl;
4740 if (!iq) nqInQ=valQ.GetD();
4741 else if(iq==1) nqInQ=valQ.GetU();
4742 else if(iq==2) nqInQ=valQ.GetS();
4745 G4cout<<
"G4Q::CHP:i="<<iq<<
",cU="<<cU<<
",cD="<<cD<<
",omi="<<oQmin<<
",oma="
4748 if(oQmax>oQmin)
for(
int oq=oQmin; oq<oQmax; oq++)
4750 G4int shift= cQPDG.GetRelCrossIndex(iq, oq);
4754 G4cout<<
"G4Q::CHP:iq="<<iq<<
",oq="<<oq<<
",QC="<<ioQC<<
",rQC="<<resQC<<
G4endl;
4757 resPDG=resQPDG.GetPDGCode();
4758 G4int resQ=resQPDG.GetQCode();
4760 if(pPrint)
G4cout<<
"G4Q::CHP:i="<<iq<<
",o="<<oq<<ioQC<<
",s="<<shift
4761 <<
",cQPDG="<<cQPDG<<
", residQC="<<resQC<<resQPDG<<
G4endl;
4767 G4bool rI=resA>0 && resU>=0 && resD>=0 &&
4768 (resU+resS>resD+resD||resD+resS>resU+resU);
4775 if (resQ > -2 && resPDG && resPDG != 10 && !rI &&
4776 (!piF || ( piF && (cPDG != 90001000 ||
G4UniformRand() < .333333) &&
4777 cPDG != 90002001 && cPDG != 90002002
4790 if(shift!=7&&is<maxC)
4798 if(pPrint)
G4cout<<
"G4Q::CHP:parentPossibility="<<possib<<
",pZ="<<charge
4799 <<
" <= envZ="<<envZ<<
", pN="<<barot-charge<<
" <= envN="
4800 <<envN<<
", cPDG="<<cPDG<<
G4endl;
4802 if(possib && charge<=envZ && barot-charge<=envN)
4807 G4int isos = barot-charge-charge;
4811 if(barot>2) pUD= pow(2.,abs(isos)-1);
4814 if(barot!=baryn)
G4cerr<<
"--Warning--G4Q::CHP:c="<<candQC<<
",p="<<parQC
4820 if(pPrint)
G4cout<<
"G4Q::CHP: dS="<<dS<<
", dI="<<dI<<
", dC="<<dC<<
", I="
4821 <<qIso<<
",i="<<isos<<
", C="<<cC<<
",c="<<charge<<
G4endl;
4869 if ( (!piF && baryn < 5 ) ||
4870 ( piF && abs(dS) < 3) )
4882 if(pPrint)
G4cout<<
"G4Q::CHP:c="<<cPDG<<
",p="<<parPDG<<
",bM="<<boundM
4883 <<
",i="<<is<<
",adr="<<parCand<<
",pPP="<<pPP<<
G4endl;
4891 if(nucBM)nDelta=(frM2-bNM2)/(nucBM+nucBM);
4893 G4int iniPDG =valQ.GetSPDGCode();
4896 G4double kCut=boundM/2.+freeE/(iniQM+boundM);
4897 if(pPrint)
G4cout<<
"G4Q::CHP:r="<<resPDG<<
",M="<<minM<<
",k="<<kLS<<
"<"
4898 <<kCut<<
",E="<<E<<
">"<<nDelta<<
",p="<<parPDG<<
G4endl;
4900 if(resPDG && minM>0.)
4903 if(pPrint)
G4cout<<
"G4Q::CHP:fM="<<frM<<
",bM="<<boundM<<
",rM="
4904 <<tmpTM<<
",tM="<<totMass<<
G4endl;
4909 G4double eDelta=(frM2-bM2)/(boundM+boundM);
4919 G4double Em=(E-nDelta-CB)*(1.-frM/totMass);
4923 if(ne>0.&&ne<1.)kf=pow(ne,cNQ);
4925 if(pPrint)
G4cout<<
"G4Q::CHP:<qi_DEF>="<<ne<<
",k="<<kf<<
",dk="<<dked
4926 <<
",dkLS="<<dkLS<<
",M="<<boundM<<
",C="<<CB<<
G4endl;
4931 if(envA-barot>bEn) rtQC+=bEnQC;
4932 else rtQC+=envQC-parQC;
4937 if(envA-barot>bEn) rtEP+=mbEn;
4938 else rtEP+=envM-boundM;
4943 if(pPrint)
G4cout<<
"G4Q::CHP:RN="<<tmpEQ<<
"="<<envM-boundM<<
"=eM="
4944 <<envM<<
"-bM="<<boundM<<
",E="<<rtE<<
",eQC="<<envQC
4945 <<
",pQC="<<parQC<<
G4endl;
4952 if ( (envA-barot <= bEn && envM > boundM) || envA-barot > bEn)
4957 if(envA-barot > bEn) minBM-=mbEn;
4958 else minBM-=envM-boundM;
4964 if(pPrint)
G4cout<<
"G4Q::CHP:M2="<<minM2<<
",R="<<rQ2<<
",m="<<mintM2
4965 <<
",RN2="<<rtQ2<<
",q="<<(minM2-rQ2)/rEP/2<<
",qN="
4966 <<(mintM2-rtQ2)/rtEP/2<<
G4endl;
4971 G4double nc=1.-(dkLS-E-E+CB+CB)/boundM;
4974 if(pPrint)
G4cout<<
"G4Q::CHP:qi_k-E="<<nc<<
",k="<<kLS<<
",E="<<E
4977 if(nc > 0. && nc < 1. && nc < ne)
4981 if(newh < kf) kf=newh;
4983 else if(nc <= 0.) kf=0.;
4987 if(pPrint)
G4cout<<
"G4Q::CHP:qi_R="<<nk<<
",kd="<<kd<<
",E="<<Em
4988 <<
",M="<<totMass<<
G4endl;
4990 if(nk > 0. && nk < 1. && nk < ne)
4994 if(newh<kf) kf=newh;
4996 else if(nk <= 0.) kf=0.;
5018 if(mix > frM) st=sqrt(mix*mix-frM2);
5019 G4double nq=1.-(dkLS-st-st)/boundM;
5021 if(pPrint)
G4cout<<
"G4Q::CHP:qi_k-st="<<nq<<
",st="<<st<<
",m="
5022 <<mix<<
",M="<<frM<<
G4endl;
5024 if(nq > 0. && nq < 1. && nq < ne)
5029 if(newh < kf) kf=newh;
5031 else if(nq<=0.)kf=0.;
5036 G4bool atrest=(eQ-mQ)/mQ<.001||k3V.
dot(rq3V)<-.999;
5043 if(atrest) nz=1.-(mintM2-rtQ2+pmk*dked)/(boundM*(rtEP+pmk));
5044 else nz=1.-(mintM2-rtQ2)/(boundM*rtEP);
5048 if(pPrint)
G4cout<<
"G4Q::CHP:q="<<nz<<
",a="<<atrest<<
",M2="
5049 <<mintM2<<
">"<<rtQ2<<
G4endl;
5051 if(nz > 0. && nz < 1. && nz < ne)
5055 if(newh < kf) kf=newh;
5057 else if(nz <= 0.) kf=0.;
5082 (piF && cPDG!=90000001 && cPDG!=90001001 && cPDG!=90001002)
5089 if(atrest) nz=1.-(minBM2-rQ2+pmk*dked)/(boundM*(rEP+pmk));
5090 else nz=1.-(minBM2-rQ2)/(boundM*rEP);
5094 if(pPrint)
G4cout<<
"G4Q::CHP:q="<<nz<<
",a="<<atrest<<
",QM2="
5095 <<minM2<<
">"<<rQ2<<
G4endl;
5097 if(nz>0.&&nz<1.&&nz<ne)
5101 if(newh<kf) kf=newh;
5103 else if(nz<=0.)kf=0.;
5122 ( (!piF && baryn > 4) ||
5124 (piF && cPDG != 90000001
5125 && cPDG != 90001001)
5133 if(atrest) nz=1.-(minM2-rQ2+pmk*dked)/(boundM*(rEP+pmk));
5134 else nz=1.-(minM2-rQ2)/(boundM*rEP);
5138 if(pPrint)
G4cout<<
"G4Q::CHP:q="<<nz<<
",a="<<atrest<<
",QM2="
5139 <<minM2<<
">"<<rQ2<<
G4endl;
5141 if(nz>0.&&nz<1.&&nz<ne)
5145 if(newh<kf) kf=newh;
5147 else if(nz<=0.)kf=0.;
5153 if(pPrint)
G4cout<<
"G4Q::CHP:"<<kf<<
",minM2="<<minM2<<
",rQ2="<<rQ2
5162 if(lz>0.&&lz<1.)low=pow(lz,cNQ);
5163 else if(lz>=1.)low=1.;
5166 if(pPrint)
G4cout<<
"G4Q::CHP:<qa_DEF>="<<lz<<
", eDel="<<eDelta
5167 <<
",nDel="<<nDelta<<
G4endl;
5174 if(pPrint)
G4cout<<
"G4Q::CHP:qa_t="<<le<<
",k="<<kLS<<
",E="<<Em
5175 <<
",bM="<<boundM<<
G4endl;
5177 if(le>0.&&le<1.&&le>lz)
5181 if(newl>low) low=newl;
5183 else if(le>=1.)low=1.;
5187 G4double lk=1.-(dkLS+E+E-CB-CB)/boundM;
5189 if(pPrint)
G4cout<<
"G4Q::CHP:qa_k+E="<<lk<<
",k="<<kLS<<
",E="<<E
5192 if(lk>0.&&lk<1.&&lk>lz)
5196 if(newl>low) low=newl;
5198 else if(lk>=1.)low=1.;
5218 if(pPrint)
G4cout<<
"G4Q::CHP:qa_k+sr="<<lp<<
",sr="<<
sr
5221 if(lp>0.&&lp<1.&&lp>lz)
5225 if(newl>low) low=newl;
5227 else if(lp>=1.)low=1.;
5235 G4double nz=1.-(qmaCB+qmaCB)/boundM;
5237 if(pPrint)
G4cout<<
"G4Q::CHP:<qa_CB>="<<nz<<
",m="<<qmaCB<<
",CB="
5243 if(newl>low) low=newl;
5245 else if(nz>1.) low=10.;
5250 if(pPrint)
G4cout<<
"G4Q::CHP:>>"<<cPDG<<
",l="<<low<<
",h="<<high
5251 <<
",ol="<<pl<<
",oh="<<ph<<
",nl="<<newl<<
",nh="
5252 <<newh<<
",kf="<<kf<<
",d="<<eDelta<<
G4endl;
5258 G4int noc=cQPDG.GetNumOfComb(iq, oq);
5261 probab=qFact*kf*nqInQ*pPP*noc/pUD;
5272 if(BarRQC==2 && !StrRQC)
5275 if(ChgRQC==1) probab/=2;
5279 if(pPrint)
G4cout<<
"G4Q::CHP:prob="<<probab<<
",qF="<<qFact<<
",iq="
5280 <<iq<<
",oq="<<oq<<
",Pho4M="<<phot4M<<
",pUD="<<pUD
5283 if(probab<0.) probab=0.;
5296 G4cout<<
"G4Q::CalcHP: FillParentClaster="<<*curParC<<
G4endl;
5301 if(pPrint)
G4cout<<
"G4Q::CHP:in="<<
index<<
",cPDG="<<cPDG<<
",pc"<<parentCounter
5302 <<parQC<<
",Env="<<theEnvironment<<
",comb="<<comb
5309 else G4cout<<
"***G4Q::CHP:dI="<<dI<<
",cC="<<cC<<
G4endl;
5316 if(pPrint)
G4cout<<
"G4Q::CHPr: probab="<<probability<<
"("<<comb<<
"),iq="<<iq
5322 else if(cPDG<80000000)
5325 G4int curnh=theQHadrons.size();
5329 for (
G4int ind=0; ind<curnh; ind++)
5331 G4int curhPDG=theQHadrons[ind]->GetPDGCode();
5332 if (curhPDG== 111) npiz++;
5333 if (curhPDG== 211) npip++;
5334 if (curhPDG==-211) npin++;
5337 comb = valQ.NOfCombinations(candQC);
5340 if ( (aPDG==111)|(aPDG==211) ) comb=1.;
5341 else if ( (aPDG==311)|(aPDG==321) ) comb=SSin2Gluons;
5343 if(cPDG== 211&&npip>0) comb*=(npip+1);
5344 if(cPDG==-211&&npip>0) comb*=(npin+1);
5345 if(cPDG==111||cPDG==221||cPDG==331||cPDG==113||cPDG==223||cPDG==333||cPDG==115||
5346 cPDG==225||cPDG==335||cPDG==117||cPDG==227||cPDG==337||cPDG==110||cPDG==220||
5352 G4double cmd=valQ.NOfCombinations(tQCd);
5353 G4double cmu=valQ.NOfCombinations(tQCu);
5354 G4double cms=valQ.NOfCombinations(tQCs);
5355 if(cPDG!=333&&cPDG!=335&&cPDG!=337) comb=(cmd+cmu)/2.;
5357 if(cPDG==331||cPDG==221) comb =(comb + cms)/4.;
5358 if(cPDG==113) comb*=4.;
5359 if(cPDG==223) comb*=2.;
5360 if(cPDG==111&&npiz>0) comb*=(npiz+1);
5362 if(abs(cPDG)<3) G4cout<<
"G4Q::CHP:comb="<<comb<<
",cmd="<<cmd<<
",cmuu="<<cmu
5363 <<
",cms="<<cms<<G4endl;
5371 G4bool priCon = aPDG < 10000 && aPDG%10 < 3;
5372 if(priCon)
G4cout<<
"G4Q::CHP:***>>cPDG="<<cPDG<<
",cQC="<<candQC<<
",comb="<<comb
5373 <<
",curQC="<<curQ<<
",mQ="<<mQ<<
",ab="<<absb<<
G4endl;
5375 if(resPDG==221 || resPDG==331)
5381 if(priCon)
G4cout<<
"G4Q::CHP:cPDG="<<cPDG<<
",c="<<comb<<
",rPDG/QC="<<resPDG<<curQ
5382 <<
",tM="<<totMass<<
">"<<frM-CB+resTM<<
"=fM="<<frM<<
"+rM="<<resTM
5385 if (comb && resPDG && totMass > frM-CB+resTM &&
5386 ((resPDG > 80000000 && resPDG != 90000000) || resPDG<10000) )
5389 if(priCon)
G4cout<<
"G4Q::CHP:ind="<<
index<<
",qQC="<<valQ<<mQ<<
",cPDG="<<cPDG
5390 <<
",rPDG="<<resPDG<<curQ<<
G4endl;
5396 if(priCon)
G4cout<<
"G4Q::CHP:rM/QC="<<resM<<curQ<<
",E="<<envPDGC<<
",rQC="
5400 if(envPDGC>80000000 && envPDGC!=90000000 && resM>0. &&
5401 resPDG!=10 && resPDG!=1114 && resPDG!=2224)
5408 if(priCon)
G4cout<<
"G4Q::CHP: **Rec**,RQMass="<<bnRQ<<
",envM="<<envM<<
",rtM="
5412 if(bnRQ<resM) resM=bnRQ;
5415 if(aPDG<10000&&aPDG%10<3)
5417 G4cout<<
"G4Q::CHP: resM="<<resM<<
", resQCode="<<resQCode<<
G4endl;
5419 if(resM>0. && resQCode>-2)
5422 G4double rndM=GetRandomMass(cPDG,limM);
5425 if(aPDG<10000&&aPDG%10<3)
5427 G4cout<<
"G4Q::CHP:rndM="<<rndM<<
",limM="<<limM<<
" > cM="<<cMass<<
" ,rM+fM="
5428 <<resM+rndM<<
" < mQ="<<mQ<<
G4endl;
5431 if(rndM>0. && resM+rndM<mQ)
5442 if(qBar) zMin= mR2/mQ/(mQ-dk);
5445 if(priCon)
G4cout<<
"G4Q::CHP:M="<<rndM<<
",ps="<<possibility<<
",zMax="<<zMax
5446 <<
",rHk="<<rHk<<
",mQ="<<mQ<<
",dk="<<dk<<
",zMin="<<zMin
5447 <<
",mR2="<<mR2<<
",rM="<<resM<<
"; "<<mQ*(mQ-dk)<<G4endl;
5452 if(rM2<resM*resM) possibility = 0.;
5454 if (possibility>0. && vap>0 && zMax>zMin)
5456 probability = vaf*(pow(zMax, vap)-pow(zMin, vap));
5458 if(priCon)
G4cout<<
"G4Q::CHP:#"<<
index<<
",mH2="<<mH2<<
",nQ="<<nOfQ<<
",p="
5459 <<probability<<
",vf="<<vaf<<
",vp="<<vap<<
",zMax="<<zMax
5460 <<
",zMin="<<zMin<<
G4endl;
5472 if(cPDG==110||cPDG==220||cPDG==330) probability*=comb;
5473 else probability*=comb*(abs(cPDG)%10);
5476 if(BarRQC==2 && !StrRQC)
5479 if(ChgRQC==1) probability/=2;
5480 else probability*=2;
5488 if(priCon)
G4cout<<
"G4Q::CHP:cM=0[cPDG"<<cPDG<<
"],mQ/QC="<<mQ<<valQ<<
",rM="
5496 if(priCon)
G4cout<<
"***G4Q::CHP: M=0, #"<<
index<<valQ<<
",cPDH="<<cPDG<<
"+rP="
5505 if(priCon)
G4cout<<
"G4Q::CHP:"<<
index<<valQ<<
",PDG="<<cPDG<<
"+r="<<resPDG<<curQ
5506 <<
":c=0(!) || tM="<<totMass<<
"<"<<frM-CB+resTM<<
" = fM="<<frM
5507 <<
"+rTM="<<resTM<<
"-CB="<<CB<<
G4endl;
5512 else probability=0.;
5514 G4int aPDG = abs(cPDG);
5515 if((cPDG>90000000 && baryn<5) || (aPDG<10000 && aPDG%10<3))
5516 G4cout<<
"G4Q::CHP:^^cPDG="<<cPDG<<
",p="<<pos<<
",rPDG="<<resPDG<<curQ<<resM<<
",p="
5517 <<probability<<
",a="<<accumulatedProbability<<
",sp="<<secondProbab<<
G4endl;
5522 accumulatedProbability += probability;
5525 secondAccumProbability += secondProbab;
5547 G4int resQPDG=valQ.GetSPDGCode();
5555 G4cout<<
"G4Q::CheckGS: EnvPDG="<<theEnvironment.
GetPDG()<<
",Quasmon="<<resQPDG<<
G4endl;
5557 if(theEnvironment.
GetPDG()!=90000000)
5559 resEMa=theEnvironment.
GetMZNS();
5561 G4cout<<
"G4Q::CheckGS: Environment Mass="<<resEMa<<
G4endl;
5571 bsCond = tmpQN.SplitBaryon();
5573 G4cout<<
"G4Q::CheckGS: No environment, theOnlyQ="<<tmpQN<<
",bsCond="<<bsCond<<
G4endl;
5584 else if(bsCond==3122)
5589 else if(bsCond==90001001)
5594 else if(bsCond==90002002)
5606 if(
G4QHadron(reTLV).DecayIn2(frag4M,quas4M))
5609 FillHadronVector(quasH);
5611 FillHadronVector(fragH);
5624 G4int nOfOUT = theQHadrons.size();
5626 G4cout<<
"G4Q::CheckGS: (totM="<<resTMa<<
" < rQM+rEM="<<resSMa<<
" || rEM="<<resEMa
5627 <<
"=0 && "<<bsCond<<
"=0) && n="<<nOfOUT<<
" >0"<<
G4endl;
5629 if ( (resTMa < resSMa || (!resEMa && !bsCond) ) && nOfOUT > 0 && corFlag)
5632 G4QHadron* theLast = theQHadrons[nOfOUT-1];
5639 G4cout<<
"G4Q::CheckGS:YES,T="<<tmpTLV<<tmpTLV.
m()<<
">rM+hM="<<resSMa+hadrMa<<
G4endl;
5641 if(tmpTLV.
m()>resSMa+hadrMa)
5646 if(!
G4QHadron(tmpTLV).DecayIn3(hadr4M,quas4M,enva4M))
5648 G4cerr<<
"---Warning---G4Q::CheckGS:DecIn Fragm+ResQ+ResEnv Error"<<
G4endl;
5656 FillHadronVector(envH);
5657 theEnvironment = vacuum;
5660 FillHadronVector(quasH);
5668 if(!
G4QHadron(tmpTLV).DecayIn2(hadr4M,quas4M))
5671 G4cerr<<
"---Warning---G4Q::CheckGS: Decay in Fragm+ResQ Error"<<
G4endl;
5679 FillHadronVector(quasH);
5685 if(nOfOUT>1 && corFlag)
5687 G4QHadron* thePrev = theQHadrons[nOfOUT-2];
5697 G4cout<<
"G4Q::CheckGS:NO, M="<<tmpTLV<<totNMa<<
">"<<totQMa+hadrMa+prevMa<<
G4endl;
5699 if(totNMa>totQMa+hadrMa+prevMa)
5702 if(!
G4QHadron(tmpTLV).DecayIn3(hadr4M,prev4M,nuc4M))
5704 G4cerr<<
"---Warning---G4Q::CheckGS:DecIn3 ResN+Last+Prev Error"<<
G4endl;
5711 FillHadronVector(envH);
5712 theEnvironment = vacuum;
5716 G4cout<<
"G4Q::CheckGS: Yes, Check D4M="<<tmpTLV-hadr4M-prev4M-nuc4M<<
G4endl;
5730 G4cout<<
"G4Quasm::CheckGS: NO, tM="<<totNMa<<
" > rp+l="<<resLastM+hadrMa
5731 <<
" || > rl+p="<<resPrevM+prevMa<<
G4endl;
5733 if (totNMa>resLastM+hadrMa)
5735 theQHadrons.pop_back();
5736 theQHadrons.pop_back();
5737 theQHadrons.push_back(theLast);
5744 else if (totNMa>resPrevM+prevMa)
5746 theQHadrons.pop_back();
5751 if(!
G4QHadron(tmpTLV).DecayIn2(prev4M,nuc4M))
5753 G4cerr<<
"---Warning---G4Q::CheckGS:DecIn2 (ResN+Last)+Prev Error"<<
G4endl;
5760 FillHadronVector(envH);
5761 theEnvironment = vacuum;
5764 G4cout<<
"G4Q::CheckGS:Yes, Check D4M="<<tmpTLV-prev4M-nuc4M<<
G4endl;
5785 if(thePDG<0) pap = 1;
5790 G4int nChan = decV.size();
5792 G4cout<<
"G4Quasm::DecQHadron: PDG="<<thePDG<<
",m="<<m_value<<
",("<<nChan<<
" channels)"<<
G4endl;
5800 for(i=0; i<nChan; i++)
5807 if(rnd<dC->GetDecayChanLimit() && m_value>dC->
GetMinMass())
break;
5809 if(i>nChan-1) i=nChan-1;
5812 G4int nPart=cV.size();
5814 G4cout<<
"G4Quasmon::DecayQHadron: resi="<<i<<
",nP="<<nPart<<
":"<<cV[0]->GetPDGCode()
5815 <<
","<<cV[1]->GetPDGCode();
5816 if(nPart>2)
G4cout<<
","<<cV[2]->GetPDGCode();
5819 if(nPart<2||nPart>3)
5821 G4cerr<<
"---Warning---G4Q::DecayQHadr:n="<<nPart<<
",ch#"<<i<<
",PDG="<<thePDG<<
G4endl;
5822 theFragments->push_back(qH);
5823 return theFragments;
5826 G4cout<<
"G4Q::DecQH:Decay("<<ElMaDecays<<
") PDG="<<thePDG<<t<<m_value<<
",nP="<<nPart<<
G4endl;
5833 G4int sPDG=cV[1]->GetPDGCode();
5835 if ( (fPDG != 22 && sPDG != 22) || ElMaDecays) {
5837 G4cout<<
"G4Q::DecQH:Yes2,fPDG="<<fPDG<<
",sPDG="<<sPDG<<
",EMF="<<ElMaDecays<<
G4endl;
5839 if(cV[0]->GetWidth()==0.)
5841 fHadr =
new G4QHadron(cV[0]->GetPDGCode());
5842 if(cV[1]->GetWidth()==0.)sHadr =
new G4QHadron(sPDG);
5864 G4cout<<
"G4Q::DQH:M="<<m_value<<
",mi="<<mim<<
",fd="<<fdm<<
",fm="<<fm<<
",sd="<<sdm
5883 G4cerr<<
"---Warning---G4Q::DecayQHadron:in2,PDGC="<<thePDG<<
", ch#"<<i<<
": 4M="
5887 theFragments->push_back(qH);
5888 return theFragments;
5899 G4int nProd=theTmpQHV->size();
5901 G4cout<<
"G4Q::DecayQHadr:(DecayIn2) nOfProdForQH1="<<nProd<<
G4endl;
5903 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
5904 else for(
G4int ip1=0; ip1<nProd; ip1++)
5907 G4int tmpS=intTmpQHV->size();
5908 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
5911 theFragments->resize(tmpS+theFragments->size());
5912 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
5915 G4cout<<
"G4Q::DecayQHadr:(DecayIn2) Copy Sec11 nProd="<<tmpS<<
G4endl;
5924 nProd=theTmpQHV->size();
5926 G4cout<<
"G4Q::DecayQHadr:(DecayIn2) nOfProdForQH2="<<nProd<<
G4endl;
5928 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
5929 else for(
G4int ip1=0; ip1<nProd; ip1++)
5932 G4int tmpS=intTmpQHV->size();
5933 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
5936 theFragments->resize(tmpS+theFragments->size());
5937 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
5940 G4cout<<
"G4Q::DecayQHadr:(DecayIn2) Copy Sec12 nProd="<<tmpS<<
G4endl;
5949 G4cout<<
"G4Q::DecQHadr: DecayIn2 is made with nH="<<theFragments->size()<<
G4endl;
5955 if(thePDG==89999003||thePDG==90002999)
G4cerr<<
"*G4Q::DQH:8999003/90002999"<<
G4endl;
5957 theFragments->push_back(qH);
5966 G4int sPDG=cV[1]->GetPDGCode();
5967 G4int tPDG=cV[2]->GetPDGCode();
5969 if ( (fPDG != 22 && sPDG != 22 && tPDG != 22) || ElMaDecays)
5972 G4cout<<
"G4Q::DQH:Y,f="<<fPDG<<
",s="<<sPDG<<
",t="<<tPDG<<
",F="<<ElMaDecays<<
G4endl;
5974 if(cV[0]->GetWidth()==0.)
5977 if(cV[1]->GetWidth()==0.)
5980 if(cV[2]->GetWidth()==0.)tHadr =
new G4QHadron(tPDG);
6009 G4double fdm = m_value - smim - tmim;
6022 G4cout<<
"G4Quasmon::DecayQHadron:3Dec. m1="<<fHadr->
GetMass()
6034 if(!qH->
DecayIn3(f4Mom,s4Mom,t4Mom))
6039 G4cerr<<
"---Warning---G4Q::DecayQHadron:in3,PDGC="<<thePDG<<
", ch#"<<i<<
G4endl;
6040 theFragments->push_back(qH);
6041 return theFragments;
6052 G4int nProd=theTmpQHV->size();
6054 G4cout<<
"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH1="<<nProd<<
G4endl;
6056 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
6057 else for(
G4int ip1=0; ip1<nProd; ip1++)
6060 G4int tmpS=intTmpQHV->size();
6061 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
6064 theFragments->resize(tmpS+theFragments->size());
6065 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
6068 G4cout<<
"G4Q::DecayQHadr:(DecayIn3) Copy Sec11 nProd="<<tmpS<<
G4endl;
6078 nProd=theTmpQHV->size();
6080 G4cout<<
"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH2="<<nProd<<
G4endl;
6082 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
6083 else for(
G4int ip1=0; ip1<nProd; ip1++)
6086 G4int tmpS=intTmpQHV->size();
6087 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
6090 theFragments->resize(tmpS+theFragments->size());
6091 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
6094 G4cout<<
"G4Q::DecayQHadr:(DecayIn3) Copy Sec12 nProd="<<tmpS<<
G4endl;
6104 nProd=theTmpQHV->size();
6106 G4cout<<
"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH3="<<nProd<<
G4endl;
6108 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
6109 else for(
G4int ip1=0; ip1<nProd; ip1++)
6112 G4int tmpS=intTmpQHV->size();
6113 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
6116 theFragments->resize(tmpS+theFragments->size());
6117 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
6120 G4cout<<
"G4Q::DecayQHadr:(DecayIn3) Copy Sec13 nProd="<<tmpS<<
G4endl;
6130 G4cout<<
"G4Q::DecQHadr: DecayIn3 is made with nH="<<theFragments->size()<<
G4endl;
6133 else theFragments->push_back(qH);
6139 G4cout<<
"G4Quas::DecQHadr:Fill PDG= "<<thePDG<<t<<m_value<<
" as it is ***0***>>"<<
G4endl;
6141 if(thePDG==89999003||thePDG==90002999)
G4cerr<<
"-War-G4Q::DQH:8999003/90002999"<<
G4endl;
6142 theFragments->push_back(qH);
6145 G4cout<<
"G4Q::DecQHadr:=-= HADRON IS DECAYED =-= with nH="<<theFragments->size()<<
G4endl;
6147 return theFragments;
6155 G4cerr<<
"---Warning---G4Q::RandomPoisson:Negative(zero) MeanValue="<<meanValue<<
G4endl;
6162 if (r<s_value)
return 0;
6165 if (r<s_value)
return 1;
6167 while ( s_value<r && i<100 )
6184 HadronizeQuasmon(nucEnviron,nQs);
6185 G4int nHadrs=theQHadrons.size();
6187 G4cout<<
"G4Quasm::Fragm:after HadronizeQuasmon nH="<<nHadrs<<
",Env="<<nucEnviron<<
G4endl;
6190 if(nHadrs)
for (
int hadron=0; hadron<nHadrs; hadron++)
6193 theFragments->push_back(curHadr);
6196 else G4cerr<<
"*******G4Quasmon::Fragment *** Nothing is in the output ***"<<
G4endl;
6198 return theFragments;
6207 q4Mom.
setE(q4Mom.
dot(boost4M)/bm);