56G4int G4QFragmentation::nCutMax=27;
59G4double G4QFragmentation::stringTension=2*GeV/fermi;
63G4double G4QFragmentation::tubeDensity =.5/fermi;
66G4double G4QFragmentation::widthOfPtSquare=-GeV*GeV;
71 static const G4double mProt2= mProt*mProt;
78 G4bool stringsInitted=
false;
87 if(cs > 0.) cM=std::sqrt(cs);
88 G4cout<<
"G4QFragmentation::Construct: *Called*, p4M="<<proj4M<<
", A="<<aNucleus<<tgMass
89 <<
",M="<<cM<<
",s="<<cs<<
",t4M="<<targ4M<<
G4endl;
93 G4int tPDG=90000000+tZ*1000+tN;
104 G4cout<<
"-EMC-G4QFragmentation::Construct:tLS4M="<<totLS4M<<
",tZLS4M="<<totZLS4M<<
G4endl;
112 G4int ModelMode=SOFT;
120 std::pair<G4double,G4double> ratios=std::make_pair(0.,0.);
121 G4int apPDG=std::abs(pPDG);
125 ratios = theQuasiElastic->
GetRatios(pMom, pPDG, tZ, tN);
126 G4double qeRat = ratios.first*ratios.second;
129 if(qeRat<1. && proj4M.
e()>thresh) difRat /= 1.-qeRat;
140 G4cout<<
">QE>G4QFragmentation::Construct:TryQE,R="<<ratios.first*ratios.second<<
",N4M="
141 <<pN4M<<
",NPDG="<<pNPDG<<
G4endl;
143 std::pair<G4LorentzVector,G4LorentzVector> QEout=theQuasiElastic->
Scatter(pNPDG,pN4M,
147 if( (pNPDG==2212 && QEout.first.e()-mProt < CB) ||
148 (pPDG==2212 && QEout.second.e()-mProt < CB) ) CoulB =
false;
149 if(QEout.first.e() > 0 && CoulB)
152 theResult->push_back(qfN);
154 theResult->push_back(qfP);
159 theResult->push_back(resNuc);
161 G4cout<<
"-->QE-->G4QFragmentation::Construct:QEDone, N4M="<<QEout.first<<
", p4M="
167 else if(rnd < difRat)
170 G4cout<<
"-->Dif-->G4QFragmentation::Construct: qe="<<qeRat<<
", dif="<<difRat-qeRat
171 <<
",P="<<proj4M.
vect().
mag()<<
", tZ="<<tZ<<
", tN="<<tN<<
G4endl;
178 G4int nSec=out->size();
179 if(nSec>1)
for(
G4int i=0; i<nSec; i++) theResult->push_back((*out)[i]);
183 G4cout<<
"-Warning-G4QFragmentation::Construct: 1 secondary in Diffractionp 4M="
184 <<proj4M<<
", s4M="<<(*out)[0]->Get4Momentum()<<
G4endl;
195 G4cout<<
"-EMC-G4QFragmentation::Construct: Nuc4M="<<sum1<<
G4endl;
201 if(projAbsB>1) nCons = projAbsB;
204 G4double pz_per_projectile = proj4M.
pz()/nCons;
206 G4double e_per_projectile = proj4M.
e()/nCons+mProt;
207 G4double vz = pz_per_projectile/e_per_projectile;
209 G4cout<<
"G4QFragmentation::Construct: Projectile4M="<<proj4M<<
", Vz="<<vz<<
", nC="
210 <<nCons<<
", pE="<<e_per_projectile<<
G4endl;
212 theCurrentVelocity.
setZ(vz);
216 G4cout<<
"-EMC-G4QFragmentation::Construct: AfterBoost, v="<<vz<<
", Nuc4M="<<sum2<<
G4endl;
219 cmProjMom.
boost(-theCurrentVelocity);
225 if(kE > 720.) maxCt=
static_cast<int>(std::log(kE/270.));
228 G4int maxCuts=std::min( 7 , std::max(1, maxCt) );
230 G4cout<<
"G4QFragmentation::Construct: Proj4MInCM="<<cmProjMom<<
", pPDG="<<pPDG<<
G4endl;
240 G4cout<<
"G4QFrag::Constr:OutR="<<outerRadius<<
",mC="<<maxCuts<<
",A="<<theNucleus<<
G4endl;
246 G4double s_value = (cmProjMom + pNuc4M).mag2();
249 G4cout<<
"G4QFrag::Construc: p4M="<<cmProjMom<<
", tgN4M="<<pNuc4M<<
", s="<<s_value<<
", ThreshM="
255 G4cerr<<
"***G4QFragmentation::Construct: s="<<s_value<<
", pN4M="<<pNuc4M<<
G4endl;
258 else if(s_value < mProt2)
268 G4double cs=(cmProjMom + cp4M).mag2();
276 pProjectile =cmProjectile;
297 Q4M.
boost(theCurrentVelocity);
300 theQuasmons.push_back(stringQuasmon);
305 if (s_value <
sqr(ThresholdMass))
308 G4cout<<
"G4QFragmentation::Construct:*OnlyDiffraction*ThM="<<ThresholdMass<<
">sqrt(s)="
309 <<std::sqrt(s_value)<<
" -> only Diffraction is possible"<<
G4endl;
311 ModelMode = DIFFRACTIVE;
315 G4cout<<
"G4QFragmentation::Construct: theIntSize="<<theInteractions.size()<<
G4endl;
318 theInteractions.clear();
325 if (proB>0) prEn-=aProjectile.
GetMass();
326 else if(proB<0) prEn+=mProt;
329 G4cout<<
"G4QFragmentation::Construct: estimated energy, prEn="<<prEn<<
G4endl;
331 while(!theInteractions.size() && ++attCnt < maxAtt)
334 G4cout<<
"G4QFragmentation::Construct: *EnterTheInteractionLOOP*, att#"<<attCnt<<
G4endl;
337 std::pair<G4double, G4double> theImpactParameter;
339 G4double impactX = theImpactParameter.first;
340 G4double impactY = theImpactParameter.second;
342 G4cout<<
"G4QFragmentation::Construct: Impact Par X="<<impactX<<
", Y="<<impactY<<
G4endl;
344 G4double impar=std::sqrt(impactX*impactX+impactY*impactY);
347 maxEn=eflen*stringTension;
348 maxNuc=
static_cast<int>(eflen*tubeDensity+0.5);
350 G4cout<<
"G4QFragment::Construct: pE="<<prEn<<
" <? mE="<<maxEn<<
", mN="<<maxNuc<<
G4endl;
360 theInteractions.push_back(anInteraction);
365 G4cout<<
"G4QFragmentation::Construct: DINR interaction is created"<<
G4endl;
371 G4int nucleonCount = 0;
373 G4cout<<
"G4QFragment::Construct:BeforeWhileOveNuc, A="<<nA<<
",p4M="<<cmProjMom<<
G4endl;
375 while( (pNucleon=theNucleus.
GetNextNucleon()) && nucleonCount<nA && totalCuts<maxCuts)
379 s_value = (cmProjMom + pNucleon->
Get4Momentum()).mag2();
381 G4cout<<
"G4QFrag::Constr:N# "<<nucleonCount<<
", s="<<s_value<<
", tgN4M="
387 G4cout<<
"G4QFragmentation::Construct: SKIP, s<.01 GeV^2, p4M="<<cmProjMom
393 G4cout<<
"G4QFragmentation::Construct:LOOPovNuc,nC="<<nucleonCount<<
", s="<<s_value<<
G4endl;
399 G4cout<<
"G4QFragmentation::Construct: s="<<s_value<<
", D2="<<Distance2<<
G4endl;
404 G4cout<<
"G4QFragmentation::Construct: Probubility="<<Probability<<
G4endl;
409 G4cout<<
"G4QFragmentation::Construct: NLOOP prob="<<Probability<<
", rndm="<<rndNumber
410 <<
", d="<<std::sqrt(Distance2)<<
G4endl;
412 if (Probability > rndNumber)
416 G4cout<<
"--->EMC-->G4QFragmentation::Construct: Target Nucleon is filled, 4M/PDG="
424 G4UniformRand() && ModelMode==SOFT ) || ModelMode==DIFFRACTIVE)
432 theInteractions.push_back(anInteraction);
436 G4cout<<
"G4QFragmentation::Construct:NLOOP DiffInteract, tC="<<totalCuts<<
G4endl;
445 for(nCut = 0; nCut < nCutMax; nCut++)
448 if(nCut) running[nCut] += running[nCut-1];
451 for(nCut = 0; nCut < nCutMax; nCut++)
if(running[nCut] > random)
break;
454 G4cout<<
"G4QFragmentation::Construct: NLOOP-Soft Chosen nCut="<<nCut<<
G4endl;
464 theInteractions.push_back(anInteraction);
467 G4cout<<
"G4QFragmentation::Construct:NLOOP SoftInteract, tC="<<totalCuts<<
G4endl;
468 impactUsed=Distance2;
475 G4cout<<
"G4QFragmentation::Construct: NUCLEONCOUNT="<<nucleonCount<<
G4endl;
478 G4int nInt=theInteractions.size();
480 G4cout<<
"G4QFrag::Con:CUT="<<totalCuts<<
",ImpPar="<<impactUsed<<
",#ofInt="<<nInt<<
G4endl;
483 if(!nInt || (nInt==1 && theInteractions[0]->GetNumberOfDINRCollisions()==1))
489 aTarget=theInteractions[0]->GetTarget();
490 pProjectile=theInteractions[0]->GetProjectile();
491 delete theInteractions[0];
492 theInteractions.clear();
499 pProjectile =cmProjectile;
521 Q4M.
boost(theCurrentVelocity);
524 theQuasmons.push_back(stringQuasmon);
533 G4cout<<
"G4QFragmentation::Construct: Before PartPairCreation nInt="<<nInt<<
G4endl;
535 for(
G4int i=0; i<nInt; i++)
537 theInteractions[i]->SplitHadrons();
539 G4QHadron* projH=theInteractions[i]->GetProjectile();
540 G4QHadron* targH=theInteractions[i]->GetTarget();
543 std::list<G4QParton*> projCP=projH->
GetColor();
545 std::list<G4QParton*> targCP=targH->
GetColor();
547 std::list<G4QParton*>::iterator picp = projCP.begin();
548 std::list<G4QParton*>::iterator pecp = projCP.end();
549 std::list<G4QParton*>::iterator piac = projAC.begin();
550 std::list<G4QParton*>::iterator peac = projAC.end();
551 std::list<G4QParton*>::iterator ticp = targCP.begin();
552 std::list<G4QParton*>::iterator tecp = targCP.end();
553 std::list<G4QParton*>::iterator tiac = targAC.begin();
554 std::list<G4QParton*>::iterator teac = targAC.end();
555 for(; picp!=pecp&& piac!=peac&& ticp!=tecp&& tiac!=teac; ++picp,++piac,++ticp,++tiac)
557 pSP+=(*picp)->Get4Momentum();
558 pSP+=(*piac)->Get4Momentum();
559 tSP+=(*ticp)->Get4Momentum();
560 tSP+=(*tiac)->Get4Momentum();
562 G4cout<<
"-EMC-G4QFragmentation::Construct: Interaction#"<<i<<
",dP4M="
569 G4QInteractionVector::iterator it;
571 G4cout<<
"G4QFragmentation::Construct: Creation ofSoftCollisionPartonPair STARTS"<<
G4endl;
574 while(rep && theInteractions.size())
576 for(it = theInteractions.begin(); it != theInteractions.end(); ++it)
582 G4cout<<
"G4QFragmentation::Construct: #0f SOFT collisions ="<<nSoftCollisions<<
G4endl;
588 for (
G4int j = 0; j < nSoftCollisions; j++)
593 thePartonPairs.push_back(aPair);
597 thePartonPairs.push_back(aPair);
599 G4cout<<
"--->G4QFragmentation::Construct: SOFT, 2 parton pairs are filled"<<
G4endl;
603 it=theInteractions.erase(it);
604 if( it != theInteractions.begin() )
609 G4cout<<
"G4QFragmentation::Construct: *** Decremented ***"<<
G4endl;
616 G4cout<<
"G4QFragmentation::Construct: *** IT Begin ***"<<
G4endl;
623 G4cout<<
"G4QFragmentation::Construct: #0fSC="<<nSoftCollisions<<
", r="<<rep<<
G4endl;
627 G4cout<<
"G4QFragmentation::Construct: *** IT While *** , r="<<rep<<
G4endl;
631 G4cout<<
"G4QFragmentation::Construct: -> Parton pairs for SOFT strings are made"<<
G4endl;
636 for(
unsigned i = 0; i < theInteractions.size(); i++)
642 G4cout<<
"G4QFragmentation::Construct: CreationOfDiffractivePartonPairs, i="<<i<<
G4endl;
652 thePartonPairs.push_back(aPartonPair);
654 G4cout<<
"G4QFragmentation::Construct: proj Diffractive PartonPair is filled"<<
G4endl;
664 thePartonPairs.push_back(aPartonPair);
666 G4cout<<
"G4QFragmentation::Constr: target Diffractive PartonPair is filled"<<
G4endl;
671 G4cout<<
"G4QFragmentation::Construct: DiffractivePartonPairs are created"<<
G4endl;
677 theInteractions.clear();
680 G4cout<<
"G4QFragmentation::Construct: Temporary objects are cleaned up"<<
G4endl;
686 G4cout<<
"--------->>G4QFragmentation::Construct: ------->> Strings are created "<<
G4endl;
690 while(thePartonPairs.size())
692 aPair = thePartonPairs.back();
693 thePartonPairs.pop_back();
701 aString->
Boost(theCurrentVelocity);
702 strings.push_back(aString);
710 G4int nStrings=strings.size();
711 G4cout<<
"-EMC-G4QFragmentation::Construct:#ofString="<<nStrings<<
",tNuc4M="<<sum<<
G4endl;
712 for(
G4int i=0; i<nStrings; i++)
725 G4cout<<
"-EMC-G4QFragmentation::Construct: String#"<<i<<
", 4M="<<strI4M<<
",LPDG="<<LPDG
726 <<LQC<<
",RPDG="<<RPDG<<RQC<<
", Ch="<<sChg<<
", BN="<<sBaN<<
G4endl;
728 G4cout<<
"-EMC-G4QFragm::Constr: r4M="<<sum-totZLS4M<<
",rC="<<rChg<<
",rB="<<rBaN<<
G4endl;
732 G4cerr<<
"******G4QFragmentation::Construct:***** No strings are created *****"<<
G4endl;
736 G4cout<<
"G4QFragmentation::Constr: BeforeRotation, #0fStrings="<<strings.size()<<
G4endl;
741 for(
unsigned astring=0; astring < strings.size(); astring++)
742 strings[astring]->LorentzRotate(toLab);
749 G4int nStrs=strings.size();
750 G4cout<<
"-EMCLS-G4QFragmentation::Constr: #ofS="<<nStrings<<
",tNuc4M(E=M)="<<sum<<
G4endl;
751 for(
G4int i=0; i<nStrs; i++)
755 G4int sChg=strings[i]->GetCharge();
757 G4int sBaN=strings[i]->GetBaryonNumber();
759 G4cout<<
"-EMCLS-G4QFragm::Construct:String#"<<i<<
",4M="<<strI4M<<strI4M.
m()<<
",Charge="
760 <<sChg<<
",BaryN="<<sBaN<<
G4endl;
762 G4cout<<
"-EMCLS-...G4QFragm::Constr:r4M="<<sm-totLS4M<<
",rC="<<rCg<<
",rB="<<rBC<<
G4endl;
770 rCg=totChg-theNucleus.
GetZ();
771 rBC=totBaN-theNucleus.
GetA();
772 nStrs=strings.size();
773 G4cout<<
"-EMCLS-G4QFrag::Constr:AfterSwap #ofS="<<nStrings<<
",tNuc4M(E=M)="<<sum<<
G4endl;
774 for(
G4int i=0; i<nStrs; i++)
778 G4int sChg=strings[i]->GetCharge();
780 G4int sBaN=strings[i]->GetBaryonNumber();
782 G4cout<<
"-EMCLS-G4QFragm::Construct:String#"<<i<<
",4M="<<strI4M<<strI4M.
m()<<
",Charge="
783 <<sChg<<
",BaryN="<<sBaN<<
G4endl;
785 G4cout<<
"-EMCLS-...G4QFragm::Constr:r4M="<<sm-totLS4M<<
",rC="<<rCg<<
",rB="<<rBC<<
G4endl;
791 G4QStringVector::iterator ist;
793 while(con && strings.size())
795 for(ist = strings.begin(); ist < strings.end(); ++ist)
800 G4QParton* cLeft=(*ist)->GetLeftParton();
801 G4QParton* cRight=(*ist)->GetRightParton();
808 if (cLPDG > 7) aLPDG= cLPDG/100;
809 else if(cLPDG <-7) aLPDG=(-cLPDG)/100;
810 if (cRPDG > 7) aRPDG= cRPDG/100;
811 else if(cRPDG <-7) aRPDG=(-cRPDG)/100;
827 if(cSM2>0.) cSM=std::sqrt(cSM2);
829 G4cout<<
"G4QFrag::Constr:NeedsFusion? cLPDG="<<cLPDG<<
",cRPDG="<<cRPDG<<
",cM(cM2If<0)="
830 <<cSM<<
",c4M"<<cS4M<<
G4endl;
838 if(L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2)
849 G4cout<<
"G4QFrag::Const:*IsItGood? realM="<<std::sqrt(cSM2)<<
" > GSM="<<miM<<
G4endl;
851 if(std::sqrt(cSM2) > miM) bad=
false;
856 G4cout<<
"G4QFrag::Const:TryFuse,L1="<<L1<<
",L2="<<L2<<
",R1="<<R1<<
",R2="<<R2<<
G4endl;
861 G4QStringVector::iterator sst;
862 G4QStringVector::iterator pst;
867 if(cSM2>0.) minC=
false;
868 for(pst = strings.begin(); pst < strings.end(); pst++)
if(pst != ist)
877 G4cout<<
"->G4QFragm::Construct: sum4M="<<pS4M<<
",M2="<<pSM2<<
",p4M="<<sS4M<<
G4endl;
881 G4QParton* pLeft=(*pst)->GetLeftParton();
882 G4QParton* pRight=(*pst)->GetRightParton();
889 if (pLPDG > 7) LPDG= pLPDG/100;
890 else if(pLPDG <-7) LPDG=(-pLPDG)/100;
891 if (pRPDG > 7) RPDG= pRPDG/100;
892 else if(pRPDG <-7) RPDG=(-pRPDG)/100;
909 G4cout<<
"G4QFragm::Construct: Partner/w pLPDG="<<pLPDG<<
", pRPDG="<<pRPDG<<
", pM2="
915 if (cST==2 && pST==2)
920 else if(cST==2 && pST==3)
924 ( (cLPDG<0 && (-cLPDG==pL1 || -cLPDG==pL2 || -cLPDG==pRPDG) ) ||
925 (cRPDG<0 && (-cRPDG==pL1 || -cRPDG==pL2 || -cRPDG==pRPDG) )
929 ( (cLPDG<0 && (-cLPDG==pR1 || -cLPDG==pR2 || -cLPDG==pLPDG) ) ||
930 (cRPDG<0 && (-cRPDG==pR1 || -cRPDG==pR2 || -cRPDG==pLPDG) )
934 ( (cLPDG>0 && ( cLPDG==pL1 || cLPDG==pL2 || cLPDG==-pRPDG) ) ||
935 (cRPDG>0 && ( cRPDG==pL1 || cRPDG==pL2 || cRPDG==-pRPDG) )
939 ( (cLPDG>0 && ( cLPDG==pR1 || cLPDG==pR2 || cLPDG==-pLPDG) ) ||
940 (cRPDG>0 && ( cRPDG==pR1 || cRPDG==pR2 || cRPDG==-pLPDG) )
944 else G4cout<<
"G4QFragmentation::Construct:2(QaQ+QDiQ/aQaDiQ) Can't fuse"<<
G4endl;
947 else if(cST==3 && pST==2)
951 ( (pLPDG<0 && (-pLPDG==L1 || -pLPDG==L2 || -pLPDG==cRPDG) ) ||
952 (pRPDG<0 && (-pRPDG==L1 || -pRPDG==L2 || -pRPDG==cRPDG) )
956 ( (pLPDG<0 && (-pLPDG==R1 || -pLPDG==R2 || -pLPDG==cLPDG) ) ||
957 (pRPDG<0 && (-pRPDG==R1 || -pRPDG==R2 || -pRPDG==cLPDG) )
961 ( (pLPDG>0 && ( pLPDG==L1 || pLPDG==L2 || pLPDG==-cRPDG) ) ||
962 (pRPDG>0 && ( pRPDG==L1 || pRPDG==L2 || pRPDG==-cRPDG) )
966 ( (pLPDG>0 && ( pLPDG==R1 || pLPDG==R2 || pLPDG==-cLPDG) ) ||
967 (pRPDG>0 && ( pRPDG==R1 || pRPDG==R2 || pRPDG==-cLPDG) )
971 else G4cout<<
"G4QFragmentation::Construct:3(QDiQ/aQaDiQ+QaQ) Can't fuse"<<
G4endl;
974 else if(cST==2 && pST==4)
979 if ( (-cLPDG==pL1 || -cLPDG==pL2) && (cRPDG==pR1 || cRPDG==pR2) ) af=1;
980 else if( (-cRPDG==pL1 || -cRPDG==pL2) && (cLPDG==pR1 || cLPDG==pR2) ) af=2;
984 if ( (-cRPDG==pR1 || -cRPDG==pR2) && (cLPDG==pL1 || cLPDG==pL2) ) af=3;
985 else if( (-cLPDG==pR1 || -cLPDG==pR2) && (cRPDG==pL1 || cRPDG==pL2) ) af=4;
988 else G4cout<<
"-G4QFragmentation::Construct: 4 (QaQ+aQDiQDiQ) Can't fuse"<<
G4endl;
991 else if(cST==4 && pST==2)
996 if ( (-pLPDG==L1 || -pLPDG==L2) && (pRPDG==R1 || pRPDG==R2) ) af=1;
997 else if( (-pRPDG==L1 || -pRPDG==L2) && (pLPDG==R1 || pLPDG==R2) ) af=2;
1001 if ( (-pRPDG==R1 || -pRPDG==R2) && (pLPDG==L1 || pLPDG==L2) ) af=3;
1002 else if( (-pLPDG==R1 || -pLPDG==R2) && (pRPDG==L1 || pRPDG==L2) ) af=4;
1005 else G4cout<<
"-G4QFragmentation::Construct: 5 (aQDiQDiQ+QaQ) Can't fuse"<<
G4endl;
1008 else if(cST==3 && pST==3)
1013 if (cLPDG<-7 && (-cRPDG==pL1 || -cRPDG==pL2) && (pRPDG==L1 || pRPDG==L2))
1015 else if(cRPDG<-7 && (-cLPDG==pL1 || -cLPDG==pL2) && (pRPDG==R1 || pRPDG==R2))
1020 if (cLPDG<-7 && (-cRPDG==pR1 || -cRPDG==pR2) && (pLPDG==L1 || pLPDG==L2))
1022 else if(cRPDG<-7 && (-cLPDG==pR1 || -cLPDG==pR2) && (pLPDG==R1 || pLPDG==R2))
1027 if (pLPDG<-7 && (-pRPDG==L1 || -pRPDG==L2) && (cRPDG==pL1 || cRPDG==pL2))
1029 else if(pRPDG<-7 && (-pLPDG==L1 || -pLPDG==L2) && (cRPDG==pR1 || cRPDG==pR2))
1034 if (pLPDG<-7 && (-pRPDG==R1 || -pRPDG==R2) && (cLPDG==pL1 || cLPDG==pL2))
1036 else if(pRPDG<-7 && (-pLPDG==R1 || -pLPDG==R2) && (cLPDG==pR1 || cLPDG==pR2))
1040 else G4cout<<
"-G4QFragmentation::Construct: 6 (QDiQ+aQaDiQ) Can't fuse"<<
G4endl;
1044 G4cout<<
"G4QFragmentation::Const: ***Possibility***, tf="<<tf<<
", af="<<af<<
G4endl;
1053 if (cLPDG > 0 && pLPDG > 0)
1056 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
1057 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
1058 else nLPDG=pLPDG*1000+cLPDG*100+3;
1059 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
1060 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
1061 else nRPDG=pRPDG*1000+cRPDG*100-3;
1063 else if(cLPDG < 0 && pLPDG < 0)
1066 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
1067 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
1068 else nRPDG=pRPDG*1000+cRPDG*100+3;
1069 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
1070 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
1071 else nLPDG=pLPDG*1000+cLPDG*100-3;
1073 else if(cRPDG > 0 && pLPDG > 0)
1076 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
1077 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
1078 else nLPDG=pLPDG*1000+cRPDG*100+3;
1079 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
1080 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
1081 else nRPDG=pRPDG*1000+cLPDG*100-3;
1083 else if(cRPDG < 0 && pLPDG < 0)
1086 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
1087 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
1088 else nRPDG=pRPDG*1000+cLPDG*100+3;
1089 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
1090 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
1091 else nLPDG=pLPDG*1000+cRPDG*100-3;
1108 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
1109 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
1110 else nRPDG=pRPDG*1000+cRPDG*100+3;
1111 if (-cLPDG == pL1) nLPDG=pL2;
1125 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
1126 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
1127 else nRPDG=pRPDG*1000+cLPDG*100+3;
1128 if (-cRPDG == pL1) nLPDG=pL2;
1144 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
1145 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
1146 else nLPDG=pLPDG*1000+cRPDG*100+3;
1147 if (-cLPDG == pR1) nRPDG=pR2;
1161 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
1162 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
1163 else nLPDG=pLPDG*1000+cLPDG*100+3;
1164 if (-cRPDG == pR1) nRPDG=pR2;
1180 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
1181 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
1182 else nRPDG=pRPDG*1000+cRPDG*100-3;
1183 if ( cLPDG == pL1) nLPDG=-pL2;
1197 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
1198 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
1199 else nRPDG=pRPDG*1000+cLPDG*100-3;
1200 if ( cRPDG == pL1) nLPDG=-pL2;
1216 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
1217 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
1218 else nLPDG=pLPDG*1000+cRPDG*100-3;
1219 if ( cLPDG == pR1) nRPDG=-pR2;
1233 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
1234 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
1235 else nLPDG=pLPDG*1000+cLPDG*100-3;
1236 if ( cRPDG == pR1) nRPDG=-pR2;
1257 if (pRPDG > cRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
1258 else if(pRPDG < cRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
1259 else nRPDG=cRPDG*1000+pRPDG*100+3;
1260 if (-pLPDG == L1) nLPDG=L2;
1274 if (pLPDG > cRPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
1275 else if(pLPDG < cRPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
1276 else nLPDG=cRPDG*1000+pLPDG*100+3;
1277 if (-pRPDG == L1) nRPDG=L2;
1293 if (pRPDG > cLPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
1294 else if(pRPDG < cLPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
1295 else nRPDG=cLPDG*1000+pRPDG*100+3;
1296 if (-pLPDG == R1) nLPDG=R2;
1310 if (pLPDG > cLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
1311 else if(pLPDG < cLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
1312 else nLPDG=cLPDG*1000+pLPDG*100+3;
1313 if (-pRPDG == R1) nRPDG=R2;
1329 if (pRPDG < cRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
1330 else if(pRPDG > cRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
1331 else nRPDG=cRPDG*1000+pRPDG*100-3;
1332 if ( pLPDG == L1) nLPDG=-L2;
1346 if (pLPDG < cRPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
1347 else if(pLPDG > cRPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
1348 else nLPDG=cRPDG*1000+pLPDG*100-3;
1349 if ( pRPDG == L1) nRPDG=-L2;
1365 if (pRPDG < cLPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
1366 else if(pRPDG > cLPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
1367 else nRPDG=cLPDG*1000+pRPDG*100-3;
1368 if ( pLPDG == R1) nLPDG=-R2;
1382 if (pLPDG < cLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
1383 else if(pLPDG > cLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
1384 else nLPDG=cLPDG*1000+pLPDG*100-3;
1385 if ( pRPDG == R1) nRPDG=-R2;
1397 if(-cLPDG == pL1) nLPDG= pL2;
1399 if( cRPDG == pR1) nRPDG=-pR2;
1404 if(-cRPDG == pL1) nLPDG= pL2;
1406 if( cLPDG == pR1) nRPDG=-pR2;
1411 if( cLPDG == pL1) nLPDG=-pL2;
1413 if(-cRPDG == pR1) nRPDG= pR2;
1418 if( cRPDG == pL1) nLPDG=-pL2;
1420 if(-cLPDG == pR1) nRPDG= pR2;
1430 if(-pLPDG == L1) nLPDG= L2;
1432 if( pRPDG == R1) nRPDG=-R2;
1437 if(-pRPDG == L1) nRPDG= L2;
1439 if( pLPDG == R1) nLPDG=-R2;
1444 if( pLPDG == L1) nLPDG=-L2;
1446 if(-pRPDG == R1) nRPDG= R2;
1451 if( pRPDG == L1) nRPDG=-L2;
1453 if(-pLPDG == R1) nLPDG= R2;
1463 if(-cRPDG == pL1) nLPDG= pL2;
1465 if( pRPDG == L1) nRPDG= -L2;
1470 if(-cLPDG == pL1) nLPDG= pL2;
1472 if( pRPDG == R1) nRPDG= -R2;
1477 if(-cRPDG == pR1) nRPDG= pR2;
1479 if( pLPDG == L1) nLPDG= -L2;
1484 if(-cLPDG == pR1) nRPDG= pR2;
1486 if( pLPDG == R1) nLPDG= -R2;
1491 if(-pRPDG == L1) nRPDG= L2;
1493 if( cRPDG == pL1) nLPDG=-pL2;
1498 if(-pLPDG == L1) nLPDG= L2;
1500 if( cRPDG == pR1) nRPDG=-pR2;
1505 if(-pRPDG == R1) nRPDG= R2;
1507 if( cLPDG == pL1) nLPDG=-pL2;
1512 if(-pLPDG == R1) nLPDG= R2;
1514 if( cLPDG == pR1) nRPDG=-pR2;
1520 if(!order)
G4cerr<<
"-Warning-G4QFrag::Constr: t="<<tf<<
", a="<<af<<
", cL="<<cLPDG
1521 <<
", cR="<<cRPDG<<
", pL="<<pLPDG<<
", pR="<<pRPDG<<
G4endl;
1526 if(std::abs(nLPDG) > 7) ++
LT;
1528 if(std::abs(nRPDG) > 7) ++RT;
1531 if(cLT==2 && cRT==2)
1549 if(nL1!=nR1 && nL1!=nR2 && nL2!=nR1 && nL2!=nR2)
1552 G4cout<<
"G4QFragmentation::Const:aLPDG="<<aLPDG<<
", aRPDG="<<aRPDG<<
G4endl;
1556 std::pair<G4int,G4int> pB=tmp.
MakeTwoBaryons(nL1, nL2, nR1, nR2);
1562 std::pair<G4int,G4int> newPair = std::make_pair(nLPDG,nRPDG);
1565 G4cout<<
"G4QFr::Con: LPDG="<<nLPDG<<
",RPDG="<<nRPDG<<
",QC="<<newStQC<<
G4endl;
1573 if(pSM2>0.) pSM=std::sqrt(pSM2);
1574 if(minC && pSM2 > maxiM2)
1579 else if(!minC || pSM > minM)
1582 if(minC || pExcess > excess)
1604 G4QParton* pLeft=(*sst)->GetLeftParton();
1605 G4QParton* pRight=(*sst)->GetRightParton();
1609 G4cout<<
"G4QFragmentation::Const:cS4M="<<cS4M<<
" fused/w pS4M="<<pL4M+pR4M<<
G4endl;
1629 G4cout<<
"G4QFragmentation::Construct:Created,4M="<<ss4M<<
",m2="<<ss4M.
m2()<<
G4endl;
1631 if( ist != strings.begin() )
1636 G4cout<<
"G4QFragmentation::Construct: *** IST Decremented ***"<<
G4endl;
1643 G4cout<<
"G4QFragmentation::Construct: *** IST Begin ***"<<
G4endl;
1651 G4cout<<
"-Warning-G4QFragm::Const:S4M="<<cS4M<<
",M2="<<cSM2<<
" Leave ASIS"<<
G4endl;
1660 G4cout<<
"G4QFragmentation::Construct: *** IST While *** , con="<<con<<
G4endl;
1668 G4int nStr=strings.size();
1669 G4cout<<
"-EMCLS-G4QFr::Const: AfterSUPPRESION #ofS="<<nStr<<
",tNuc4M(E=M)="<<sum<<
G4endl;
1670 for(
G4int i=0; i<nStr; i++)
1674 G4int sChg=strings[i]->GetCharge();
1676 G4int sBaN=strings[i]->GetBaryonNumber();
1678 G4cout<<
"-EMCLS-G4QFragm::Construct: St#"<<i<<
", 4M="<<strI4M<<
", M="<<strI4M.
m()
1679 <<
", C="<<sChg<<
", B="<<sBaN<<
G4endl;
1681 G4cout<<
"-EMCLS-G4QFragm::Construct:r4M="<<t4M-totLS4M<<
",rC="<<rC<<
",rB="<<rB<<
G4endl;
1687 G4cout<<
"G4QFragmentation::Construct: problem="<<problem<<
G4endl;
1691 G4int nOfStr=strings.size();
1693 G4cout<<
"G4QFragmentation::Construct:SecurityDiQaDiQReduction,#OfStr="<<nOfStr<<
G4endl;
1695 for (
G4int astring=0; astring < nOfStr; astring++)
1707 G4cout<<
"G4QFragmentation::Constr:TrySelfReduString,L="<<sPDG<<
",R="<<nPDG<<
G4endl;
1714 G4cout<<
"+G4QFragm::Const:#"<<astring<<
" Reduced, L="<<sPDG<<
",R="<<nPDG<<
G4endl;
1718 else G4cout<<
"--*--G4QFragm::Const:#"<<astring<<
" DQ-aDQ reduction Failed"<<
G4endl;
1721 else if(sPDG==3 && nPDG==-3)
1728 else if(sPDG==-3 && nPDG==3)
1742 G4int nStri=strings.size();
1743 G4cout<<
"-EMCLS-G4QFr::Const: FinalConstruct, #ofSt="<<nStri<<
",tN4M(E=M)="<<t4M<<
G4endl;
1744 for(
G4int i=0; i<nStri; i++)
1748 G4int sChg=strings[i]->GetCharge();
1750 G4int sBaN=strings[i]->GetBaryonNumber();
1752 G4cout<<
"-EMCLS-G4QFragm::Construct: St#"<<i<<
", 4M="<<strI4M<<
", M="<<strI4M.
m()
1753 <<
", C="<<sChg<<
", B="<<sBaN<<
G4endl;
1755 G4cout<<
"-EMCLS-G4QFragm::Construct:r4M="<<u4M-totLS4M<<
",rC="<<rCh<<
",rB="<<rBa<<
G4endl;
1761 std::for_each(strings.begin(), strings.end(),
DeleteQString() );
1775 G4cout<<
"*******>G4QFragmentation::Fragment: ***Called***, Res="<<theResult<<
G4endl;
1777 G4int striNum=strings.size();
1778 G4int hadrNum=theResult->size();
1780 G4int nQm=theQuasmons.size();
1784 G4cout<<
"-EMCLS-G4QF::Fragment: CHECKRecovery, #ofS="<<striNum<<
", #Nuc4M(E=M)="<<totLS4M
1785 <<
",#Q="<<nQm<<
",#H="<<hadrNum<<
G4endl;
1786 for(
G4int i=0; i < striNum; i++)
1790 G4int sChg=strings[i]->GetCharge();
1792 G4int sBaN=strings[i]->GetBaryonNumber();
1794 G4cout<<
"-EMCLS-G4QFragm::Fragment: String#"<<i<<
", 4M="<<strI4M<<
", M="<<strI4M.
m()
1795 <<
", C="<<sChg<<
", B="<<sBaN<<
G4endl;
1797 for(
G4int i=0; i < nQm; i++)
1801 G4int hChg=theQuasmons[i]->GetCharge();
1803 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
1805 G4cout<<
"-EMCLS-G4QFragmentation::Fragment: Quasmon#"<<i<<
", 4M="<<hI4M<<
", C="<<hChg
1808 for(
G4int i=0; i < hadrNum; i++)
1812 G4int hChg=(*theResult)[i]->GetCharge();
1814 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
1816 G4cout<<
"-EMCLS-G4QFr::Fragment:H#"<<i<<
",4M="<<hI4M<<
",C="<<hChg<<
",B="<<hBaN<<
G4endl;
1820 G4cout<<
"***>G4QFragmentation::Fragment: #OfStr="<<striNum<<
", #OfRes="<<hadrNum<<
G4endl;
1822 if(!striNum && hadrNum)
1825 G4cout<<
"***>G4QFragmentation::Fragment:**Quasi-Elastic**,#OfResult="<<hadrNum<<
G4endl;
1834 for(
G4int ih=0; ih<hadrNum; ih++)
delete (*theResult)[ih];
1840 theResult->push_back(resNuc);
1842 G4int nQuas=theQuasmons.size();
1843 G4int theRS=theResult->size();
1845 G4cout<<
"***>G4QFragmentation::Fragment:beforeEnv,#OfQ="<<nQuas<<
",#OfR="<<theRS<<
G4endl;
1849 G4QHadron* resNuc = (*theResult)[theRS-1];
1852 if(resNucPDG==90000000 || resNuc4M.
m2()<800000.)
1855 if(resNucPDG == 90000000) resNuc->
Set4Momentum(resNuc4M);
1863 theResult->pop_back();
1866 G4cout<<
"G4QFragmentation::Fragment:#OfRemainingHadron="<<theRS<<
",A="<<theEnv<<
G4endl;
1869 for(
G4int j=theRS-1; j>-2; --j)
1874 G4cout<<
"G4QFragmentation::Fragm:rN4M"<<qsum4M<<qsum4M.
m()<<
",rNQC="<<qsumQC<<
G4endl;
1879 for(
G4int i=0; i<nQuas; ++i)
1887 G4cout<<
"G4QFr::Fr:Q#"<<i<<
",Q4M="<<cur4M<<
",QQC="<<curQC<<
",sQC="<<qsumQC<<
G4endl;
1899 G4cout<<
"G4QFragmentation::Fragment: QC="<<qsumQC<<
",PDG="<<miPDG<<
G4endl;
1914 else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2)
1919 G4cout<<
"G4QFragmentation::Fragment: PDG="<<miPDG<<
",rM="<<reM<<
",GSM="<<gsM<<
G4endl;
1921 if(reM > gsM)
break;
1928 firstQ->
SetQC(firstQC+hQC);
1930 theResult->pop_back();
1937 G4cerr<<
"*G4QFr::Fr:PDG="<<miPDG<<
",M="<<reM<<
",GSM="<<gsM<<
",QC="<<qsumQC<<
G4endl;
1942 if(nucE < 1.E-12) nucE=0.;
1943 else if(resNucPDG==22 && nQuas==1)
1953 G4cout<<
"G4QFragm::Fragm: nucE="<<nucE<<
",nQ="<<nQuas<<
G4endl;
1955 if(nucE) nucVel=resNuc4M.
vect()/nucE;
1958 G4int sqCg=rnChg-totChg;
1959 G4int sqBN=rnBaN-totBaN;
1961 for(
G4int i=0; i<nQuas; ++i)
1966 <<resNuc4M<<resNucPDG<<
G4endl;
1968 if(nucE) curQuasm->
Boost(-nucVel);
1972 G4cout<<
"G4QFragmentation::Fragment: Quasmon# "<<i<<
" added, 4M="<<cQ4M<<
G4endl;
1979 G4cout<<
"-EMCLS-G4QFrag::Fragm: r4M="<<sq4M<<
", rC="<<sqCg<<
", rB="<<sqBN<<
G4endl;
1984 G4cout<<
"G4QFrag::Fragm: *** Before Del Output ***"<<
G4endl;
1988 G4cout<<
"G4QFrag::Fragm: *** After Del Output ***"<<
G4endl;
1994 G4cerr<<
"***G4QFragmentation::Fragment: G4QE Exception is catched"<<
G4endl;
1996 G4Exception(
"G4QFragmentation::Fragment()",
"HAD_CHPS_0027",
2000 G4cout<<
"G4QFrag::Fragm: *** Before Del Pan ***"<<
G4endl;
2004 G4cout<<
"G4QFrag::Fragm: *** After Del Pan ***"<<
G4endl;
2008 G4int nOut=output->size();
2009 for(
G4int j=0; j<nOut; j++)
2012 if(nucE) curHadr->
Boost(nucVel);
2015 if((hPDG!=90000000 && hPDG!=22) || h4M!=nul4M) theResult->push_back(curHadr);
2016 else delete curHadr;
2021 else if(!striNum)
G4cout<<
"-Warning-G4QFragmentation::Fragment:Nothing was done"<<
G4endl;
2023 G4cout<<
"=--=>G4QFragmentation::Fragment: Final #OfResult="<<theResult->size()<<
G4endl;
2025 G4int nQ =theQuasmons.size();
2026 if(nQ) theQuasmons.clear();
2027 G4int nHd=theResult->size();
2032 G4cout<<
"-EMCLS-G4QFragmentation::Fragment: #ofHadr="<<nHd<<
", #OfQuasm="<<nQ<<
G4endl;
2033 for(
G4int i=0; i<nHd; i++)
2037 G4int hChg=(*theResult)[i]->GetCharge();
2039 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
2041 G4cout<<
"-EMCLS-G4QFragmentation::Fragment: Hadron#"<<i<<
", 4M="<<hI4M<<
", PDG="
2042 <<(*theResult)[i]->GetPDGCode()<<
", C="<<hChg<<
", B="<<hBaN<<
G4endl;
2044 G4cout<<
"-EMCLS-G4QFrag::Fragm: r4M="<<f4M-totLS4M<<
", rC="<<fCh<<
", rB="<<fBN<<
G4endl;
2047 G4QHadron* resNuc = (*theResult)[nHd-1];
2050 if(rnBn==1 && (rnCg==-2 || rnCg==3 || rnCg==-1 || rnCg==2))
2064 if(rnCg==-2 || rnCg==3)
2067 if(!
G4QHadron(tot4M).DecayIn3(m14M,m24M,n4M))
2069 G4cerr<<
"***G4QFrag::Frag: tM="<<tot4M.
m()<<
" -> m1="<<mPiCh<<
" + m2="<<mPiCh
2070 <<
" + nM="<<nM<<
" = "<<2*mPiCh+nM<<
G4endl;
2073 theResult->pop_back();
2076 theResult->push_back(m1H);
2078 G4cout<<
"G4QFragment::Fragment:DecayIn3, M1="<<mPDG<<m14M<<
G4endl;
2081 theResult->push_back(m2H);
2083 G4cout<<
"G4QFragment::Fragment:DecayIn3, M2="<<mPDG<<m24M<<
G4endl;
2086 theResult->push_back(nH);
2088 G4cout<<
"G4QFragment::Fragment:DecayIn3, Nucleon="<<nPDG<<n4M<<
G4endl;
2093 if(!
G4QHadron(tot4M).DecayIn2(m14M,n4M))
2095 G4cerr<<
"***G4QFrag::Frag: tM="<<tot4M.
m()<<
" -> m1="<<mPiCh
2096 <<
" + nM="<<nM<<
" = "<<mPiCh+nM<<
G4endl;
2099 theResult->pop_back();
2102 theResult->push_back(m1H);
2104 G4cout<<
"G4QFragment::Fragment:DecayIn2, M1="<<mPDG<<m14M<<
G4endl;
2107 theResult->push_back(nH);
2109 G4cout<<
"G4QFragment::Fragment:DecayIn2, Nucleon="<<nPDG<<n4M<<
G4endl;
2120 if(!
G4QHadron(tot4M).DecayIn2(n14M,n24M))
2122 G4cerr<<
"***G4QFrag::Frag: tM="<<tot4M.
m()<<
" -> n*2="<<2*mNeut<<
G4endl;
2125 theResult->pop_back();
2128 theResult->push_back(n1H);
2130 G4cout<<
"G4QFragment::Fragment:DecayIn2, Neutron1="<<n14M<<
G4endl;
2133 theResult->push_back(n2H);
2135 G4cout<<
"G4QFragment::Fragment:DecayIn2, Neutron2="<<n24M<<
G4endl;
2143 if(!
G4QHadron(tot4M).DecayIn2(n14M,n24M))
2145 G4cerr<<
"***G4QFrag::Frag: tM="<<tot4M.
m()<<
" -> n*2="<<2*mProt<<
G4endl;
2148 theResult->pop_back();
2151 theResult->push_back(n1H);
2153 G4cout<<
"G4QFragment::Fragment:DecayIn2, Proton1="<<n14M<<
G4endl;
2156 theResult->push_back(n2H);
2158 G4cout<<
"G4QFragment::Fragment:DecayIn2, Proton2="<<n24M<<
G4endl;
2163 nHd=theResult->size();
2169 for(
G4int i=0; i<nHd; ++i)
2174 if(hPDG==90000000 || hPDG==22)
2179 else if( curh4M == nul4M)
2181 G4QHadron* theLast = (*theResult)[nHd-1];
2182 if(theLast != curHadr)
2189 theResult->pop_back();
2192 if(i == nHd-1)
break;
2198 for(
G4int j=i+1; j<nHd; ++j)
2207 if(hPDG==2212) E -= mProt+mProt;
2208 else E -= mNeut+mNeut;
2212 if(hPDG==2212) piPDG=-211;
2213 for(
G4int k=0; k<nHd; ++k)
2215 G4int mPDG=(*theResult)[k]->GetPDGCode();
2218 if(mPDG==111 || mPDG==piPDG)
2244 if(S > D2) found= 1;
2250 if(S > D2) found=-1;
2263 if(p2 < 0.)
G4cout<<
"-Warning-G4QFragment::Fragment: P2="<<p2<<
G4endl;
2267 if(pc2 < .00000000000001)
2268 G4cout<<
"-Warning-G4QFragment::Fragment: PC2="<<pc2<<m4M<<
G4endl;
2269 else r=std::sqrt(p2/pc2);
2271 m4M.
setE(std::sqrt(mPi2+p2));
2275 n4M.
setE(std::sqrt(mN2+p2));
2278 (*theResult)[k]->SetPDGCode(sPDG);
2279 (*theResult)[k]->Set4Momentum(m4M);
2282 (*theResult)[i]->SetPDGCode(nPDG);
2283 (*theResult)[i]->Set4Momentum(n4M);
2287 (*theResult)[j]->SetPDGCode(nPDG);
2288 (*theResult)[j]->Set4Momentum(n4M);
2304 maxQC = cHadr->
GetQC();
2312 G4cout<<
"G4QFragmentation::Fra: ResNucEnv with A="<<maxBN<<
", Z="<<maxChg<<
G4endl;
2318 G4int nHadr=theResult->size();
2320 if(nHadr>2)
for(
unsigned f=0; f<theResult->size(); f++)
2322 G4int fBN=(*theResult)[f]->GetBaryonNumber();
2324 G4int fPDG=(*theResult)[f]->GetPDGCode();
2326 G4cout<<
"G4QFragmentation::Fra:"<<f<<
",PDG="<<fPDG<<
",fBN="<<fBN<<
",f4M="<<fLV<<
G4endl;
2335 G4cout<<
"G4QFrag::Frag:=>Before Gamma Suppression<=, nH="<<nHadr<<
",frag="<<frag<<
G4endl;
2337 if(nHadr>2 && frag)
for(
G4int h=nHadr-1; h>=0; h--)
2343 if(hPDG==89999003||hPDG==90002999)
G4cout<<
"-Warning-G4QFr::Fr:nD-/pD++="<<hPDG<<
G4endl;
2345 G4cout<<
"G4QFragmentation::Fragm: h#"<<h<<
", hPDG="<<hPDG<<
", hNFrag="<<hF<<
G4endl;
2349 if(hChg > 0 && hPDG!=321)
2353 if(h4M.
e()-h4M.
m() < hCB) UCB =
true;
2355 if(hF || hPDG==22 || UCB)
2357 G4int last=theResult->size()-1;
2358 G4QHadron* theLast = (*theResult)[last];
2363 if(UCB) sumQC+=curHadr->
GetQC();
2365 G4cout<<
"G4QFragmentation::Frag: gam4M="<<h4M<<
" is added to s4M="<<sum<<
G4endl;
2368 nHadr =
static_cast<G4int>(theResult->size())-1;
2377 G4cout<<
"G4QFragmentation::Fragment: Exchange with the last is done"<<
G4endl;
2380 theResult->pop_back();
2383 G4cout<<
"G4QFragmentation::Fragment: The last is compessed"<<
G4endl;
2388 G4cout<<
"G4QFragment::Frag: nH="<<nHadr<<
"="<<theResult->size()<<
", sum="<<sum<<
G4endl;
2390 if(nHadr > 1)
for(
unsigned hdr=0; hdr<theResult->size()-1; hdr++)
2396 G4QHadron* theLast = (*theResult)[theResult->size()-1];
2414 nHadr=theResult->size();
2417 G4QHadron* theLast = (*theResult)[nHadr-1];
2425 theResult->pop_back();
2431 G4cout<<
"G4QFra::Fra:gSum4M="<<sum<<
" is added to "<<new4M<<
", QC="<<newQC<<
G4endl;
2436 theNew->
SetQC(exResQC);
2441 nHadr=theResult->size();
2443 else G4cout<<
"-Warning-G4QFragmentation::Fragment:RA="<<nucEnvBN<<
",E/M cons?"<<
G4endl;
2461 G4int nStri=strings.size();
2462 G4cout<<
"-EMCLS-G4QFr::Breed: CHECKRecovery #ofS="<<nStri<<
",N4M(E=M)="<<totLS4M<<
G4endl;
2463 for(
G4int i=0; i<nStri; i++)
2467 G4int sChg=strings[i]->GetCharge();
2469 G4int sBaN=strings[i]->GetBaryonNumber();
2471 G4cout<<
"-EMCLS-G4QFragm::Breeder: St#"<<i<<
", 4M="<<strI4M<<
", M="<<strI4M.
m()
2472 <<
", C="<<sChg<<
", B="<<sBaN<<
G4endl;
2475 G4int nOfStr=strings.size();
2477 G4cout<<
"G4QFragmentation::Breeder: BeforeFragmentation, #OfStr="<<nOfStr<<
G4endl;
2482 for(
G4int i=0; i < nOfStr; ++i)
2489 if(pS4M.
m2() < 0.) ftBad=
true;
2491 G4cout<<
">G4QFrag::Br:1stTest,S#"<<i<<
",P="<<crStr<<
",4M="<<pS4M<<
",QC="<<pSQC<<
G4endl;
2498 G4cout<<
"->G4QFragmentation::Breeder:*TotQ*,QC="<<ftQC<<
",4M="<<ft4M<<ft4M.
m()<<
G4endl;
2500 theQuasmons.push_back(stringQuasmon);
2504 theResult->push_back(resNuc);
2507 for (
G4int astring=0; astring < nOfStr; astring++)
2513 G4int nOfHadr=theResult->size();
2514 G4cout<<
"-EMCLS-G4QFragmentation::Breeder:#ofSt="<<nOfStr<<
",#ofHad="<<nOfHadr<<
G4endl;
2515 for(
G4int i=astring; i<nOfStr; i++)
2519 G4int sChg=strings[i]->GetCharge();
2521 G4int sBaN=strings[i]->GetBaryonNumber();
2523 G4cout<<
"-EMCLS-G4QF::Breed:S#"<<i<<
",4M="<<strI4M<<
",C="<<sChg<<
",B="<<sBaN<<
G4endl;
2525 for(
G4int i=0; i<nOfHadr; i++)
2529 G4int hChg=(*theResult)[i]->GetCharge();
2531 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
2533 G4cout<<
"-EMCLS-G4QFr::Breed: H#"<<i<<
",4M="<<hI4M<<
",C="<<hChg<<
",B="<<hBaN<<
G4endl;
2535 G4cout<<
"....-EMCLS-G4QFrag::Br:r4M="<<sum-totLS4M<<
",rC="<<rChg<<
",rB="<<rBaN<<
G4endl;
2543 G4cout<<
"=--=>G4QFragmentation::Breeder: String#"<<astring<<
",s4M/m="<<curString4M
2562 G4cout<<
"G4QFragmentation::Breeder:TryReduceString, L="<<sPDG<<
",R="<<nPDG<<
G4endl;
2572 G4cout<<
"G4QFragmentation::Breeder:AfterReduction,L="<<sPDG<<
",R="<<nPDG<<
G4endl;
2583 else G4cout<<
"^G4QFragmentation::Breeder: DQ-aDQ reduction to Q-aQ Failed"<<
G4endl;
2587 G4cout<<
"G4QFrag::Breed:AfterRedAttempt, theH="<<theHadrons<<
", L4M="
2590 unsigned next=astring+1;
2594 if(next < strings.size())
2601 if(dPDG<-99 || (dPDG>0&&dPDG<7) || qPDG>99 || (qPDG<0 && qPDG>-7))
2607 if(dPDG>99) dPDG/=100;
2608 if(qPDG<-99) qPDG=-(-qPDG)/100;
2610 G4cout<<
"G4QFrag::Breed:TryFuseStringS, q="<<qPDG<<
", a="<<dPDG<<
", n="<<next
2617 for (restr=next; restr < nOfStr; restr++)
2629 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
2635 if(aPDG > 99) aPDG/=100;
2636 if(rPDG <-99) rPDG=-(-rPDG)/100;
2639 G4cout<<
"G4QFragm::Breed: TryReduce #"<<restr<<
",q="<<rPDG<<
",a="<<aPDG<<
G4endl;
2641 if(
LT==2 && RT==2 && PLT==2 && PRT==2)
2643 G4int cQ1=(-qPDG)/10;
2644 G4int cQ2=(-qPDG)%10;
2647 G4int pQ1=(-rPDG)/10;
2648 G4int pQ2=(-rPDG)%10;
2652 G4cout<<
"G4QFragment::Breeder: cQ="<<cQ1<<
","<<cQ2<<
", cA="<<cA1<<
","<<cA2
2653 <<
", pQ="<<pQ1<<
","<<pQ2<<
", pA="<<pA1<<
","<<pA2<<
G4endl;
2655 G4bool iQA = (cA1==pQ1 || cA1==pQ2 || cA2==pQ1 || cA2==pQ2);
2656 G4bool iAQ = (cQ1==pA1 || cQ1==pA2 || cQ2==pA1 || cQ2==pA2);
2661 if(sPDG>0 && uPDG<0)
2663 std::pair<G4int,G4int> resLL=
ReducePair(sPDG/100, (-uPDG)/100);
2664 G4int newCL=resLL.first;
2665 G4int newPL=resLL.second;
2666 if(!newCL || !newPL)
2668 G4cerr<<
"*G4QFragmentation::Breeder:CL="<<newCL<<
",PL="<<newPL<<
G4endl;
2671 std::pair<G4int,G4int> resRR=
ReducePair((-nPDG)/100, mPDG/100);
2672 G4int newCR=resRR.first;
2673 G4int newPR=resRR.second;
2674 if(!newCR || !newPR)
2676 G4cerr<<
"*G4QFragmentation::Breeder:CR="<<newCR<<
",PR="<<newPR<<
G4endl;
2685 else if(sPDG<0 && uPDG>0)
2687 std::pair<G4int,G4int> resLL=
ReducePair((-sPDG)/100, uPDG/100);
2688 G4int newCL=resLL.first;
2689 G4int newPL=resLL.second;
2690 if(!newCL || !newPL)
2692 G4cerr<<
"*G4QFragmentation::Breeder:CL="<<newCL<<
",PL="<<newPL<<
G4endl;
2695 std::pair<G4int,G4int> resRR=
ReducePair(nPDG/100, (-mPDG)/100);
2696 G4int newCR=resRR.first;
2697 G4int newPR=resRR.second;
2698 if(!newCR || !newPR)
2700 G4cerr<<
"*G4QFragmentation::Breeder:CR="<<newCR<<
",PR="<<newPR<<
G4endl;
2709 else if(sPDG>0 && mPDG<0)
2711 std::pair<G4int,G4int> resLL=
ReducePair(sPDG/100, (-mPDG)/100);
2712 G4int newCL=resLL.first;
2713 G4int newPR=resLL.second;
2714 if(!newCL || !newPR)
2716 G4cerr<<
"*G4QFragmentation::Breeder:CL="<<newCL<<
",PR="<<newPR<<
G4endl;
2719 std::pair<G4int,G4int> resRR=
ReducePair((-nPDG)/100, uPDG/100);
2720 G4int newCR=resRR.first;
2721 G4int newPL=resRR.second;
2722 if(!newCR || !newPR)
2724 G4cerr<<
"*G4QFragmentation::Breeder:CR="<<newCR<<
",PL="<<newPL<<
G4endl;
2735 std::pair<G4int,G4int> resLL=
ReducePair((-sPDG)/100, mPDG/100);
2736 G4int newCL=resLL.first;
2737 G4int newPR=resLL.second;
2738 if(!newCL || !newPR)
2740 G4cerr<<
"*G4QFragmentation::Breeder:CL="<<newCL<<
",PR="<<newPR<<
G4endl;
2743 std::pair<G4int,G4int> resRR=
ReducePair(nPDG/100, (-uPDG)/100);
2744 G4int newCR=resRR.first;
2745 G4int newPL=resRR.second;
2746 if(!newCR || !newPR)
2748 G4cerr<<
"*G4QFragmentation::Breeder:CR="<<newCR<<
",PL="<<newPL<<
G4endl;
2762 G4cout<<
"G4QFragm::Breed:TryFuse/w #"<<restr<<
",q="<<rPDG<<
",a="<<aPDG<<
G4endl;
2765 if( (LS==2 && PLS==2) ||
2766 ( ( (LS==2 && PLS==3) || (LS==3 && PLS==2) ) &&
2767 ( (aPDG> 7 && (-dPDG==aPDG/10 || -dPDG==aPDG%10) ) ||
2768 (dPDG> 7 && (-aPDG==dPDG/10 || -aPDG==dPDG%10) ) ||
2769 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) ) ||
2770 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )
2774 ( ( LS==2 && PLS==4 &&
2775 (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) ) &&
2776 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) )
2778 ( LS==4 && PLS==2 &&
2779 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) ) &&
2780 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )
2783 ( LS==3 && PLS==3 &&
2784 ( (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) &&
2785 qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10)
2787 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) &&
2788 rPDG<-7 && (dPDG==(-rPDG)/10 || dPDG==(-rPDG)%10)
2798 G4cout<<
"G4QFragmentation::Breeder: StringCand#"<<restr<<
", q="<<rPDG<<
", a="
2799 <<aPDG<<
", L="<<uPDG<<
", R="<<mPDG<<
",dV="<<dV<<
G4endl;
2812 G4cout<<
"-G4QFragmentation::Breeder:Reduced #"<<astring<<
" & #"<<restr<<
G4endl;
2820 G4cout<<
"G4QFragm::Breeder: StPartner#"<<fustr<<
",LT="<<
LT<<
",RT="<<RT<<
G4endl;
2829 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
2835 if(aPDG > 99) aPDG/=100;
2836 if(rPDG <-99) rPDG=-(-rPDG)/100;
2843 G4cout<<
"G4QFragmentation::Breeder:BeforeFuDir,sL="<<sPDG<<
",nR="<<nPDG<<
",uL="
2844 <<uPDG<<
",mR="<<mPDG<<
",L4M="<<L4M<<
",R4M="<<R4M<<
G4endl;
2853 if ( (uPDG<0 || nPDG<0) && -uPDG==nPDG ) Left->
SetPDGCode(sPDG);
2854 else if( (mPDG<0 || sPDG<0) && -mPDG==sPDG ) Right->
SetPDGCode(nPDG);
2855 else if( (uPDG<0 || sPDG<0) && -uPDG==sPDG ) Left->
SetPDGCode(nPDG);
2856 else if( (mPDG<0 || nPDG<0) && -mPDG==nPDG ) Right->
SetPDGCode(sPDG);
2866 <<
",S="<<L4M+PL4M+R4M+PR4M<<
", L="<<Left->
Get4Momentum()<<
", R="
2870 else if(fusionDONE<0)
2878 <<
",S="<<L4M+PL4M+R4M+PR4M<<
", L="<<Left->
Get4Momentum()<<
", R="
2883 else G4cout<<
"-Warning-G4QFragmentation::Breeder: WrongStringFusion"<<
G4endl;
2886 G4cout<<
"#EMC#G4QFragmentation::Breeder:StringFused,F="<<fusionDONE<<
",L="<<L4M
2887 <<
",R="<<R4M<<
",pL="<<PL4M<<
",pR="<<PR4M<<
",nL="<<Left->
Get4Momentum()
2893 G4cout<<
"###G4QFragmentation::Breeder: Str#"<<astring<<
" fused/w Str#"<<fustr
2911 G4int nHadr=theResult->size();
2918 while( (nHadr=theResult->size()) && !theHadrons)
2921 for(
G4int i=0; i<nHadr; i++)
2924 G4int hPDG=(*theResult)[i]->GetPDGCode();
2926 G4cout<<
"-EMC-G4QFrag::Breed:H#"<<i<<
",4M="<<h4M<<hQC<<
",PDG="<<hPDG<<
G4endl;
2944 for (
G4int reh=0; reh < nHadr; reh++)
2949 if(curPDG==331 && sPDG!=3 && nPDG!=3 && sPDG!=-3 && nPDG!=-3)
2951 if(sPDG==2 || sPDG==-2 || nPDG==2 || nPDG==-2)
2955 else if(curPDG==221 && sPDG!=2 && nPDG!=2 && sPDG!=-2 && nPDG!=-2)
2957 else if(curPDG==111 && sPDG!=1 && nPDG!=1 && sPDG!=-1 && nPDG!=-1)
2962 G4cout<<
"G4QFragmentation::Breeder: Hadron # "<<reh<<
", QC="<<curQC
2963 <<
", P1="<<partPDG1<<
", P2="<<partPDG2<<
G4endl;
2965 if(partPDG1 || partPDG2)
2968 if(sumT>3 && partPDG1 && partPDG2) cCur=2;
2970 G4double M2=(cString4M+curHadr4M).m2();
2972 G4cout<<
"G4QFragmentation::Breeder:*IN*Hadron#"<<reh<<
",M2="<<M2<<
G4endl;
2974 if( (sumT<4 || cCur>=cMax) && M2 > maM2)
2994 G4cout<<
"G4QFrag::Br:*Selected*,P1="<<partPDG1<<
",P2="<<partPDG2<<
G4endl;
3002 G4cout<<
"G4QFragmentation::Breeder: fuh="<<fuhad<<
",fus="<<fusDONE<<
G4endl;
3011 if(secPDG && cMax>1)
3014 G4cout<<
"G4QFragm::Br:TryReduce,nPDG="<<newPDG<<
",sPDG="<<secPDG<<
G4endl;
3020 <<
", L4M="<<newLeft<<
", R4M="<<cRight->
Get4Momentum()<<
", h4M="
3031 <<
", L4M="<<cLeft->
Get4Momentum()<<
", R4M="<<newRight<<
", h4M="
3036 else G4cout<<
"-G4QFragmentation::Breeder: Wrong String+HadronFusion"<<
G4endl;
3039 if(fusDONE)
G4cout<<
"####G4QFragmentation::Breeder: String #"<<astring
3040 <<
" is fused with Hadron #"<<fuhad
3059 G4cout<<
"G4QFragmentation::Breeder: before HR, nH="<<theResult->size()<<
G4endl;
3062 G4QHadronVector::iterator ih;
3066 for(ih = theResult->begin(); ih != theResult->end(); ih++)
3069 G4cout<<
"G4QFrag::Breeder:#"<<icon<<
", i="<<(*ih)<<
", sH="<<selHP<<
G4endl;
3085 G4cout<<
"-EMC->>G4QFragmentation::Breeder: S+=H, 4M="<<curString4M<<
", M="
3086 <<curString4M.
m()<<
", Charge="<<curStrChg<<
", B="<<curStrBaN<<
G4endl;
3089 theResult->erase(ih);
3097 if(!found)
G4cout<<
"*G4QFragmentation::Breeder:nH="<<theResult->size()<<
G4endl;
3102 G4cout<<
"G4QFrag::Breeder: tH="<<theHadrons<<
",nH="<<theResult->size()<<
G4endl;
3106 G4cout<<
"*G4QFragmentation::Breeder: *CanTryToDecay?* nH="<<theHadrons<<
", next="
3107 <<next<<
" =? nS="<<strings.size()<<
", nR="<<theResult->size()<<
G4endl;
3109 if(!theHadrons && next == strings.size() && !(theResult->size()))
3125 if(std::fabs(ttM-h1M-h2M)<=eps)
3128 h14M=part1*curString4M;
3129 h24M=curString4M-h14M;
3135 G4cerr<<
"***G4QFragmentation::Breeder: tM="<<ttM<<
"->h1="<<h1QPDG<<
"("
3136 <<h1M<<
")+h2="<<h1QPDG<<
"("<<h2M<<
")="<<h1M+h2M<<
G4endl;
3141 theResult->push_back(h1H);
3147 G4cout<<
"-EMC->>G4QFragment::Breeder: String=Hadr ChiPro1 is filled, f4M="
3148 <<f4M<<
", fPDG="<<fPD<<
", fCg="<<fCg<<
", fBN="<<fBN<<
G4endl;
3151 theResult->push_back(h2H);
3157 G4cout<<
"-EMC->>G4QFragmentation::Breeder: String=Hadr ChiPro2 is filled, s4M="
3158 <<s4M<<
", sPDG="<<sPD<<
", sCg="<<sCg<<
", sBN="<<sBN<<
G4endl;
3161 G4cout<<
"-EMC-..Chi..G4QFragmentation::Breeder: DecayCHECK, Ch4M="
3162 <<curString4M<<
", d4M="<<curString4M-h14M-h24M<<
G4endl;
3169 theQuasmons.push_back(stringQuasmon);
3175 if (miPDG>0 && miPDG%10 < 3) miPDG+=2;
3176 else if(miPDG<0 && (-miPDG)%10< 3) miPDG-=2;
3180 G4int tmpN=tmpQHadVec->size();
3182 G4cout<<
"G4QFragmentation::Breeder: Decay the Last, Res#H="<<tmpN<<
G4endl;
3186 for(
G4int aH=0; aH < tmpN; aH++)
3188 theResult->push_back((*tmpQHadVec)[aH]);
3195 G4cout<<
"-EMC->>G4QFragment::Breeder:String=Hadr,H#"<<aH<<
" filled, 4M="
3196 <<p4M<<
", PDG="<<PDG<<
", Chg="<<Chg<<
", BaN="<<BaN<<
G4endl;
3204 G4cout<<
"G4QFragmentat::Breeder:==> to Quasm="<<miQC<<curString4M<<
", Nuc="
3205 <<theNucleus<<theNucleus.
Get4Momentum()<<
", NString="<<strings.size()
3206 <<
", nR="<<theResult->size()<<
", nQ="<<theQuasmons.size()<<
G4endl;
3208 theQuasmons.push_back(stringQuasmon);
3210 tmpQHadVec->clear();
3214 tmpQHadVec->clear();
3223 G4cout<<
"G4QFragmentation::Breeder: theH="<<theHadrons<<
"?=0, next="<<next<<
G4endl;
3225 if(!theHadrons && next < strings.size())
3231 G4cout<<
"---->>G4QFragmentation::Breeder: SQC="<<miQC<<
", miSPDG="<<miPDG<<
G4endl;
3242 G4cout<<
"---->>G4QFragmentation::Breeder: minMass="<<miM<<
", realM2="<<cM2<<
G4endl;
3248 if(std::fabs(cM-miM) < eps)
3253 theResult->push_back(sHad);
3255 G4cout<<
"----->>G4QFragmentation::Breeder:S->H="<<miPDG<<curString4M<<
G4endl;
3269 theResult->push_back(h1H);
3275 G4cout<<
"-EMC->>G4QFragmentation::Breeder:Str=2HadrAR Prod-F is filled, f4M="
3276 <<f4M<<
", fPDG="<<fPD<<
", fCg="<<fCg<<
", fBN="<<fBN<<
G4endl;
3279 theResult->push_back(h2H);
3285 G4cout<<
"-EMC->>G4QFragmentation::Breeder:Str=2HadrAR Prod-S is filled, s4M="
3286 <<s4M<<
", sPDG="<<sPD<<
", sCg="<<sCg<<
", sBN="<<sBN<<
G4endl;
3304 G4cout<<
"G4QFr::Breed:TryRecover,V="<<curV<<
",cM2="<<cM2<<
",miM="<<miM<<
G4endl;
3306 nOfStr=strings.size();
3307 for(restr=next; restr < nOfStr; ++restr)
if(restr != astring)
3310 G4cout<<
"G4QFragmentation::Breeder: rS="<<restr<<
", nS="<<nOfStr<<
G4endl;
3314 G4cout<<
"G4QFragmentation::Breeder: pString="<<pString<<
G4endl;
3318 G4cout<<
"G4QFragmentation::Breeder: p4M="<<p4M<<
G4endl;
3325 G4cout<<
"G4QFrag::Breeder: pM2="<<pM2<<
",miM2="<<miM2<<
",cM2="<<cM2<<
G4endl;
3330 G4cout<<
"G4QFragmentation::Breeder: D="<<D<<
",dM4="<<dM4<<
",D2="<<D2<<
G4endl;
3333 if(D > 0. && pM2>.01) x=(std::sqrt(D)-D2)/pM2;
3335 else G4cout<<
"G4QFragment::Breeder: D="<<D<<
",D2="<<D2<<
",pM4="<<dM4<<
G4endl;
3336 G4cout<<
"G4QFragmentation::Breeder: pM2="<<pM2<<
",D2="<<D2<<
",x="<<x<<
G4endl;
3338 if(x > 0. && x < 1.)
3351 if(delta > 0. && delta > maD)
3355 G4cout<<
"G4QFragmentation::Breeder: Subtr,S#"<<restr<<
",d="<<maD<<
G4endl;
3369 G4cout<<
"G4QFragmentation::Breeder: FreeAdd,S#"<<restr<<
",x="<<x<<
G4endl;
3378 G4cout<<
"G4QFragmentation::Breeder:EndOfLOOP r="<<restr<<
"<"<<nOfStr<<
G4endl;
3382 G4cout<<
"G4QFragmentation::Breeder: AfterLOOP fustr="<<fustr<<
G4endl;
3388 G4cout<<
"G4QFragmentation::Breeder: Found Sum4M="<<sum4M<<
G4endl;
3391 curString4M+=selX*s4M;
3392 if(std::abs(miPDG)%10 > 2)
3398 G4cout<<
"G4QFragmentation::Breeder:DecStH,nH="<<tmpQHadVec->size()<<
G4endl;
3400 for(
unsigned aH=0; aH < tmpQHadVec->size(); aH++)
3402 theResult->push_back((*tmpQHadVec)[aH]);
3409 G4cout<<
"-EMC->>G4QFragmentation::Breeder:St=Had,pH#"<<aH<<
" filled, 4M="
3410 <<p4M<<
", PDG="<<PDG<<
", Chg="<<Chg<<
", BaN="<<BaN<<
G4endl;
3413 tmpQHadVec->clear();
3416 else if(miPDG == 10)
3428 if(std::fabs(ttM-h1M-h2M)<=eps)
3431 h14M=part1*curString4M;
3432 h24M=curString4M-h14M;
3438 G4cerr<<
"***G4QFragmentation::Breeder: tM="<<ttM<<
"->h1="<<h1QPDG
3439 <<
"("<<h1M<<
")+h2="<<h1QPDG<<
"("<<h2M<<
")="<<h1M+h2M<<
G4endl;
3444 theResult->push_back(h1H);
3450 G4cout<<
"-EMC->>G4QFragmentation::Breeder:Str=Hadr Prod-F's filled, f4M="
3451 <<f4M<<
", fPDG="<<fPD<<
", fCg="<<fCg<<
", fBN="<<fBN<<
G4endl;
3454 theResult->push_back(h2H);
3460 G4cout<<
"-EMC->>G4QFragmentation::Breeder:Str=Hadr Prod-S's filled, s4M="
3461 <<s4M<<
", sPDG="<<sPD<<
", sCg="<<sCg<<
", sBN="<<sBN<<
G4endl;
3464 G4cout<<
"-EMC-Chipo.G4QFragmentation::Breeder:DecCHECK,c4M="<<curString4M
3465 <<
", ChQC="<<miQC<<
", d4M="<<curString4M-h14M-h24M<<
G4endl;
3470 G4cerr<<
"***G4QFragm::Breeder:tM="<<ttM<<miQC<<
"->h1="<<h1QPDG<<
"(" <<h1M
3471 <<
")+h2="<<h1QPDG<<
"("<<h2M<<
") = "<<h1M+h2M<<
G4endl;
3478 theResult->push_back(sHad);
3480 G4cout<<
"-EMC->>G4QFragmentation::Breeder:Str=Hadr Filled, 4M="
3481 <<curString4M<<
", PDG="<<miPDG<<
G4endl;
3490 G4cout<<
"-EMC-...Cor...G4QFragmentation::Breeder:CorCHECK Sum="<<sum4M
3491 <<
" =? "<<curString4M+pString->
Get4Momentum()<<
", M="<<miM<<
" ?= "
3495 G4cout<<
"---->>G4QFragmentation::Breeder:*Corrected* String->Hadr="<<miPDG
3496 <<curString4M<<
" by String #"<<fustr<<
G4endl;
3504 else G4cout<<
"***G4QFragmentation::Breeder: **No SSCorrection**,next="<<next<<
G4endl;
3511 if (lPDGcS==3 && rPDGcS==-3)
3516 else if(lPDGcS==-3 && rPDGcS==3)
3522 G4int nofRH=theResult->size();
3524 G4cout<<
"G4QFragmentation::Breeder: theH="<<theHadrons<<
", #OfH="<<nofRH<<
G4endl;
3526 if(!theHadrons && nofRH)
3529 G4cout<<
"!G4QFragmentation::Breeder:Can try SHCor, nH="<<theResult->size()<<
G4endl;
3550 for (reha=0; reha < nofRH; reha++)
3557 if(tM2 >=
sqr(pM+miM+eps))
3571 else G4cout<<
"G4QFragmentation::Breeder:H# "<<reha<<
",tM="<<std::sqrt(tM2)<<
" < "
3572 <<
" mS="<<miM<<
" + mH="<<pM<<
" = "<<pM+miM<<
G4endl;
3576 G4cout<<
"G4QFragment::Breeder: fuha="<<fuha<<
", dMmin="<<dMmin<<
G4endl;
3585 if(std::fabs(sM-miM-spM)<=eps)
3595 G4cerr<<
"***G4QFragmentation::Breeder: *SH*, tM="<<sM<<
"->h1=("<<miPDG<<
")"
3596 <<miM<<
" + h2="<<spM<<
" = "<<miM+spM<<
G4endl;
3602 G4cout<<
"-EMC->...>G4QFragmentation::Breeder: H# "<<fuha<<
" is updated, new4M="
3608 G4cerr<<
"***G4QFragm::Breeder: HS Failed, tM="<<sM<<
"->h1M=("<<miPDG<<
")"<<miM
3609 <<
"+h2M="<<spM<<
" = "<<miM+spM<<
G4endl;
3612 if(std::abs(miPDG)%10 > 2)
3618 G4cout<<
"G4QFragment::Breeder:*HS* DecStrHad, nH="<<tmpQHadVec->size()<<
G4endl;
3620 for(
unsigned aH=0; aH < tmpQHadVec->size(); aH++)
3622 theResult->push_back((*tmpQHadVec)[aH]);
3629 G4cout<<
"-EMC->>G4QFragmentation::Breeder: Str+Hadr PrH#"<<aH<<
" filled, 4M="
3630 <<p4M<<
", PDG="<<PDG<<
", Chg="<<Chg<<
", BaN="<<BaN<<
G4endl;
3633 tmpQHadVec->clear();
3636 else if(miPDG == 10)
3648 if(std::fabs(ttM-h1M-h2M)<=eps)
3658 G4cerr<<
"***G4QFragmentation::Breeder: HS tM="<<ttM<<
"->h1="<<h1QPDG<<
"("
3659 <<h1M<<
")+h2="<<h1QPDG<<
"("<<h2M<<
")="<<h1M+h2M<<
G4endl;
3664 theResult->push_back(h1H);
3670 G4cout<<
"-EMC->>G4QFragmentation::Breeder: CorStrHadr Prod-1 is filled, f4M="
3671 <<f4M<<
", fPDG="<<fPD<<
", fCg="<<fCg<<
", fBN="<<fBN<<
G4endl;
3674 theResult->push_back(h2H);
3680 G4cout<<
"-EMC->>G4QFragmentation::Breeder: CorStrHadr Prod-2 is filled, n4M="
3681 <<n4M<<
", nPDG="<<nPD<<
", nCg="<<nCg<<
", nBN="<<nBN<<
G4endl;
3684 G4cout<<
"-EMC-...HS-Chipo...G4QFragmentation::Breeder:DecCHECK, Ch4M="<<mi4M
3685 <<
", ChQC="<<miQC<<
", d4M="<<mi4M-h14M-h24M<<
G4endl;
3692 theResult->push_back(sHad);
3694 G4cout<<
"----->>G4QFragmentation::Breeder: CorStr=Hadr is Filled, 4M="
3695 <<curString4M<<
", StPDG="<<miPDG<<
G4endl;
3699 G4cout<<
"-EMC-...Cor...G4QFragmentation::Breeder:StHadCor CHECK Sum="<<s4M
3703 G4cout<<
"------->>G4QFragmentation::Breeder: *Corrected* String+Hadr="<<miPDG
3704 <<mi4M<<
" by Hadron #"<<reha<<
G4endl;
3711 G4cout<<
"G4QFragmentation::Breeder: Str+Hadr Failed, 4M="<<curString4M
3712 <<
", PDG="<<miPDG<<
" -> Now try to recover the string as a hadron"<<
G4endl;
3739 theQuasmons.push_back(stringQuasmon);
3745 if(theHadrons) nHfin=theHadrons->size();
3750 G4int tnSt=strings.size();
3751 for(
G4int i=astring; i < tnSt; ++i)
3758 G4cout<<
"=--=>G4QFragmentation::Breeder:S#"<<i<<
",4M="<<pS4M<<
",QC="<<pSQC<<
G4endl;
3762 G4cout<<
"==>G4QFragmentation::Breeder:AllStrings are summed up in a Quasmon"<<
G4endl;
3765 theQuasmons.push_back(stringQuasmon);
3769 G4cout<<
"G4QFragmentation::Breeder: Trying to decay hadrons #ofHRes="<<nHfin<<
G4endl;
3771 for(
G4int aTrack=0; aTrack<nHfin; aTrack++)
3773 G4QHadron* curHadron=(*theHadrons)[aTrack];
3779 G4cout<<
"----->>G4QFragmentation::Breeder:S#"<<astring<<
",H#"<<aTrack<<
",PDG="<<hPDG
3782 if(std::abs(hPDG)%10 > 2)
3786 G4cout<<
"G4QFragmentation::Breeder:-DECAY'S DONE-,nH="<<tmpQHadVec->size()<<
G4endl;
3788 for(
unsigned aH=0; aH < tmpQHadVec->size(); aH++)
3790 theResult->push_back((*tmpQHadVec)[aH]);
3804 G4cout<<
"-EMC->>G4QFragmentation::Breeder:String*Filled, 4M="<<p4M<<
", PDG="<<PDG
3805 <<
", Chg="<<Chg<<
", BaN="<<BaN<<
G4endl;
3809 G4cout<<
"-EMC-.G4QFr::Br:Dec,r4M="<<curH4M<<
",rC="<<curHCh<<
",rB="<<curHBN<<
G4endl;
3811 tmpQHadVec->clear();
3816 theResult->push_back(curHadron);
3818 curString4M-=curH4M;
3824 G4cout<<
"-EMC->>-->>G4QFragmentation::Breeder: curH filled 4M="<<curH4M<<
",PDG="
3830 if(theHadrons)
delete theHadrons;
3832 G4cout<<
"-EMC-.........G4QFragmentation::Breeder: StringDecay CHECK, r4M="<<curString4M
3833 <<
", rChg="<<curStrChg<<
", rBaN="<<curStrBaN<<
G4endl;
3836 if(curStrChg || curStrBaN || curString4M.
t() > eps || std::fabs(curString4M.
x()) > eps
3837 || std::fabs(curString4M.
y()) > eps || std::fabs(curString4M.
z()) > eps )
3843 G4int nHadr=theResult->size();
3856 for(
G4int i=0; i<nHadr; i++)
3872 G4cout<<
"G4QFragmentation::Breeder: H#"<<i<<
", d4M="<<curString4M+hI4M
3873 <<
",dCh="<<hCh+curStrChg<<
",dBN="<<hBN+curStrBaN<<
G4endl;
3874 if( !(hCh+curStrChg) && !(hBN+curStrBaN) && std::fabs(dEn+hEn)<eps &&
3875 std::fabs(dPx+hPx)<eps && std::fabs(dPy+hPy)<eps && std::fabs(dPz+hPz)<eps )
3877 G4cout<<
"G4QFragmentation::Breeder: ***Cured*** Redundent Hadron # "<<i<<
G4endl;
3878 G4QHadron* theLast = (*theResult)[nHadr-1];
3883 theResult->pop_back();
3887 if( !(hCh+mCh+curStrChg) && !(hBN+mBN+curStrBaN) && std::fabs(dEn+hEn+mEn)<eps &&
3888 std::fabs(dPx+hPx+mPx)<eps && std::fabs(dPy+hPy+mPy)<eps &&
3889 std::fabs(dPz+hPz+mPz)<eps && i>0)
3891 G4cout<<
"G4QFragmentation::Breeder:***Cured*** Redundent 2Hadrons i="<<i<<
G4endl;
3893 G4QHadron* theLast = (*theResult)[nHadr-1];
3901 theResult->pop_back();
3903 theLast = (*theResult)[nHadr-2];
3911 theResult->pop_back();
3913 nHadr=theResult->size();
3917 G4cout<<
"*Warning*G4QFragmentation::Breeder: Nonconservation isn't cured!"<<
G4endl;
3925 theResult->push_back(resNuc);
3930 G4int nHadr=theResult->size();
3931 G4int nQuasm=theQuasmons.size();
3932 G4cout<<
"-EMCLS-G4QFragmentation::Breeder:#ofHadr="<<nHadr<<
", #OfQuasm="<<nQuasm<<
",rN="
3934 for(
G4int i=0; i<nHadr; i++)
3938 G4int hChg=(*theResult)[i]->GetCharge();
3940 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
3942 G4cout<<
"-EMCLS-G4QFragmentation::Breeder:(1) Hadron#"<<i<<
", 4M="<<hI4M<<
", PDG="
3943 <<(*theResult)[i]->GetPDGCode()<<
", C="<<hChg<<
", B="<<hBaN<<
G4endl;
3945 for(
G4int i=0; i<nQuasm; i++)
3949 G4int hChg=theQuasmons[i]->GetCharge();
3951 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
3953 G4cout<<
"-EMCLS-G4QFragmentation::Breeder:(1) Quasmon#"<<i<<
", 4M="<<hI4M<<
", C="<<hChg
3956 G4cout<<
"-EMCLS-G4QFragm::Breed: LS r4M="<<s4M-totLS4M<<
",rC="<<rCh<<
",rB="<<rBN<<
G4endl;
3959 G4int nRes=theResult->size();
3961 G4cout<<
"G4QFragmentation::Breeder: Strings4M="<<ft4M<<
", nRes="<<nRes<<
G4endl;
3965 if(nRes > 2 && maxEn > 0.)
3967 std::list<std::pair<G4double, G4QHadron*> > theSorted;
3968 std::list<std::pair<G4double, G4QHadron*> >::iterator current;
3969 for(
G4int secondary = 0; secondary<nRes-1; ++secondary)
3971 G4QHadron* ih =theResult->operator[](secondary);
3976 if(hM2>0.00001) toSort=h4M.
e()+h3M.
dot(LSDir)/std::sqrt(hM2);
3978 G4cout<<
"G4QFragmentation::Breeder:#"<<secondary<<
",M2="<<hM2<<
",s="<<toSort<<
G4endl;
3980 std::pair<G4double, G4QHadron*> it;
3984 for(current = theSorted.begin(); current!=theSorted.end(); ++current)
3986 if((*current).first > toSort)
3988 theSorted.insert(current, it);
3993 if(!inserted) theSorted.push_back(it);
3998 for(current = theSorted.begin(); current!=theSorted.end(); ++current)
4005 if (hPDG> 1111 && hPDG< 3335) dE=h4M.
e()-ih->
GetMass();
4006 else if(hPDG>-1111 && hPDG<1111 && hPDG!=22) dE=h4M.
e();
4010 G4cout<<
"G4QFragmentation::Breeder:dE="<<dE<<
",mE="<<maxEn<<
",t="<<tested<<
G4endl;
4012 if(tested && dE < maxEn)
4018 G4cout<<
"%->G4QFragmentation::Breeder:Exclude,4M="<<h4M<<
",dE="<<maxEn<<
G4endl;
4022 else theResult->push_back(ih);
4026 G4cout<<
"%->G4QFragmentation::Breeder:QuasmonIsFilled,4M="<<q4M<<
",QC="<<qQC<<
G4endl;
4028 if(q4M != vac4M) theQuasmons.push_back(softQuasmon);
4029 else delete softQuasmon;
4030 theResult->push_back(resNuc);
4035 G4int nHd=theResult->size();
4036 G4int nQm=theQuasmons.size();
4037 G4cout<<
"-EMCLS-G4QFragmentation::Breeder:#ofHadr="<<nHd<<
", #OfQuasm="<<nQm<<
",rN="
4039 for(
G4int i=0; i<nHd; i++)
4043 G4int hChg=(*theResult)[i]->GetCharge();
4045 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
4047 G4cout<<
"-EMCLS-G4QFragmentation::Breeder:(2) Hadron#"<<i<<
", 4M="<<hI4M<<
", PDG="
4048 <<(*theResult)[i]->GetPDGCode()<<
", C="<<hChg<<
", B="<<hBaN<<
G4endl;
4050 for(
G4int i=0; i<nQm; i++)
4054 G4int hChg=theQuasmons[i]->GetCharge();
4056 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
4058 G4cout<<
"-EMCLS-G4QFragmentation::Breeder:(2) Quasmon#"<<i<<
", 4M="<<hI4M<<
", C="
4059 <<hChg<<
", B="<<hBaN<<
G4endl;
4061 G4cout<<
"-EMCLS-G4QFragm::Breed:, r4M="<<f4M-totLS4M<<
",rC="<<fCh<<
",rB="<<fBN<<
G4endl;
4073 G4double Mprojectile2=Mprojectile*Mprojectile;
4078 G4cout<<
"G4QFragm::ExciteDiffPartici:Ep="<<Pprojectile.
e()<<
",Et="<<Ptarget.
e()<<
G4endl;
4087 G4cout<<
"G4QFragmentation::ExciteDiffParticipants: *1* abort Collision!! *1*"<<
G4endl;
4094 G4cout<<
"G4QFragment::ExciteDiffParticipantts:BeforBoost Pproj="<<Pprojectile<<
", Ptarg="
4101 G4cout<<
"G4QFragment::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<
"Ptarg="
4102 <<Ptarget<<
", cms4M="<<Pprojectile+Ptarget<<
G4endl;
4103 G4cout<<
"G4QFragment::ExciteDiffParticipants: ProjX+="<<Pprojectile.
plus()<<
", ProjX-="
4104 <<Pprojectile.
minus()<<
", TargX+="<< Ptarget.
plus()<<
", TargX-="<<Ptarget.
minus()
4110 G4cout<<
"G4QFragmentation::ExciteDiffParticipants: Before DO"<<
G4endl;
4117 G4cout<<
"G4QFragmentation::ExciteDiffParticipants: maxPtSq="<<maxPtSquare<<
G4endl;
4118 if(whilecount++>=500 && whilecount%100==0)
4119 G4cout<<
"G4QFragmentation::ExciteDiffParticipantts: can loop, loopCount="<<whilecount
4120 <<
", maxPtSquare="<<maxPtSquare<<
G4endl;
4125 G4cout<<
"G4QFragmentation::ExciteDiffParticipants: *2* abort Loop!! *2*"<<
G4endl;
4131 G4cout<<
"G4QFragment::ExciteDiffParticipants: generated Pt="<<Qmomentum<<
", ProjPt="
4132 <<Pprojectile+Qmomentum<<
", TargPt="<<Ptarget-Qmomentum<<
G4endl;
4140 G4cout<<
"G4QFragment::ExciteDiffParticip: X-plus="<<Xplus<<
",X-minus="<<Xminus<<
G4endl;
4145 Qmomentum.
setPz((Qplus-Qminus)/2);
4146 Qmomentum.
setE( (Qplus+Qminus)/2);
4148 G4cout<<
"G4QFragment::ExciteDiffParticip: Qplus="<<Qplus<<
", Qminus="<<Qminus<<
", pt2="
4149 <<pt2<<
", Qmomentum="<<Qmomentum<<
", ProjM="<<(Pprojectile+Qmomentum).mag()
4150 <<
", TargM="<<(Ptarget-Qmomentum).mag()<<
G4endl;
4152 }
while((Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
4153 (Ptarget-Qmomentum).mag2()<=Mtarget2);
4154 Pprojectile += Qmomentum;
4155 Ptarget -= Qmomentum;
4157 G4cout<<
"G4QFragment::ExciteDiffParticipan: Proj(Q)="<<Pprojectile<<
", Targ(Q)="<<Ptarget
4158 <<
", Proj(back)="<<toLab*Pprojectile<<
", Targ(bac)="<< toLab*Ptarget <<
G4endl;
4164 G4cout<<
"G4QFragmentation::ExciteDiffParticipants: TargetMass="<<Ptarget.
mag()<<
G4endl;
4168 G4cout<<
"G4QFragment::ExciteDiffParticipants:ProjectileMass="<<Pprojectile.
mag()<<
G4endl;
4181 G4double Mprojectile2=Mprojectile*Mprojectile;
4186 G4cout<<
"G4QFragm::ExSingDiffPartici:Ep="<<Pprojectile.
e()<<
",Et="<<Ptarget.
e()<<
G4endl;
4193 G4cout<<
"--1/2--G4QFragmentation::ExSingDiffParticipants: Projectile is fixed"<<
G4endl;
4195 Mprojectile2 = projectile->
GetMass2()*(1.-perCent);
4200 G4cout<<
"---1/2---G4QFragmentation::ExSingDiffParticipants: Target is fixed"<<
G4endl;
4202 Mtarget2 = target->
GetMass2()*(1.-perCent);
4212 G4cout<<
"G4QFragment::ExciteSingDiffParticipants: *1* abort Collision!! *1*"<<
G4endl;
4219 G4cout<<
"G4QFragm::ExciteSingDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<
",Ptarg="
4226 G4cout<<
"G4QFragment::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<
"Ptarg="
4227 <<Ptarget<<
", cms4M="<<Pprojectile+Ptarget<<
G4endl;
4229 G4cout<<
"G4QFragment::ExciteDiffParticipantts: ProjX+="<<Pprojectile.
plus()<<
", ProjX-="
4230 <<Pprojectile.
minus()<<
", TargX+="<< Ptarget.
plus()<<
", TargX-="<<Ptarget.
minus()
4239 if(whilecount++>=500 && whilecount%100==0)
4241 G4cout<<
"G4QFragment::ExciteSingDiffParticipantts: can loop, loopCount="<<whilecount
4242 <<
", maxPtSquare="<<maxPtSquare<<
G4endl;
4247 G4cout<<
"G4QFragmentation::ExciteSingDiffParticipants: *2* abort Loop!! *2*"<<
G4endl;
4253 G4cout<<
"G4QFragm::ExciteSingDiffParticipants: generated Pt="<<Qmomentum<<
", ProjPt="
4254 <<Pprojectile+Qmomentum<<
", TargPt="<<Ptarget-Qmomentum<<
G4endl;
4262 G4cout<<
"G4QFragm::ExciteSingDiffPartici:X-plus="<<Xplus<<
",X-minus="<<Xminus<<
G4endl;
4268 Qminus=(projectile->
GetMass2()+pt2)/(Pprojectile.
plus()+Qplus) - Pprojectile.
minus();
4270 Qmomentum.
setPz((Qplus-Qminus)/2);
4271 Qmomentum.
setE( (Qplus+Qminus)/2);
4273 G4cout<<
"G4QFragm::ExciteDiffParticip: Qplus="<<Qplus<<
", Qminus="<<Qminus<<
", pt2="
4274 <<pt2<<
", Qmomentum="<<Qmomentum<<
", ProjM="<<(Pprojectile+Qmomentum).mag()
4275 <<
", TargM="<<(Ptarget-Qmomentum).mag()<<
G4endl;
4280 }
while((Ptarget-Qmomentum).mag2()<=Mtarget2 ||
4281 (Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
4282 (Ptarget-Qmomentum).e() < 0. || (Pprojectile+Qmomentum).e() < 0.);
4283 Pprojectile += Qmomentum;
4284 Ptarget -= Qmomentum;
4286 G4cout<<
"G4QFragmentation::ExciteSingDiffParticipan: Proj(Q)="<<Pprojectile<<
"(E="
4287 <<Pprojectile.
e()<<
"), Targ(Q)="<<Ptarget<<
"(E="<<Ptarget.
e()
4288 <<
"), Proj(back)="<<toLab*Pprojectile<<
", Targ(bac)="<< toLab*Ptarget <<
G4endl;
4294 G4cout<<
"G4QFragm::ExciteSingleDiffParticipants: TargetMass="<<Ptarget.
mag()<<
G4endl;
4298 G4cout<<
"G4QFragm::ExciteSingleParticipants:ProjectileMass="<<Pprojectile.
mag()<<
G4endl;
4307 stringTension = stT;
4309 widthOfPtSquare = -2*SigPt*SigPt;
4316 if(Xmax == Xmin)
return Xmin;
4317 if( Xmin < 0. || Xmax < Xmin)
4319 G4cerr<<
"***G4QFragmentation::ChooseX: Xmin="<<Xmin<<
", Xmax="<<Xmax<<
G4endl;
4327 G4cout<<
"G4QFragmentation::ChooseX: DiffractiveX="<<x<<
G4endl;
4336 G4cout<<
"G4QFragmentation::GaussianPt:width2="<<widthSq<<
",maxPt2="<<maxPtSquare<<
G4endl;
4340 if(rm>-.01) pt2=widthSq*(std::sqrt(1.-
G4UniformRand()*(1.-
sqr(1.+rm)))-1.);
4341 else pt2=widthSq*std::log(1.-
G4UniformRand()*(1.-std::exp(rm)));
4344 return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.);
4349 if (PDG1 < 7 && PDG1 > 0 && PDG2 < 7 && PDG2 > 0)
4351 if(PDG1 > PDG2)
return PDG1*1000+PDG2*100+1;
4352 else return PDG2*1000+PDG1*100+1;
4354 else if (PDG1 >-7 && PDG1 < 0 && PDG2 >-7 && PDG2 < 0)
4356 if(-PDG1 > -PDG2)
return PDG1*1000+PDG2*100-1;
4357 else return PDG2*1000+PDG1*100-1;
4359 else if (PDG1 <-99 && PDG2 < 7 && PDG2 > 0)
4361 G4int PDG=-PDG1/100;
4362 if(PDG2==PDG/10)
return -PDG%10;
4363 if(PDG2==PDG%10)
return -PDG/10;
4366 G4cerr<<
"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<
", PDG2="<<PDG2<<
G4endl;
4370 else if (PDG2 <-99 && PDG1 < 7 && PDG1 > 0)
4372 G4int PDG=-PDG2/100;
4373 if(PDG1==PDG/10)
return -PDG%10;
4374 if(PDG1==PDG%10)
return -PDG/10;
4377 G4cerr<<
"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<
", PDG2="<<PDG2<<
G4endl;
4381 else if (PDG1 > 99 && PDG2 >-7 && PDG2 < 0)
4384 if(PDG2==-PDG/10)
return PDG%10;
4385 if(PDG2==-PDG%10)
return PDG/10;
4388 G4cerr<<
"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<
", PDG2="<<PDG2<<
G4endl;
4392 else if (PDG2 > 99 && PDG1 >-7 && PDG1 < 0)
4395 if(PDG1==-PDG/10)
return PDG%10;
4396 if(PDG1==-PDG%10)
return PDG/10;
4399 G4cerr<<
"***G4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<
", PDG2="<<PDG2<<
G4endl;
4405 G4cerr<<
"***G4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<
", PDG2="<<PDG2<<
G4endl;
4415 G4cout<<
"G4QFragmentation::ReducePair: **Called** P1="<<P1<<
", P2="<<P2<<
G4endl;
4421 if (P11==P21)
return std::make_pair(P12,P22);
4422 else if(P11==P22)
return std::make_pair(P12,P21);
4423 else if(P12==P21)
return std::make_pair(P11,P22);
4424 else if(P12==P22)
return std::make_pair(P11,P21);
4426 G4cout<<
"-Warning-G4QFragmentation::ReducePair:**Failed**, P1="<<P1<<
", P2="<<P2<<
G4endl;
4428 return std::make_pair(0,0);
4437 if (LS==2 && MPS==2 )
4440 G4cout<<
"G4QFragmentation::AnnihilationOrder: QaQ/QaQ->DiQ-aDiQ, uPDG="<<uPDG<<
G4endl;
4442 if ( (uPDG>0 && sPDG>0 && mPDG<0 && nPDG<0) ||
4443 (uPDG<0 && sPDG<0 && mPDG>0 && nPDG>0) ) Ord= 1;
4444 else if( (uPDG>0 && nPDG>0 && mPDG<0 && sPDG<0) ||
4445 (uPDG<0 && nPDG<0 && mPDG>0 && sPDG>0) ) Ord=-1;
4446 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 22 fusion, L="<<uPDG
4447 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4449 else if ( LS==2 && MPS==3 )
4454 G4cout<<
"G4QFragmentation::AnnihOrder: pLDiQ, sPDG="<<sPDG<<
", nPDG="<<nPDG<<
G4endl;
4456 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) ) Ord= 1;
4457 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) ) Ord=-1;
4458 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong pLDiQ, L="<<uPDG
4459 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4464 G4cout<<
"G4QFragmentation::AnnihOrder: pRDiQ, sPDG="<<sPDG<<
", nPDG="<<nPDG<<
G4endl;
4466 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) ) Ord=-1;
4467 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) ) Ord= 1;
4468 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong pRDiQ, L="<<uPDG
4469 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4474 G4cout<<
"G4QFragmentation::AnnihOrder: pLaDiQ, sPDG="<<sPDG<<
", nPDG="<<nPDG<<
G4endl;
4476 if ( sPDG>0 && (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1;
4477 else if( nPDG>0 && (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
4478 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong pLaDiQ, L="<<uPDG
4479 <<
", R="<<mPDG<<
", cL="<<sPDG<<
", cR="<<nPDG<<
G4endl;
4484 G4cout<<
"G4QFragmentation::AnnihOrder: pRaDiQ, sPDG="<<sPDG<<
", nPDG="<<nPDG<<
G4endl;
4486 if ( sPDG>0 && (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
4487 else if( nPDG>0 && (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
4488 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong pRaDiQ, L="<<uPDG
4489 <<
", R="<<mPDG<<
", cL="<<sPDG<<
", cR="<<nPDG<<
G4endl;
4491 else if( (sPDG<0 && (-sPDG==mPDG || -sPDG==uPDG) ) ||
4492 (nPDG<0 && (-nPDG==mPDG || -nPDG==uPDG) ) ) Ord= 2;
4494 else G4cout<<
"-Warning-G4QFragmentation::AnnihilatOrder: Wrong 23StringFusion"<<
G4endl;
4495 G4cout<<
"G4QFragmentation::AnnihilationOrder: Ord="<<Ord<<
",sPDG="<<sPDG<<
",nPDG="
4496 <<nPDG<<
", uPDG="<<uPDG<<
",mPDG="<<mPDG<<
G4endl;
4499 else if ( LS==3 && MPS==2 )
4504 G4cout<<
"G4QFragmentation::AnnihOrder: cLDiQ, uPDG="<<uPDG<<
", mPDG="<<mPDG<<
G4endl;
4506 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) ) Ord= 1;
4507 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) ) Ord=-1;
4508 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong cLDiQ, L="<<uPDG
4509 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4514 G4cout<<
"G4QFragmentation::AnnihOrder: cRDiQ, uPDG="<<uPDG<<
", mPDG="<<mPDG<<
G4endl;
4516 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) ) Ord=-1;
4517 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) ) Ord= 1;
4518 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong cRDiQ, L="<<uPDG
4519 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4524 G4cout<<
"G4QFragmentation::AnnihOrder: cLaDiQ, uPDG="<<uPDG<<
", mPDG="<<mPDG<<
G4endl;
4526 if ( uPDG>0 && (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
4527 else if( mPDG>0 && (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
4528 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong cLaDiQ, L="<<uPDG
4529 <<
", R="<<mPDG<<
", cL="<<sPDG<<
", cR="<<nPDG<<
G4endl;
4534 G4cout<<
"G4QFragmentation::AnnihOrder: cRaDiQ, uPDG="<<uPDG<<
", mPDG="<<mPDG<<
G4endl;
4536 if ( uPDG>0 && (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
4537 else if( mPDG>0 && (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
4538 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong cRaDiQ, L="<<uPDG
4539 <<
", R="<<mPDG<<
", cL="<<sPDG<<
", cR="<<nPDG<<
G4endl;
4541 else if( (uPDG<0 && (-uPDG==sPDG || -uPDG==nPDG) ) ||
4542 (mPDG<0 && (-mPDG==sPDG || -mPDG==nPDG) ) ) Ord=2;
4544 else G4cout<<
"-Warning-G4QFragmentation::AnnihilatOrder: Wrong 32StringFusion"<<
G4endl;
4545 G4cout<<
"G4QFragmentation::AnnihilationOrder: Ord="<<Ord<<
",sPDG="<<sPDG<<
",nPDG="
4546 <<nPDG<<
", uPDG="<<uPDG<<
",mPDG="<<mPDG<<
G4endl;
4549 else if ( (LS==2 && MPS==4) || (LS==4 && MPS==2) )
4554 G4cout<<
"G4QFragmentation::AnnihilOrder:pLDiQ, sPDG="<<sPDG<<
",nPDG="<<nPDG<<
G4endl;
4556 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
4557 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
4558 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
4559 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
4560 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 24 pLDiQ, L="<<uPDG
4561 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4566 G4cout<<
"G4QFragmentation::AnnihilOrder:PRDiQ, sPDG="<<sPDG<<
",nPDG="<<nPDG<<
G4endl;
4568 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
4569 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
4570 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
4571 (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1;
4572 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 24 pLDiQ, L="<<uPDG
4573 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4578 G4cout<<
"G4QFragmentation::AnnihilOrder:cLDiQ, uPDG="<<uPDG<<
",mPDG="<<mPDG<<
G4endl;
4580 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
4581 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
4582 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
4583 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
4584 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 24 cLDiQ, L="<<uPDG
4585 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4590 G4cout<<
"G4QFragmentation::AnnihilOrder:cRDiQ, uPDG="<<uPDG<<
",mPDG="<<mPDG<<
G4endl;
4592 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) &&
4593 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
4594 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
4595 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
4596 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 24 cRDiQ, L="<<uPDG
4597 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4600 else G4cout<<
"-Warning-G4QFragmentation::AnnihilOrder: Wrong 24 StringFusion"<<
G4endl;
4601 G4cout<<
"G4QFragmentation::AnnihilationOrder: Ord="<<Ord<<
",sPDG="<<sPDG<<
",nPDG="
4602 <<nPDG<<
", uPDG="<<uPDG<<
",mPDG="<<mPDG<<
G4endl;
4605 else if ( LS==3 && MPS==3 )
4610 G4cout<<
"G4QFragmentation::AnnihilOrder: pLDiQ, sPDG="<<sPDG<<
",nPDG="<<nPDG<<
G4endl;
4612 if ( sPDG<-7 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
4613 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
4614 else if( nPDG<-7 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
4615 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
4616 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 33 pLDiQ, L="<<uPDG
4617 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4622 G4cout<<
"G4QFragmentation::AnnihilOrder: pRDiQ, sPDG="<<sPDG<<
",nPDG="<<nPDG<<
G4endl;
4624 if ( sPDG<-7 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
4625 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
4626 else if( nPDG<-7 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
4627 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
4628 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 33 pRDiQ, L="<<uPDG
4629 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4634 G4cout<<
"G4QFragmentation::AnnihilOrder: cLDiQ, uPDG="<<uPDG<<
",mPDG="<<mPDG<<
G4endl;
4636 if ( uPDG<-7 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
4637 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
4638 else if( mPDG<-7 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
4639 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
4640 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 33 cLDiQ, L="<<uPDG
4641 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4646 G4cout<<
"G4QFragmentation::AnnihilOrder: cRDiQ, uPDG="<<uPDG<<
",mPDG="<<mPDG<<
G4endl;
4648 if ( uPDG<-7 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
4649 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord= 1;
4650 else if( mPDG<-7 && (-uPDG==nPDG/1000 || -sPDG==(nPDG/100)%10) &&
4651 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
4652 else G4cerr<<
"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 33 cRDiQ, L="<<uPDG
4653 <<
",R="<<mPDG<<
",cL="<<sPDG<<
",cR="<<nPDG<<
G4endl;
4656 else G4cout<<
"-Warning-G4QFragmentation::AnnihilOrder: Wrong 33 StringFusion"<<
G4endl;
4657 G4cout<<
"G4QFragmentation::AnnihilationOrder: Ord="<<Ord<<
",sPDG="<<sPDG<<
",nPDG="
4658 <<nPDG<<
", uPDG="<<uPDG<<
",mPDG="<<mPDG<<
G4endl;
4667 G4QStringVector::iterator ist;
4668 for(ist = strings.begin(); ist < strings.end(); ist++)
4670 G4QParton* cLeft=(*ist)->GetLeftParton();
4671 G4QParton* cRight=(*ist)->GetRightParton();
4676 if(std::fabs(cSM2)<.01)
4693 G4QStringVector::iterator sst;
4694 G4QStringVector::iterator pst;
4698 G4cout<<
"G4QFragmentation::SwapPartons: M2=="<<cSM2<<
", 4M="<<cS4M<<
",LPDG="<<cLPDG
4699 <<
",RPDG="<<cRPDG<<
G4endl;
4701 for(pst = strings.begin(); pst < strings.end(); pst++)
if(pst != ist)
4703 G4QParton* pLeft=(*pst)->GetLeftParton();
4704 G4QParton* pRight=(*pst)->GetRightParton();
4712 G4cout<<
"G4QFragmentation::SwapPartons: p4M="<<cS4M<<
",pM2="<<cS4M.
m2()<<
",LPDG="
4713 <<pLPDG<<
",RPDG="<<pRPDG<<
G4endl;
4717 if(((cLPDG<-7 || (cLPDG>0 && cLPDG< 7) ) && (pLPDG<-7 || (pLPDG>0 && pLPDG< 7) ))||
4718 ((cLPDG> 7 || (cLPDG<0 && cLPDG>-7) ) && (pLPDG> 7 || (pLPDG<0 && cLPDG>-7) )))
4722 if(pLM2>0. && cLM2>0.)
4725 if(cLT+pRT==3) pLM-=baryM;
4727 if(cRT+pLT==3) cLM-=baryM;
4728 LM=std::min(pLM2,cLM2);
4731 if(((cRPDG<-7 || (cRPDG>0 && cRPDG< 7) ) && (pRPDG<-7 || (pRPDG>0 && pRPDG< 7) ))||
4732 ((cRPDG> 7 || (cRPDG<0 && cRPDG>-7) ) && (pRPDG> 7 || (pRPDG<0 && cRPDG>-7) )) )
4736 if(pRM2>0. && cRM2>0.)
4739 if(cRT+pLT==3) pRM-=baryM;
4741 if(cLT+pRT==3) cRM-=baryM;
4742 RM=std::min(pRM,cRM);
4766 G4QParton* pLeft=(*sst)->GetLeftParton();
4767 G4QParton* pRight=(*sst)->GetRightParton();
4771 (*sst)->SetLeftParton(cLeft);
4772 (*ist)->SetLeftParton(swap);
4777 (*sst)->SetRightParton(cRight);
4778 (*ist)->SetRightParton(swap);
4782 else G4cout<<
"***G4QFragmentation::SwapPartons:**Failed**,cLPDG="<<cLPDG<<
",cRPDG="
4783 <<cRPDG<<
",-->cM2="<<cSM2<<
G4endl;
4796 static const G4double mAlPr = mAlph+mProt;
4797 static const G4double mAlNt = mAlph+mNeut;
4798 static const G4double dProt = mProt+mProt;
4799 static const G4double dNeut = mNeut+mNeut;
4800 static const G4double dAlph = mAlph+mAlph;
4809 <<
",QC="<<theQC<<
", BN="<<theBN<<
G4endl;
4824 if(
sqr(h1M+h2M) < chM2 )
4828 if(!
G4QHadron(ch4M).DecayIn2(h14M,h24M))
4830 G4cerr<<
"***G4QFrag::EvaporateResid: CM="<<std::sqrt(chM2)<<
" -> h1="<<h1QPDG<<
"("
4831 <<h1M<<
") + h2="<<h1QPDG<<
"("<<h2M<<
") = "<<h1M+h2M<<
" **Failed**"<<
G4endl;
4833 G4Exception(
"G4QFragmentation::EvaporateResidual()",
"HAD_CHPS_0000",
4838 theResult->push_back(h1H);
4840 G4cout<<
"G4QFragm::EvaporateResidual: Chipolino -> H1="<<h1QPDG<<h14M<<
G4endl;
4843 theResult->push_back(qH);
4845 G4cout<<
"G4QE::EvaporateResidual: Chipolino -> H2="<<h2QPDG<<h24M<<
G4endl;
4851 <<
", chipoM="<<std::sqrt(chM2)<<
" < m1="<<h1M<<
"("<<h1QPDG<<
") + m2="<<h2M
4852 <<
"("<<h2QPDG<<
") = "<<h1M+h2M<<
G4endl;
4854 G4Exception(
"G4QFragmentation::EvaporateResidual()",
"HAD_CHPS_0001",
4879 G4cout<<
"G4QFragmentation::EvaRes:(!)Meson(!) PDG="<<thePDG<<
",4M="<<mesLV<<mesLV.
m()
4887 G4cout<<
"G4QFragment::EvaRes: qH.Charge = "<<theC<<
G4endl;
4890 if( thePDG == 10 && theBN > 0 ) thePDG=theQC.
GetZNSPDGCode();
4891 if(theS>0) thePDG-=theS*999999;
4893 G4cout<<
"G4QFragment::EvaRes: S="<<theS<<
", newPDG="<<thePDG<<
G4endl;
4897 G4cout<<
"G4QFragment::EvaRes: totGSM="<<totGSM<<
G4endl;
4901 if(!theC) totGSM=dNeut;
4902 else if(theC==2) totGSM=dProt;
4907 if (theC==3) totGSM=mAlPr;
4908 else if(theC==2) totGSM=mAlNt;
4910 else if(theBN==8) totGSM=dAlph;
4915 G4cout<<
"G4QFragment::EvaRes: Excitation = "<<totMass-totGSM<<
G4endl;
4917 if(std::fabs(totMass-totGSM) < eps)
4919 theResult->push_back(qH);
4921 else if(totMass > totGSM)
4924 G4cout<<
"G4QFragment::EvaRes: try Evaporate Nucleus PDG="<<thePDG<<
G4endl;
4928 G4cout<<
"G4QFragment::EvaRes: ** Evaporation is done **"<<
G4endl;
4936 G4cout<<
"-War-G4QFr::EvaRes:*Must correct* "<<theQC<<q4M<<totMass<<
"<"<<totGSM<<
G4endl;
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4QHadron * > G4QHadronVector
std::vector< G4QInteraction * > G4QInteractionVector
std::vector< G4QPartonPair * > G4QPartonPairVector
CLHEP::Hep3Vector G4ThreeVector
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
double dot(const Hep3Vector &) const
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
Hep3Vector boostVector() const
HepLorentzVector & boost(double, double, double)
void setV(const Hep3Vector &)
HepLorentzVector & transform(const HepRotation &)
G4QParticle * GetQParticle(G4int PDG) const
static G4QCHIPSWorld * Get()
G4int GetStrangeness() const
G4int GetSPDGCode() const
G4int GetZNSPDGCode() const
G4int AddParton(G4int pPDG) const
G4QHadronVector * ProjFragment(G4int pPDG, G4LorentzVector p4M, G4int tgZ, G4int tgN)
static G4QDiffractionRatio * GetPointer()
G4double GetRatio(G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
void DecayMeson(G4QHadron *dB, G4QHadronVector *HV)
void DecayAntistrange(G4QHadron *aS, G4QHadronVector *HV)
void AddQuasmon(G4Quasmon *Q)
void DecayBaryon(G4QHadron *dB, G4QHadronVector *HV)
G4QHadronVector * Fragment()
G4int SumPartonPDG(G4int PDG1, G4int PFG2) const
std::pair< G4int, G4int > ReducePair(G4int P1, G4int P2) const
static void SetParameters(G4int nC, G4double strTens, G4double tubeDens, G4double SigPt)
G4bool IsSingleDiffractive()
G4double ChooseX(G4double Xmin, G4double Xmax) const
G4QFragmentation(const G4QNucleus &aNucleus, const G4QHadron &aPrimary)
G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare) const
G4QHadronVector * Fragment()
G4bool ExciteDiffParticipants(G4QHadron *aPartner, G4QHadron *bPartner) const
G4bool ExciteSingDiffParticipants(G4QHadron *aPartner, G4QHadron *bPartner) const
void EvaporateResidual(G4QHadron *hadrNuc)
G4int AnnihilationOrder(G4int LS, G4int MS, G4int uP, G4int mP, G4int sP, G4int nP)
G4LorentzVector Get4Momentum() const
void SetNFragments(const G4int &nf)
G4int GetBaryonNumber() const
void IncrementCollisionCount(G4int aCount)
G4QParton * GetNextParton()
void SetPDGCode(const G4QPDGCode &PDG)
G4bool DecayIn2(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom)
void Boost(const G4LorentzVector &theBoost)
const G4ThreeVector & GetPosition() const
void SetQC(const G4QContent &newQC)
std::list< G4QParton * > GetColor()
void SetQPDG(const G4QPDGCode &QPDG)
std::list< G4QParton * > GetAntiColor()
G4QParton * GetNextAntiParton()
G4double GetMass2() const
G4int GetNFragments() const
void Set4Momentum(const G4LorentzVector &aMom)
G4QPDGCode GetQPDG() const
G4int GetNumberOfSoftCollisions()
void SetNumberOfSoftCollisions(G4int nofSoft)
G4QHadron * GetProjectile() const
G4QHadron * GetTarget() const
void SetTarget(G4QHadron *aTarget)
void SetNumberOfDINRCollisions(G4int nofDINR)
void SetNumberOfDiffractiveCollisions(G4int nofDiff)
void InitByPDG(G4int newPDG)
std::pair< G4double, G4double > ChooseImpactXandY(G4double maxImpact)
G4double GetOuterRadius()
void SubtractNucleon(G4QHadron *pNucleon)
void DoLorentzBoost(const G4LorentzVector &theBoost)
void DoLorentzRotation(const G4LorentzRotation &theLoRot)
G4QContent GetQCZNS() const
G4double GetThickness(G4double b)
G4LorentzVector GetNucleons4Momentum()
G4QHadron * GetNextNucleon()
G4double CoulombBarrier(const G4double &cZ=1, const G4double &cA=1, G4double dZ=0., G4double dA=0.)
void EvaporateNucleus(G4QHadron *hA, G4QHadronVector *oHV)
G4double GetGSMass() const
G4double GetNuclMass(G4int Z, G4int N, G4int S)
std::pair< G4int, G4int > MakeTwoBaryons(G4int L1, G4int L2, G4int R1, G4int R2)
G4double MinMassOfFragm()
void SetPDGCode(G4int aPDG)
const G4int & GetType() const
const G4LorentzVector & Get4Momentum() const
void Set4Momentum(const G4LorentzVector &aMomentum)
G4bool ReduceDiQADiQ(G4QParton *d1, G4QParton *d2)
G4double GetPomDiffProbability(const G4double s_value, const G4double imp2)
G4double GetCutPomProbability(const G4double s, const G4double ip2, const G4int nPom)
G4double GetPomInelProbability(const G4double s_value, const G4double imp2)
G4QHadronVector * FragmentString(G4bool QL)
G4QParton * GetLeftParton() const
G4int GetDirection() const
G4LorentzVector Get4Momentum() const
G4int GetBaryonNumber() const
G4QParton * GetRightParton() const
void Boost(G4ThreeVector &Velocity)
std::pair< G4LorentzVector, G4LorentzVector > Scatter(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
static G4QuasiFreeRatios * GetPointer()
std::pair< G4double, G4double > GetRatios(G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
void Set4Momentum(G4LorentzVector Q4M)
G4QHadronVector * DecayQHadron(G4QHadron *hadron)
void Boost(const G4LorentzVector &theBoost)
G4LorentzVector Get4Momentum() const
void SetQC(G4QContent QQC)
G4int GetBaryonNumber() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)