69{
71 static const G4double mProt2= mProt*mProt;
78 G4bool stringsInitted=
false;
82#ifdef edebug
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;
90#endif
93 G4int tPDG=90000000+tZ*1000+tN;
98#ifdef edebug
104 G4cout<<
"-EMC-G4QFragmentation::Construct:tLS4M="<<totLS4M<<
",tZLS4M="<<totZLS4M<<
G4endl;
105
106#endif
112 G4int ModelMode=SOFT;
115#ifdef debug
117#endif
119
120 std::pair<G4double,G4double> ratios=std::make_pair(0.,0.);
121 G4int apPDG=std::abs(pPDG);
122 if(apPDG>99)
123 {
125 ratios = theQuasiElastic->
GetRatios(pMom, pPDG, tZ, tN);
126 G4double qeRat = ratios.first*ratios.second;
128
129 if(qeRat<1. && proj4M.e()>thresh) difRat /= 1.-qeRat;
130 else difRat=0.;
131 difRat += qeRat;
133 if(rnd < qeRat)
134 {
139#ifdef debug
140 G4cout<<
">QE>G4QFragmentation::Construct:TryQE,R="<<ratios.first*ratios.second<<
",N4M="
141 <<pN4M<<
",NPDG="<<pNPDG<<
G4endl;
142#endif
143 std::pair<G4LorentzVector,G4LorentzVector> QEout=theQuasiElastic->
Scatter(pNPDG,pN4M,
144 pPDG,proj4M);
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)
150 {
152 theResult->push_back(qfN);
154 theResult->push_back(qfP);
159 theResult->push_back(resNuc);
160#ifdef debug
161 G4cout<<
"-->QE-->G4QFragmentation::Construct:QEDone, N4M="<<QEout.first<<
", p4M="
163#endif
164 return;
165 }
166 }
167 else if(rnd < difRat)
168 {
169#ifdef debug
170 G4cout<<
"-->Dif-->G4QFragmentation::Construct: qe="<<qeRat<<
", dif="<<difRat-qeRat
171 <<
",P="<<proj4M.vect().mag()<<
", tZ="<<tZ<<
", tN="<<tN<<
G4endl;
172#endif
174
175
176
178 G4int nSec=out->size();
179 if(nSec>1)
for(
G4int i=0; i<nSec; i++) theResult->push_back((*out)[i]);
180 else if(nSec>0)
181 {
182#ifdef debug
183 G4cout<<
"-Warning-G4QFragmentation::Construct: 1 secondary in Diffractionp 4M="
184 <<proj4M<<
", s4M="<<(*out)[0]->Get4Momentum()<<
G4endl;
185#endif
186 delete (*out)[0];
187 }
188 out->clear();
189 delete out;
190 if(nSec>1) return;
191 }
192 }
193#ifdef edebug
195 G4cout<<
"-EMC-G4QFragmentation::Construct: Nuc4M="<<sum1<<
G4endl;
196#endif
198
201 if(projAbsB>1) nCons = projAbsB;
202
204 G4double pz_per_projectile = proj4M.
pz()/nCons;
205
206 G4double e_per_projectile = proj4M.e()/nCons+mProt;
207 G4double vz = pz_per_projectile/e_per_projectile;
208#ifdef debug
209 G4cout<<
"G4QFragmentation::Construct: Projectile4M="<<proj4M<<
", Vz="<<vz<<
", nC="
210 <<nCons<<
", pE="<<e_per_projectile<<
G4endl;
211#endif
212 theCurrentVelocity.setZ(vz);
214#ifdef edebug
216 G4cout<<
"-EMC-G4QFragmentation::Construct: AfterBoost, v="<<vz<<
", Nuc4M="<<sum2<<
G4endl;
217#endif
219 cmProjMom.
boost(-theCurrentVelocity);
221#ifdef debug
223#endif
225 if(kE > 720.) maxCt=static_cast<int>(std::log(kE/270.));
226
227
228 G4int maxCuts=std::min( 7 , std::max(1, maxCt) );
229#ifdef debug
230 G4cout<<
"G4QFragmentation::Construct: Proj4MInCM="<<cmProjMom<<
", pPDG="<<pPDG<<
G4endl;
231#endif
232
233
234
236
239#ifdef debug
240 G4cout<<
"G4QFrag::Constr:OutR="<<outerRadius<<
",mC="<<maxCuts<<
",A="<<theNucleus<<
G4endl;
241#endif
243
246 G4double s_value = (cmProjMom + pNuc4M).mag2();
248#ifdef debug
249 G4cout<<
"G4QFrag::Construc: p4M="<<cmProjMom<<
", tgN4M="<<pNuc4M<<
", s="<<s_value<<
", ThreshM="
251#endif
252 ModelMode = SOFT;
253 if (s_value < 0.)
254 {
255 G4cerr<<
"***G4QFragmentation::Construct: s="<<s_value<<
", pN4M="<<pNuc4M<<
G4endl;
257 }
258 else if(s_value < mProt2)
259 {
266 {
268 G4double cs=(cmProjMom + cp4M).mag2();
269 if(cs > maxS)
270 {
271 maxS=cs;
272 pNucleon=bNucleon;
273 }
274 }
276 pProjectile =cmProjectile;
282 delete aTarget;
283 delete pProjectile;
284
285
286
287
288
289
290
291
292
293
294
295
296
297 Q4M.
boost(theCurrentVelocity);
298 Q4M=toLab*Q4M;
300 theQuasmons.push_back(stringQuasmon);
303 return;
304 }
305 if (s_value <
sqr(ThresholdMass))
306 {
307#ifdef debug
308 G4cout<<
"G4QFragmentation::Construct:*OnlyDiffraction*ThM="<<ThresholdMass<<
">sqrt(s)="
309 <<std::sqrt(s_value)<<
" -> only Diffraction is possible"<<
G4endl;
310#endif
311 ModelMode = DIFFRACTIVE;
312 }
313
314#ifdef debug
315 G4cout<<
"G4QFragmentation::Construct: theIntSize="<<theInteractions.size()<<
G4endl;
316#endif
318 theInteractions.clear();
321
325 if (proB>0) prEn-=aProjectile.
GetMass();
326 else if(proB<0) prEn+=mProt;
327#ifdef debug
329 G4cout<<
"G4QFragmentation::Construct: estimated energy, prEn="<<prEn<<
G4endl;
330#endif
331 while(!theInteractions.size() && ++attCnt < maxAtt)
332 {
333#ifdef debug
334 G4cout<<
"G4QFragmentation::Construct: *EnterTheInteractionLOOP*, att#"<<attCnt<<
G4endl;
335#endif
336
337 std::pair<G4double, G4double> theImpactParameter;
339 G4double impactX = theImpactParameter.first;
340 G4double impactY = theImpactParameter.second;
341#ifdef debug
342 G4cout<<
"G4QFragmentation::Construct: Impact Par X="<<impactX<<
", Y="<<impactY<<
G4endl;
343#endif
344 G4double impar=std::sqrt(impactX*impactX+impactY*impactY);
347 maxEn=eflen*stringTension;
348 maxNuc=static_cast<int>(eflen*tubeDensity+0.5);
349#ifdef debug
350 G4cout<<
"G4QFragment::Construct: pE="<<prEn<<
" <? mE="<<maxEn<<
", mN="<<maxNuc<<
G4endl;
351#endif
352 if(prEn < maxEn)
353 {
360 theInteractions.push_back(anInteraction);
364#ifdef debug
365 G4cout<<
"G4QFragmentation::Construct: DINR interaction is created"<<
G4endl;
366#endif
367 break;
368 }
369
371 G4int nucleonCount = 0;
372#ifdef debug
373 G4cout<<
"G4QFragment::Construct:BeforeWhileOveNuc, A="<<nA<<
",p4M="<<cmProjMom<<
G4endl;
374#endif
375 while( (pNucleon=theNucleus.
GetNextNucleon()) && nucleonCount<nA && totalCuts<maxCuts)
376 {
377 ++nucleonCount;
378
379 s_value = (cmProjMom + pNucleon->
Get4Momentum()).mag2();
380#ifdef debug
381 G4cout<<
"G4QFrag::Constr:N# "<<nucleonCount<<
", s="<<s_value<<
", tgN4M="
383#endif
384 if(s_value<=10000.)
385 {
386#ifdef debug
387 G4cout<<
"G4QFragmentation::Construct: SKIP, s<.01 GeV^2, p4M="<<cmProjMom
389#endif
390 continue;
391 }
392#ifdef sdebug
393 G4cout<<
"G4QFragmentation::Construct:LOOPovNuc,nC="<<nucleonCount<<
", s="<<s_value<<
G4endl;
395#endif
398#ifdef sdebug
399 G4cout<<
"G4QFragmentation::Construct: s="<<s_value<<
", D2="<<Distance2<<
G4endl;
400#endif
401 G4double Probability = theProbability.GetPomInelProbability(s_value, Distance2);
402
403#ifdef sdebug
404 G4cout<<
"G4QFragmentation::Construct: Probubility="<<Probability<<
G4endl;
405#endif
407
408#ifdef sdebug
409 G4cout<<
"G4QFragmentation::Construct: NLOOP prob="<<Probability<<
", rndm="<<rndNumber
410 <<
", d="<<std::sqrt(Distance2)<<
G4endl;
411#endif
412 if (Probability > rndNumber)
413 {
415#ifdef edebug
416 G4cout<<
"--->EMC-->G4QFragmentation::Construct: Target Nucleon is filled, 4M/PDG="
418#endif
419
423 if((theProbability.GetPomDiffProbability(s_value,Distance2)/Probability >
424 G4UniformRand() && ModelMode==SOFT ) || ModelMode==DIFFRACTIVE)
425 {
426
432 theInteractions.push_back(anInteraction);
433
434 totalCuts++;
435#ifdef debug
436 G4cout<<
"G4QFragmentation::Construct:NLOOP DiffInteract, tC="<<totalCuts<<
G4endl;
437#endif
438 }
439 else
440 {
441
442
445 for(nCut = 0; nCut < nCutMax; nCut++)
446 {
447 running[nCut]= theProbability.GetCutPomProbability(s_value, Distance2, nCut+1);
448 if(nCut) running[nCut] += running[nCut-1];
449 }
451 for(nCut = 0; nCut < nCutMax; nCut++) if(running[nCut] > random) break;
452 delete [] running;
453#ifdef debug
454 G4cout<<
"G4QFragmentation::Construct: NLOOP-Soft Chosen nCut="<<nCut<<
G4endl;
455#endif
456
457
458
464 theInteractions.push_back(anInteraction);
465 totalCuts += nCut+1;
466#ifdef debug
467 G4cout<<
"G4QFragmentation::Construct:NLOOP SoftInteract, tC="<<totalCuts<<
G4endl;
468 impactUsed=Distance2;
469#endif
470 }
471 }
472 }
473
474#ifdef debug
475 G4cout<<
"G4QFragmentation::Construct: NUCLEONCOUNT="<<nucleonCount<<
G4endl;
476#endif
477 }
478 G4int nInt=theInteractions.size();
479#ifdef debug
480 G4cout<<
"G4QFrag::Con:CUT="<<totalCuts<<
",ImpPar="<<impactUsed<<
",#ofInt="<<nInt<<
G4endl;
481#endif
482
483 if(!nInt || (nInt==1 && theInteractions[0]->GetNumberOfDINRCollisions()==1))
484 {
487 if(nInt)
488 {
489 aTarget=theInteractions[0]->GetTarget();
490 pProjectile=theInteractions[0]->GetProjectile();
491 delete theInteractions[0];
492 theInteractions.clear();
493 }
494 else
495 {
499 pProjectile =cmProjectile;
503 }
506 delete aTarget;
507 delete pProjectile;
508
509
510
511
512
513
514
515
516
517
518
519
520
521 Q4M.
boost(theCurrentVelocity);
522 Q4M=toLab*Q4M;
524 theQuasmons.push_back(stringQuasmon);
527 return;
528 }
529
530
531
532#ifdef debug
533 G4cout<<
"G4QFragmentation::Construct: Before PartPairCreation nInt="<<nInt<<
G4endl;
534#endif
535 for(
G4int i=0; i<nInt; i++)
536 {
537 theInteractions[i]->SplitHadrons();
538#ifdef edebug
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)
556 {
557 pSP+=(*picp)->Get4Momentum();
558 pSP+=(*piac)->Get4Momentum();
559 tSP+=(*ticp)->Get4Momentum();
560 tSP+=(*tiac)->Get4Momentum();
561 }
562 G4cout<<
"-EMC-G4QFragmentation::Construct: Interaction#"<<i<<
",dP4M="
564#endif
565 }
566
567
568
569 G4QInteractionVector::iterator it;
570#ifdef debug
571 G4cout<<
"G4QFragmentation::Construct: Creation ofSoftCollisionPartonPair STARTS"<<
G4endl;
572#endif
574 while(rep && theInteractions.size())
575 {
576 for(it = theInteractions.begin(); it != theInteractions.end(); ++it)
577 {
581#ifdef debug
582 G4cout<<
"G4QFragmentation::Construct: #0f SOFT collisions ="<<nSoftCollisions<<
G4endl;
583#endif
584 if (nSoftCollisions)
585 {
588 for (
G4int j = 0; j < nSoftCollisions; j++)
589 {
593 thePartonPairs.push_back(aPair);
597 thePartonPairs.push_back(aPair);
598#ifdef debug
599 G4cout<<
"--->G4QFragmentation::Construct: SOFT, 2 parton pairs are filled"<<
G4endl;
600#endif
601 }
602 delete *it;
603 it=theInteractions.erase(it);
604 if( it != theInteractions.begin() )
605 {
606 it--;
607 rep=false;
608#ifdef debug
609 G4cout<<
"G4QFragmentation::Construct: *** Decremented ***"<<
G4endl;
610#endif
611 }
612 else
613 {
614 rep=true;
615#ifdef debug
616 G4cout<<
"G4QFragmentation::Construct: *** IT Begin ***"<<
G4endl;
617#endif
618 break;
619 }
620 }
621 else rep=false;
622#ifdef debug
623 G4cout<<
"G4QFragmentation::Construct: #0fSC="<<nSoftCollisions<<
", r="<<rep<<
G4endl;
624#endif
625 }
626#ifdef debug
627 G4cout<<
"G4QFragmentation::Construct: *** IT While *** , r="<<rep<<
G4endl;
628#endif
629 }
630#ifdef debug
631 G4cout<<
"G4QFragmentation::Construct: -> Parton pairs for SOFT strings are made"<<
G4endl;
632#endif
633
634
635
636 for(unsigned i = 0; i < theInteractions.size(); i++)
637 {
638
641#ifdef debug
642 G4cout<<
"G4QFragmentation::Construct: CreationOfDiffractivePartonPairs, i="<<i<<
G4endl;
643#endif
644
647 if (aParton)
648 {
652 thePartonPairs.push_back(aPartonPair);
653#ifdef debug
654 G4cout<<
"G4QFragmentation::Construct: proj Diffractive PartonPair is filled"<<
G4endl;
655#endif
656 }
657
660 if (aParton)
661 {
664 thePartonPairs.push_back(aPartonPair);
665#ifdef debug
666 G4cout<<
"G4QFragmentation::Constr: target Diffractive PartonPair is filled"<<
G4endl;
667#endif
668 }
669 }
670#ifdef debug
671 G4cout<<
"G4QFragmentation::Construct: DiffractivePartonPairs are created"<<
G4endl;
672#endif
673
674
675
677 theInteractions.clear();
678 delete cmProjectile;
679#ifdef debug
680 G4cout<<
"G4QFragmentation::Construct: Temporary objects are cleaned up"<<
G4endl;
681#endif
682
684
685#ifdef debug
686 G4cout<<
"--------->>G4QFragmentation::Construct: ------->> Strings are created "<<
G4endl;
687#endif
690 while(thePartonPairs.size())
691 {
692 aPair = thePartonPairs.back();
693 thePartonPairs.pop_back();
694#ifdef debug
696#endif
698#ifdef debug
700#endif
701 aString->
Boost(theCurrentVelocity);
702 strings.push_back(aString);
703 stringsInitted=true;
704 delete aPair;
705 }
706#ifdef edebug
710 G4int nStrings=strings.size();
711 G4cout<<
"-EMC-G4QFragmentation::Construct:#ofString="<<nStrings<<
",tNuc4M="<<sum<<
G4endl;
712 for(
G4int i=0; i<nStrings; i++)
713 {
716 sum+=strI4M;
723 rChg-=sChg;
724 rBaN-=sBaN;
725 G4cout<<
"-EMC-G4QFragmentation::Construct: String#"<<i<<
", 4M="<<strI4M<<
",LPDG="<<LPDG
726 <<LQC<<
",RPDG="<<RPDG<<RQC<<
", Ch="<<sChg<<
", BN="<<sBaN<<
G4endl;
727 }
728 G4cout<<
"-EMC-G4QFragm::Constr: r4M="<<sum-totZLS4M<<
",rC="<<rChg<<
",rB="<<rBaN<<
G4endl;
729#endif
730 if(!stringsInitted)
731 {
732 G4cerr<<
"******G4QFragmentation::Construct:***** No strings are created *****"<<
G4endl;
734 }
735#ifdef debug
736 G4cout<<
"G4QFragmentation::Constr: BeforeRotation, #0fStrings="<<strings.size()<<
G4endl;
737#endif
738
739
740
741 for(unsigned astring=0; astring < strings.size(); astring++)
742 strings[astring]->LorentzRotate(toLab);
744
745#ifdef edebug
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++)
752 {
755 G4int sChg=strings[i]->GetCharge();
756 rCg-=sChg;
757 G4int sBaN=strings[i]->GetBaryonNumber();
758 rBC-=sBaN;
759 G4cout<<
"-EMCLS-G4QFragm::Construct:String#"<<i<<
",4M="<<strI4M<<strI4M.
m()<<
",Charge="
760 <<sChg<<
",BaryN="<<sBaN<<
G4endl;
761 }
762 G4cout<<
"-EMCLS-...G4QFragm::Constr:r4M="<<
sm-totLS4M<<
",rC="<<rCg<<
",rB="<<rBC<<
G4endl;
763#endif
764
765
766
768#ifdef edebug
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++)
775 {
778 G4int sChg=strings[i]->GetCharge();
779 rCg-=sChg;
780 G4int sBaN=strings[i]->GetBaryonNumber();
781 rBC-=sBaN;
782 G4cout<<
"-EMCLS-G4QFragm::Construct:String#"<<i<<
",4M="<<strI4M<<strI4M.
m()<<
",Charge="
783 <<sChg<<
",BaryN="<<sBaN<<
G4endl;
784 }
785 G4cout<<
"-EMCLS-...G4QFragm::Constr:r4M="<<
sm-totLS4M<<
",rC="<<rCg<<
",rB="<<rBC<<
G4endl;
786#endif
787
788
789
791 G4QStringVector::iterator ist;
793 while(con && strings.size())
794 {
795 for(ist = strings.begin(); ist < strings.end(); ++ist)
796 {
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;
814 if(aLPDG)
815 {
816 L1=aLPDG/10;
817 L2=aLPDG%10;
818 }
821 if(aRPDG)
822 {
823 R1=aRPDG/10;
824 R2=aRPDG%10;
825 }
827 if(cSM2>0.) cSM=std::sqrt(cSM2);
828#ifdef debug
829 G4cout<<
"G4QFrag::Constr:NeedsFusion? cLPDG="<<cLPDG<<
",cRPDG="<<cRPDG<<
",cM(cM2If<0)="
830 <<cSM<<
",c4M"<<cS4M<<
G4endl;
831#endif
832 if(cSM>0.)
833 {
836 if(cLT==2 && cRT==2)
837 {
838 if(L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2)
839 {
840 single=false;
844 }
845 }
847
848#ifdef debug
849 G4cout<<
"G4QFrag::Const:*IsItGood? realM="<<std::sqrt(cSM2)<<
" > GSM="<<miM<<
G4endl;
850#endif
851 if(std::sqrt(cSM2) > miM) bad=false;
852 }
853 if(bad)
854 {
855#ifdef debug
856 G4cout<<
"G4QFrag::Const:TryFuse,L1="<<L1<<
",L2="<<L2<<
",R1="<<R1<<
",R2="<<R2<<
G4endl;
857#endif
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)
869 {
876#ifdef debug
877 G4cout<<
"->G4QFragm::Construct: sum4M="<<pS4M<<
",M2="<<pSM2<<
",p4M="<<sS4M<<
G4endl;
878#endif
879
880
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;
895 if(LPDG)
896 {
897 pL1=LPDG/10;
898 pL2=LPDG%10;
899 }
902 if(RPDG)
903 {
904 pR1=RPDG/10;
905 pR2=RPDG%10;
906 }
908#ifdef debug
909 G4cout<<
"G4QFragm::Construct: Partner/w pLPDG="<<pLPDG<<
", pRPDG="<<pRPDG<<
", pM2="
911#endif
912
915 if (cST==2 && pST==2)
916 {
917 tf=1;
918 af=1;
919 }
920 else if(cST==2 && pST==3)
921 {
922 tf=2;
923 if (pLPDG > 7 &&
924 ( (cLPDG<0 && (-cLPDG==pL1 || -cLPDG==pL2 || -cLPDG==pRPDG) ) ||
925 (cRPDG<0 && (-cRPDG==pL1 || -cRPDG==pL2 || -cRPDG==pRPDG) )
926 )
927 ) af=1;
928 else if(pRPDG > 7 &&
929 ( (cLPDG<0 && (-cLPDG==pR1 || -cLPDG==pR2 || -cLPDG==pLPDG) ) ||
930 (cRPDG<0 && (-cRPDG==pR1 || -cRPDG==pR2 || -cRPDG==pLPDG) )
931 )
932 ) af=2;
933 else if(pLPDG <-7 &&
934 ( (cLPDG>0 && ( cLPDG==pL1 || cLPDG==pL2 || cLPDG==-pRPDG) ) ||
935 (cRPDG>0 && ( cRPDG==pL1 || cRPDG==pL2 || cRPDG==-pRPDG) )
936 )
937 ) af=3;
938 else if(pRPDG <-7 &&
939 ( (cLPDG>0 && ( cLPDG==pR1 || cLPDG==pR2 || cLPDG==-pLPDG) ) ||
940 (cRPDG>0 && ( cRPDG==pR1 || cRPDG==pR2 || cRPDG==-pLPDG) )
941 )
942 ) af=4;
943#ifdef debug
944 else G4cout<<
"G4QFragmentation::Construct:2(QaQ+QDiQ/aQaDiQ) Can't fuse"<<
G4endl;
945#endif
946 }
947 else if(cST==3 && pST==2)
948 {
949 tf=3;
950 if (cLPDG > 7 &&
951 ( (pLPDG<0 && (-pLPDG==L1 || -pLPDG==L2 || -pLPDG==cRPDG) ) ||
952 (pRPDG<0 && (-pRPDG==L1 || -pRPDG==L2 || -pRPDG==cRPDG) )
953 )
954 ) af=1;
955 else if(cRPDG > 7 &&
956 ( (pLPDG<0 && (-pLPDG==R1 || -pLPDG==R2 || -pLPDG==cLPDG) ) ||
957 (pRPDG<0 && (-pRPDG==R1 || -pRPDG==R2 || -pRPDG==cLPDG) )
958 )
959 ) af=2;
960 else if(cLPDG <-7 &&
961 ( (pLPDG>0 && ( pLPDG==L1 || pLPDG==L2 || pLPDG==-cRPDG) ) ||
962 (pRPDG>0 && ( pRPDG==L1 || pRPDG==L2 || pRPDG==-cRPDG) )
963 )
964 ) af=3;
965 else if(cRPDG <-7 &&
966 ( (pLPDG>0 && ( pLPDG==R1 || pLPDG==R2 || pLPDG==-cLPDG) ) ||
967 (pRPDG>0 && ( pRPDG==R1 || pRPDG==R2 || pRPDG==-cLPDG) )
968 )
969 ) af=4;
970#ifdef debug
971 else G4cout<<
"G4QFragmentation::Construct:3(QDiQ/aQaDiQ+QaQ) Can't fuse"<<
G4endl;
972#endif
973 }
974 else if(cST==2 && pST==4)
975 {
976 tf=4;
977 if (pLPDG > 7)
978 {
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;
981 }
982 else if(pRPDG > 7)
983 {
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;
986 }
987#ifdef debug
988 else G4cout<<
"-G4QFragmentation::Construct: 4 (QaQ+aQDiQDiQ) Can't fuse"<<
G4endl;
989#endif
990 }
991 else if(cST==4 && pST==2)
992 {
993 tf=5;
994 if (cLPDG > 7)
995 {
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;
998 }
999 else if(cRPDG > 7)
1000 {
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;
1003 }
1004#ifdef debug
1005 else G4cout<<
"-G4QFragmentation::Construct: 5 (aQDiQDiQ+QaQ) Can't fuse"<<
G4endl;
1006#endif
1007 }
1008 else if(cST==3 && pST==3)
1009 {
1010 tf=6;
1011 if(pLPDG > 7)
1012 {
1013 if (cLPDG<-7 && (-cRPDG==pL1 || -cRPDG==pL2) && (pRPDG==L1 || pRPDG==L2))
1014 af=1;
1015 else if(cRPDG<-7 && (-cLPDG==pL1 || -cLPDG==pL2) && (pRPDG==R1 || pRPDG==R2))
1016 af=2;
1017 }
1018 else if(pRPDG > 7)
1019 {
1020 if (cLPDG<-7 && (-cRPDG==pR1 || -cRPDG==pR2) && (pLPDG==L1 || pLPDG==L2))
1021 af=3;
1022 else if(cRPDG<-7 && (-cLPDG==pR1 || -cLPDG==pR2) && (pLPDG==R1 || pLPDG==R2))
1023 af=4;
1024 }
1025 else if(cLPDG > 7)
1026 {
1027 if (pLPDG<-7 && (-pRPDG==L1 || -pRPDG==L2) && (cRPDG==pL1 || cRPDG==pL2))
1028 af=5;
1029 else if(pRPDG<-7 && (-pLPDG==L1 || -pLPDG==L2) && (cRPDG==pR1 || cRPDG==pR2))
1030 af=6;
1031 }
1032 else if(cRPDG > 7)
1033 {
1034 if (pLPDG<-7 && (-pRPDG==R1 || -pRPDG==R2) && (cLPDG==pL1 || cLPDG==pL2))
1035 af=7;
1036 else if(pRPDG<-7 && (-pLPDG==R1 || -pLPDG==R2) && (cLPDG==pR1 || cLPDG==pR2))
1037 af=8;
1038 }
1039#ifdef debug
1040 else G4cout<<
"-G4QFragmentation::Construct: 6 (QDiQ+aQaDiQ) Can't fuse"<<
G4endl;
1041#endif
1042 }
1043#ifdef debug
1044 G4cout<<
"G4QFragmentation::Const: ***Possibility***, tf="<<tf<<
", af="<<af<<
G4endl;
1045#endif
1046 if(tf && af)
1047 {
1048
1050 switch (tf)
1051 {
1052 case 1:
1053 if (cLPDG > 0 && pLPDG > 0)
1054 {
1055 order= 1;
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;
1062 }
1063 else if(cLPDG < 0 && pLPDG < 0)
1064 {
1065 order= 1;
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;
1072 }
1073 else if(cRPDG > 0 && pLPDG > 0)
1074 {
1075 order=-1;
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;
1082 }
1083 else if(cRPDG < 0 && pLPDG < 0)
1084 {
1085 order=-1;
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;
1092 }
1093 break;
1094 case 2:
1095 switch (af)
1096 {
1097 case 1:
1098 if(cLPDG < 0)
1099 {
1100 order= 1;
1101 if(-cLPDG==pRPDG)
1102 {
1103 nLPDG=pLPDG;
1104 nRPDG=cRPDG;
1105 }
1106 else
1107 {
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;
1112 else nLPDG=pL1;
1113 }
1114 }
1115 else
1116 {
1117 order=-1;
1118 if(-cRPDG==pRPDG)
1119 {
1120 nLPDG=pLPDG;
1121 nRPDG=cLPDG;
1122 }
1123 else
1124 {
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;
1129 else nLPDG=pL1;
1130 }
1131 }
1132 break;
1133 case 2:
1134 if(cLPDG < 0)
1135 {
1136 order=-1;
1137 if(-cLPDG==pLPDG)
1138 {
1139 nLPDG=cRPDG;
1140 nRPDG=pRPDG;
1141 }
1142 else
1143 {
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;
1148 else nRPDG=pR1;
1149 }
1150 }
1151 else
1152 {
1153 order= 1;
1154 if(-cRPDG==pLPDG)
1155 {
1156 nLPDG=cLPDG;
1157 nRPDG=pRPDG;
1158 }
1159 else
1160 {
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;
1165 else nRPDG=pR1;
1166 }
1167 }
1168 break;
1169 case 3:
1170 if(cLPDG > 0)
1171 {
1172 order= 1;
1173 if(cLPDG==-pRPDG)
1174 {
1175 nLPDG=pLPDG;
1176 nRPDG=cRPDG;
1177 }
1178 else
1179 {
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;
1184 else nLPDG=-pL1;
1185 }
1186 }
1187 else
1188 {
1189 order=-1;
1190 if(cRPDG==-pRPDG)
1191 {
1192 nLPDG=pLPDG;
1193 nRPDG=cLPDG;
1194 }
1195 else
1196 {
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;
1201 else nLPDG=-pL1;
1202 }
1203 }
1204 break;
1205 case 4:
1206 if(cLPDG > 0)
1207 {
1208 order=-1;
1209 if(cLPDG==-pLPDG)
1210 {
1211 nLPDG=cRPDG;
1212 nRPDG=pRPDG;
1213 }
1214 else
1215 {
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;
1220 else nRPDG=-pR1;
1221 }
1222 }
1223 else
1224 {
1225 order= 1;
1226 if(cRPDG==-pLPDG)
1227 {
1228 nLPDG=cLPDG;
1229 nRPDG=pRPDG;
1230 }
1231 else
1232 {
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;
1237 else nRPDG=-pR1;
1238 }
1239 }
1240 break;
1241 }
1242 break;
1243 case 3:
1244 switch (af)
1245 {
1246 case 1:
1247 if(pLPDG < 0)
1248 {
1249 order= 1;
1250 if(-pLPDG==cRPDG)
1251 {
1252 nLPDG=cLPDG;
1253 nRPDG=pRPDG;
1254 }
1255 else
1256 {
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;
1261 else nLPDG=L1;
1262 }
1263 }
1264 else
1265 {
1266 order=-1;
1267 if(-pRPDG==cRPDG)
1268 {
1269 nLPDG=cLPDG;
1270 nRPDG=pLPDG;
1271 }
1272 else
1273 {
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;
1278 else nRPDG=L1;
1279 }
1280 }
1281 break;
1282 case 2:
1283 if(pLPDG < 0)
1284 {
1285 order=-1;
1286 if(-pLPDG==cLPDG)
1287 {
1288 nLPDG=pRPDG;
1289 nRPDG=cRPDG;
1290 }
1291 else
1292 {
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;
1297 else nLPDG=R1;
1298 }
1299 }
1300 else
1301 {
1302 order= 1;
1303 if(-pRPDG==cLPDG)
1304 {
1305 nLPDG=pLPDG;
1306 nRPDG=cRPDG;
1307 }
1308 else
1309 {
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;
1314 else nRPDG=R1;
1315 }
1316 }
1317 break;
1318 case 3:
1319 if(pLPDG > 0)
1320 {
1321 order= 1;
1322 if(pLPDG==-cRPDG)
1323 {
1324 nLPDG=cLPDG;
1325 nRPDG=pRPDG;
1326 }
1327 else
1328 {
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;
1333 else nLPDG=-L1;
1334 }
1335 }
1336 else
1337 {
1338 order=-1;
1339 if(pRPDG==-cRPDG)
1340 {
1341 nLPDG=cLPDG;
1342 nRPDG=pLPDG;
1343 }
1344 else
1345 {
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;
1350 else nRPDG=-L1;
1351 }
1352 }
1353 break;
1354 case 4:
1355 if(pLPDG > 0)
1356 {
1357 order=-1;
1358 if(pLPDG==-cLPDG)
1359 {
1360 nLPDG=pRPDG;
1361 nRPDG=cRPDG;
1362 }
1363 else
1364 {
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;
1369 else nLPDG=-R1;
1370 }
1371 }
1372 else
1373 {
1374 order= 1;
1375 if(pRPDG==-cLPDG)
1376 {
1377 nLPDG=pLPDG;
1378 nRPDG=cRPDG;
1379 }
1380 else
1381 {
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;
1386 else nRPDG=-R1;
1387 }
1388 }
1389 break;
1390 }
1391 break;
1392 case 4:
1393 switch (af)
1394 {
1395 case 1:
1396 order= 1;
1397 if(-cLPDG == pL1) nLPDG= pL2;
1398 else nLPDG= pL1;
1399 if( cRPDG == pR1) nRPDG=-pR2;
1400 else nRPDG=-pR1;
1401 break;
1402 case 2:
1403 order=-1;
1404 if(-cRPDG == pL1) nLPDG= pL2;
1405 else nLPDG= pL1;
1406 if( cLPDG == pR1) nRPDG=-pR2;
1407 else nRPDG=-pR1;
1408 break;
1409 case 3:
1410 order= 1;
1411 if( cLPDG == pL1) nLPDG=-pL2;
1412 else nLPDG=-pL1;
1413 if(-cRPDG == pR1) nRPDG= pR2;
1414 else nRPDG= pR1;
1415 break;
1416 case 4:
1417 order=-1;
1418 if( cRPDG == pL1) nLPDG=-pL2;
1419 else nLPDG=-pL1;
1420 if(-cLPDG == pR1) nRPDG= pR2;
1421 else nRPDG= pR1;
1422 break;
1423 }
1424 break;
1425 case 5:
1426 switch (af)
1427 {
1428 case 1:
1429 order= 1;
1430 if(-pLPDG == L1) nLPDG= L2;
1431 else nLPDG= L1;
1432 if( pRPDG == R1) nRPDG=-R2;
1433 else nRPDG=-R1;
1434 break;
1435 case 2:
1436 order=-1;
1437 if(-pRPDG == L1) nRPDG= L2;
1438 else nRPDG= L1;
1439 if( pLPDG == R1) nLPDG=-R2;
1440 else nLPDG=-R1;
1441 break;
1442 case 3:
1443 order= 1;
1444 if( pLPDG == L1) nLPDG=-L2;
1445 else nLPDG=-L1;
1446 if(-pRPDG == R1) nRPDG= R2;
1447 else nRPDG= R1;
1448 break;
1449 case 4:
1450 order=-1;
1451 if( pRPDG == L1) nRPDG=-L2;
1452 else nRPDG=-L1;
1453 if(-pLPDG == R1) nLPDG= R2;
1454 else nLPDG= R1;
1455 break;
1456 }
1457 break;
1458 case 6:
1459 switch (af)
1460 {
1461 case 1:
1462 order=-1;
1463 if(-cRPDG == pL1) nLPDG= pL2;
1464 else nLPDG= pL1;
1465 if( pRPDG == L1) nRPDG= -L2;
1466 else nRPDG= -L1;
1467 break;
1468 case 2:
1469 order= 1;
1470 if(-cLPDG == pL1) nLPDG= pL2;
1471 else nLPDG= pL1;
1472 if( pRPDG == R1) nRPDG= -R2;
1473 else nRPDG= -R1;
1474 break;
1475 case 3:
1476 order= 1;
1477 if(-cRPDG == pR1) nRPDG= pR2;
1478 else nRPDG= pR1;
1479 if( pLPDG == L1) nLPDG= -L2;
1480 else nLPDG= -L1;
1481 break;
1482 case 4:
1483 order=-1;
1484 if(-cLPDG == pR1) nRPDG= pR2;
1485 else nRPDG= pR1;
1486 if( pLPDG == R1) nLPDG= -R2;
1487 else nLPDG= -R1;
1488 break;
1489 case 5:
1490 order=-1;
1491 if(-pRPDG == L1) nRPDG= L2;
1492 else nRPDG= L1;
1493 if( cRPDG == pL1) nLPDG=-pL2;
1494 else nLPDG=-pL1;
1495 break;
1496 case 6:
1497 order= 1;
1498 if(-pLPDG == L1) nLPDG= L2;
1499 else nLPDG= L1;
1500 if( cRPDG == pR1) nRPDG=-pR2;
1501 else nRPDG=-pR1;
1502 break;
1503 case 7:
1504 order= 1;
1505 if(-pRPDG == R1) nRPDG= R2;
1506 else nRPDG= R1;
1507 if( cLPDG == pL1) nLPDG=-pL2;
1508 else nLPDG=-pL1;
1509 break;
1510 case 8:
1511 order=-1;
1512 if(-pLPDG == R1) nLPDG= R2;
1513 else nLPDG= R1;
1514 if( cLPDG == pR1) nRPDG=-pR2;
1515 else nRPDG=-pR1;
1516 break;
1517 }
1518 break;
1519 }
1520 if(!order)
G4cerr<<
"-Warning-G4QFrag::Constr: t="<<tf<<
", a="<<af<<
", cL="<<cLPDG
1521 <<
", cR="<<cRPDG<<
", pL="<<pLPDG<<
", pR="<<pRPDG<<
G4endl;
1522 else
1523 {
1524
1526 if(std::abs(nLPDG) > 7) ++
LT;
1528 if(std::abs(nRPDG) > 7) ++RT;
1531 if(cLT==2 && cRT==2)
1532 {
1533 aLPDG=0;
1534 aRPDG=0;
1535 if(cLPDG>0)
1536 {
1537 aLPDG=nLPDG/100;
1538 aRPDG=(-nRPDG)/100;
1539 }
1540 else
1541 {
1542 aRPDG=nRPDG/100;
1543 aLPDG=(-nLPDG)/100;
1544 }
1549 if(nL1!=nR1 && nL1!=nR2 && nL2!=nR1 && nL2!=nR2)
1550 {
1551#ifdef debug
1552 G4cout<<
"G4QFragmentation::Const:aLPDG="<<aLPDG<<
", aRPDG="<<aRPDG<<
G4endl;
1553#endif
1554 sing=false;
1556 std::pair<G4int,G4int> pB=tmp.
MakeTwoBaryons(nL1, nL2, nR1, nR2);
1558 }
1559 }
1560 if(sing)
1561 {
1562 std::pair<G4int,G4int> newPair = std::make_pair(nLPDG,nRPDG);
1564#ifdef debug
1565 G4cout<<
"G4QFr::Con: LPDG="<<nLPDG<<
",RPDG="<<nRPDG<<
",QC="<<newStQC<<
G4endl;
1566#endif
1567 G4int minPDG=newStQC.GetSPDGCode();
1569 }
1570
1573 if(pSM2>0.) pSM=std::sqrt(pSM2);
1574 if(minC && pSM2 > maxiM2)
1575 {
1576 maxiM2=pSM2;
1577 win=true;
1578 }
1579 else if(!minC || pSM > minM)
1580 {
1581 pExcess=pSM-minM;
1582 if(minC || pExcess > excess)
1583 {
1584 minC=false;
1585 excess=pExcess;
1586 win=true;
1587 }
1588 }
1589 if(win)
1590 {
1591 sst=pst;
1592 sLPDG=nLPDG;
1593 sRPDG=nRPDG;
1594 sOrd=order;
1595 }
1596 }
1597 }
1598
1599 }
1600 if(sOrd)
1601 {
1604 G4QParton* pLeft=(*sst)->GetLeftParton();
1605 G4QParton* pRight=(*sst)->GetRightParton();
1608#ifdef debug
1609 G4cout<<
"G4QFragmentation::Const:cS4M="<<cS4M<<
" fused/w pS4M="<<pL4M+pR4M<<
G4endl;
1610#endif
1611 if(sOrd>0)
1612 {
1613 pL4M+=cL4M;
1614 pR4M+=cR4M;
1615 }
1616 else
1617 {
1618 pL4M+=cR4M;
1619 pR4M+=cL4M;
1620 }
1625 delete (*ist);
1626 strings.erase(ist);
1627#ifdef debug
1629 G4cout<<
"G4QFragmentation::Construct:Created,4M="<<ss4M<<
",m2="<<ss4M.
m2()<<
G4endl;
1630#endif
1631 if( ist != strings.begin() )
1632 {
1633 ist--;
1634 con=false;
1635#ifdef debug
1636 G4cout<<
"G4QFragmentation::Construct: *** IST Decremented ***"<<
G4endl;
1637#endif
1638 }
1639 else
1640 {
1641 con=true;
1642#ifdef debug
1643 G4cout<<
"G4QFragmentation::Construct: *** IST Begin ***"<<
G4endl;
1644#endif
1645 break;
1646 }
1647 }
1648 else
1649 {
1650#ifdef debug
1651 G4cout<<
"-Warning-G4QFragm::Const:S4M="<<cS4M<<
",M2="<<cSM2<<
" Leave ASIS"<<
G4endl;
1652#endif
1653 ++problem;
1654 con=false;
1655 }
1656 }
1657 else con=false;
1658 }
1659#ifdef debug
1660 G4cout<<
"G4QFragmentation::Construct: *** IST While *** , con="<<con<<
G4endl;
1661#endif
1662 }
1663#ifdef edebug
1664
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++)
1671 {
1673 t4M+=strI4M;
1674 G4int sChg=strings[i]->GetCharge();
1675 rC-=sChg;
1676 G4int sBaN=strings[i]->GetBaryonNumber();
1677 rB-=sBaN;
1678 G4cout<<
"-EMCLS-G4QFragm::Construct: St#"<<i<<
", 4M="<<strI4M<<
", M="<<strI4M.
m()
1679 <<
", C="<<sChg<<
", B="<<sBaN<<
G4endl;
1680 }
1681 G4cout<<
"-EMCLS-G4QFragm::Construct:r4M="<<t4M-totLS4M<<
",rC="<<rC<<
",rB="<<rB<<
G4endl;
1682#endif
1683
1684
1685
1686#ifdef debug
1687 G4cout<<
"G4QFragmentation::Construct: problem="<<problem<<
G4endl;
1688#endif
1689 if(problem)
1690 {
1691 G4int nOfStr=strings.size();
1692#ifdef debug
1693 G4cout<<
"G4QFragmentation::Construct:SecurityDiQaDiQReduction,#OfStr="<<nOfStr<<
G4endl;
1694#endif
1695 for (
G4int astring=0; astring < nOfStr; astring++)
1696 {
1705 {
1706#ifdef debug
1707 G4cout<<
"G4QFragmentation::Constr:TrySelfReduString,L="<<sPDG<<
",R="<<nPDG<<
G4endl;
1708#endif
1710 {
1713#ifdef debug
1714 G4cout<<
"+G4QFragm::Const:#"<<astring<<
" Reduced, L="<<sPDG<<
",R="<<nPDG<<
G4endl;
1715#endif
1716 }
1717#ifdef debug
1718 else G4cout<<
"--*--G4QFragm::Const:#"<<astring<<
" DQ-aDQ reduction Failed"<<
G4endl;
1719#endif
1720 }
1721 else if(sPDG==3 && nPDG==-3)
1722 {
1723 sPDG= 1;
1724 nPDG=-1;
1727 }
1728 else if(sPDG==-3 && nPDG==3)
1729 {
1730 sPDG=-1;
1731 nPDG= 1;
1734 }
1735 }
1737 }
1738#ifdef edebug
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++)
1745 {
1747 u4M+=strI4M;
1748 G4int sChg=strings[i]->GetCharge();
1749 rCh-=sChg;
1750 G4int sBaN=strings[i]->GetBaryonNumber();
1751 rBa-=sBaN;
1752 G4cout<<
"-EMCLS-G4QFragm::Construct: St#"<<i<<
", 4M="<<strI4M<<
", M="<<strI4M.
m()
1753 <<
", C="<<sChg<<
", B="<<sBaN<<
G4endl;
1754 }
1755 G4cout<<
"-EMCLS-G4QFragm::Construct:r4M="<<u4M-totLS4M<<
",rC="<<rCh<<
",rB="<<rBa<<
G4endl;
1756#endif
1757}
std::vector< G4QHadron * > G4QHadronVector
std::vector< G4QInteraction * > G4QInteractionVector
std::vector< G4QPartonPair * > G4QPartonPairVector
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
HepLorentzVector & boost(double, double, double)
static G4QCHIPSWorld * Get()
G4QHadronVector * ProjFragment(G4int pPDG, G4LorentzVector p4M, G4int tgZ, G4int tgN)
static G4QDiffractionRatio * GetPointer()
G4double GetRatio(G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
G4bool IsSingleDiffractive()
G4bool ExciteDiffParticipants(G4QHadron *aPartner, G4QHadron *bPartner) const
G4bool ExciteSingDiffParticipants(G4QHadron *aPartner, G4QHadron *bPartner) const
G4LorentzVector Get4Momentum() const
G4int GetBaryonNumber() const
void IncrementCollisionCount(G4int aCount)
G4QParton * GetNextParton()
const G4ThreeVector & GetPosition() const
std::list< G4QParton * > GetColor()
std::list< G4QParton * > GetAntiColor()
G4QParton * GetNextAntiParton()
void Set4Momentum(const G4LorentzVector &aMom)
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)
G4double GetThickness(G4double b)
G4LorentzVector GetNucleons4Momentum()
G4QHadron * GetNextNucleon()
G4double CoulombBarrier(const G4double &cZ=1, const G4double &cA=1, G4double dZ=0., G4double dA=0.)
G4double GetGSMass() const
std::pair< G4int, G4int > MakeTwoBaryons(G4int L1, G4int L2, G4int R1, G4int R2)
void SetPDGCode(G4int aPDG)
const G4int & GetType() const
const G4LorentzVector & Get4Momentum() const
void Set4Momentum(const G4LorentzVector &aMomentum)
G4bool ReduceDiQADiQ(G4QParton *d1, G4QParton *d2)
G4QParton * GetLeftParton() 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 G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)