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());
179G4double G4Quasmon::Temperature=180.;
182G4bool G4Quasmon::WeakDecays=
false;
183G4bool 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);
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;
941 for (
G4int index=0; index<totCand; index++)
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)
964 G4cout<<
"G4Q::HQ:cfM="<<cfM<<
",cMs="<<cMs<<
",ind="<<index<<
G4endl;
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="
1586 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;
3792void 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;
4520void 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="
4593void 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();
4648 if(nHC)
for (
unsigned index=0; index<nHC; index++)
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
4789 G4int is=index+shift;
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
4815 <<
",s="<<shift<<
",i="<<index<<
",s="<<is<<
G4endl;
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.;
5216 G4double lp=1.-(dkLS+sr+sr)/boundM;
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
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);
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4QDecayChan * > G4QDecayChanVector
std::vector< G4QHadron * > G4QHadronVector
std::vector< G4QPDGCode * > G4QPDGCodeVector
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
double dot(const Hep3Vector &) const
Hep3Vector boostVector() const
double dot(const HepLorentzVector &) const
void setVect(const Hep3Vector &)
G4QParticle * GetQParticle(G4int PDG) const
static G4QCHIPSWorld * Get()
G4int GetQPEntries() const
G4double GetNBMass() const
G4bool GetParPossibility() const
void SetSecondIntProb(G4double intP)
G4bool GetPossibility() const
void SetPossibility(G4bool choice)
void SetIntegProbability(G4double intP)
G4double GetPreProbability() const
void SetParPossibility(G4bool choice)
void FillPClustVec(G4QParentCluster *pCl)
void SetEBMass(G4double newMass)
G4double GetEBMass() const
void SetRelProbability(G4double relP)
G4int GetPClustEntries() const
G4QParentCluster * TakeParClust(G4int nPC)
void SetSecondRelProb(G4double relP)
void ClearParClustVector()
void SetNBMass(G4double newMass)
G4int DecQAQ(const G4int &nQAQ=1)
G4int GetBaryonNumber() const
G4int GetStrangeness() const
G4int GetSPDGCode() const
G4QContent SplitChipo(G4double mQ)
G4double GetDecayChanLimit() const
G4double GetMinMass() const
G4LorentzVector Get4Momentum() const
G4bool CorMDecayIn2(G4double corM, G4LorentzVector &fr4Mom)
G4bool DecayIn2(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom)
G4bool DecayIn3(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom, G4LorentzVector &t4Mom)
void SetQPDG(const G4QPDGCode &QPDG)
G4int GetNFragments() const
void Set4Momentum(const G4LorentzVector &aMom)
G4QPDGCode GetQPDG() const
G4int GetNDefMesonC() const
G4double GetProbability(G4int bn=0) const
void InitCandidateVector(G4QCandidateVector &theQCandidates, G4int nM=45, G4int nB=72, G4int nC=117)
G4QContent GetQCZNS() const
G4int GetMaxClust() const
G4double CoulBarPenProb(const G4double &CB, const G4double &E, const G4int &C, const G4int &B)
G4int GetNDefBaryonC() const
void PrepareCandidates(G4QCandidateVector &theQCandidates, G4bool piF=false, G4bool gaF=false, G4LorentzVector LV=G4LorentzVector(0., 0., 0., 0.))
G4int UpdateClusters(G4bool din)
G4double CoulombBarrier(const G4double &cZ=1, const G4double &cA=1, G4double dZ=0., G4double dA=0.)
G4double GetGSMass() const
G4QContent GetQuarkContent() const
G4double GetNuclMass(G4int Z, G4int N, G4int S)
void SetNBind(G4double bEn)
void SetLow(G4double loLim)
G4QContent GetTransQC() const
G4double GetNBind() const
G4double GetNBMass() const
G4double GetEBind() const
G4double GetProbability() const
void SetTransQC(G4QContent newTrans)
void SetEBind(G4double bEn)
void SetNQPart2(G4int nm2)
void SetEBMass(G4double bMass)
void SetHigh(G4double hiLim)
void SetNBMass(G4double bMass)
G4double GetEBMass() const
G4QDecayChanVector GetDecayVector()
G4double MinMassOfFragm()
G4QHadronVector * DecayQHadron(G4QHadron *hadron)
void Boost(const G4LorentzVector &theBoost)
G4Quasmon(G4QContent qQCont=G4QContent(0, 0, 0, 0, 0, 0), G4LorentzVector q4M=G4LorentzVector(0., 0., 0., 0.), G4LorentzVector ph4M=G4LorentzVector(0., 0., 0., 0.))
static void OpenElectromagneticDecays()
static void SetEtaSup(G4double etaetap)
static void SetSOverU(G4double ssin2g)
G4QHadronVector * Fragment(G4QNucleus &nucEnviron, G4int nQ=1)
static void CloseElectromagneticDecays()
G4int CalculateNumberOfQPartons(G4double qMass)
static void SetParameters(G4double temper=180., G4double ssin2g=.3, G4double etaetap=.3)
G4QHadronVector * DecayQuasmon()
static void SetTemper(G4double temperature)
const G4Quasmon & operator=(const G4Quasmon &right)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
std::ostringstream G4ExceptionDescription