84 for (
G4int i = 0; i < 250; i++ ) {
140 if ( aNucleon )
delete aNucleon;
148 if ( aNucleon )
delete aNucleon;
168 <<
"FTF init Target A Z " << aNucleus.
GetA_asInt()
308 G4cout <<
"FTF PutOnMassShell Success? " << Success <<
G4endl;
320 G4cout <<
"FTF ExciteParticipants Success? " << Success <<
G4endl;
326 G4cout <<
"FTF BuildStrings ";
332 G4cout <<
"FTF BuildStrings " << theStrings <<
" OK" << G4endl
333 <<
"FTF GetResiduals of Nuclei " <<
G4endl;
346 std::vector< G4VSplitableHadron* > primaries;
351 if ( primaries.end() ==
352 std::find( primaries.begin(), primaries.end(), interaction.
GetProjectile() ) ) {
366 if ( aNucleon )
delete aNucleon;
368 NumberOfInvolvedNucleonsOfProjectile = 0;
373 if ( aNucleon )
delete aNucleon;
375 NumberOfInvolvedNucleonsOfTarget = 0;
378 G4cout <<
"End of FTF. Go to fragmentation" << G4endl
379 <<
"To continue - enter 1, to stop - ^C" <<
G4endl;
408 G4cout <<
"G4FTFModel::StoreInvolvedNucleon -------------" <<
G4endl;
423 while ( ( aProjectileNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
444 #ifdef debugReggeonCascade
445 G4cout <<
"G4FTFModel::ReggeonCascade -----------" <<
G4endl
454 for (
G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
481 Neighbour->
Hit( targetSplitable );
489 #ifdef debugReggeonCascade
490 G4cout <<
"Final NumberOfInvolvedNucleonsOfTarget "
500 for (
G4int InvPN = 0; InvPN < InitNINp; InvPN++ ) {
512 while ( ( Neighbour = theProjectileNucleus->
GetNextNucleon() ) ) {
527 Neighbour->
Hit( projectileSplitable );
535 #ifdef debugReggeonCascade
536 G4cout <<
"NumberOfInvolvedNucleonsOfProjectile "
546 G4bool isProjectileNucleus =
false;
548 isProjectileNucleus =
true;
551 #ifdef debugPutOnMassShell
553 if ( isProjectileNucleus ) {
554 G4cout <<
"PutOnMassShell for Nucleus_Nucleus " <<
G4endl;
559 if ( Pprojectile.z() < 0.0 ) {
571 #ifdef debugPutOnMassShell
577 if ( ! isOk )
return false;
586 if ( ! isProjectileNucleus ) {
587 Mprojectile = Pprojectile.mag();
588 M2projectile = Pprojectile.mag2();
589 SumMasses += Mprojectile + 20.0*
MeV;
591 #ifdef debugPutOnMassShell
592 G4cout <<
"Projectile : ";
597 if ( ! isOk )
return false;
604 #ifdef debugPutOnMassShell
605 G4cout <<
"Psum " << Psum/
GeV <<
" GeV" << G4endl <<
"SqrtS " << SqrtS/
GeV <<
" GeV" << G4endl
606 <<
"SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/
GeV <<
" "
607 << PrResidualMass/
GeV <<
" " << TargetResidualMass/
GeV <<
" GeV" <<
G4endl;
610 if ( SqrtS < SumMasses ) {
616 G4double savedSumMasses = SumMasses;
617 if ( isProjectileNucleus ) {
618 SumMasses -= std::sqrt(
sqr( PrResidualMass ) + PprojResidual.
perp2() );
620 + PprojResidual.
perp2() );
622 SumMasses -= std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.
perp2() );
624 + PtargetResidual.
perp2() );
626 if ( SqrtS < SumMasses ) {
627 SumMasses = savedSumMasses;
628 if ( isProjectileNucleus ) {
635 if ( isProjectileNucleus ) {
639 #ifdef debugPutOnMassShell
640 if ( isProjectileNucleus ) {
641 G4cout <<
"PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/
GeV <<
" "
644 G4cout <<
"TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/
GeV <<
" "
646 <<
"Sum masses " << SumMasses/
GeV <<
G4endl;
650 if ( isProjectileNucleus && thePrNucleus->
GetMassNumber() != 1 ) {
659 if ( ! isOk )
return false;
670 if ( Ptmp.pz() <= 0.0 ) {
677 if ( isProjectileNucleus ) {
679 YprojectileNucleus = Ptmp.
rapidity();
681 Ptmp = toCms*Ptarget;
682 G4double YtargetNucleus = Ptmp.rapidity();
686 if ( isProjectileNucleus ) {
693 #ifdef debugPutOnMassShell
694 if ( isProjectileNucleus ) {
695 G4cout <<
"Y projectileNucleus " << YprojectileNucleus <<
G4endl;
697 G4cout <<
"Y targetNucleus " << YtargetNucleus << G4endl
699 <<
" DcorP DcorT " << DcorP <<
" " << DcorT <<
" AveragePt2 " << AveragePt2 <<
G4endl;
706 G4int NumberOfTries = 0;
708 G4bool OuterSuccess =
true;
710 const G4int maxNumberOfLoops = 1000;
711 G4int loopCounter = 0;
714 const G4int maxNumberOfInnerLoops = 10000;
717 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
723 DcorP *= ScaleFactor;
724 DcorT *= ScaleFactor;
725 AveragePt2 *= ScaleFactor;
727 if ( isProjectileNucleus ) {
730 thePrNucleus, PprojResidual,
738 theTargetNucleus, PtargetResidual,
743 #ifdef debugPutOnMassShell
744 G4cout <<
"SqrtS, Mp+Mt, Mp, Mt " << SqrtS/
GeV <<
" "
745 << ( std::sqrt( M2proj ) + std::sqrt( M2target) )/
GeV <<
" "
746 << std::sqrt( M2proj )/
GeV <<
" " << std::sqrt( M2target )/
GeV <<
G4endl;
749 if ( ! isOk )
return false;
750 }
while ( ( SqrtS < std::sqrt( M2proj ) + std::sqrt( M2target ) ) &&
751 NumberOfTries < maxNumberOfInnerLoops );
752 if ( NumberOfTries >= maxNumberOfInnerLoops ) {
753 #ifdef debugPutOnMassShell
754 G4cout <<
"BAD situation: forced exit of the inner while loop!" <<
G4endl;
758 if ( isProjectileNucleus ) {
759 isOk =
CheckKinematics( S, SqrtS, M2proj, M2target, YprojectileNucleus,
true,
762 WminusTarget, WplusProjectile, OuterSuccess );
767 WminusTarget, WplusProjectile, OuterSuccess );
768 if ( ! isOk )
return false;
769 }
while ( ( ! OuterSuccess ) &&
770 ++loopCounter < maxNumberOfLoops );
771 if ( loopCounter >= maxNumberOfLoops ) {
772 #ifdef debugPutOnMassShell
773 G4cout <<
"BAD situation: forced exit of the while loop!" <<
G4endl;
784 if ( ! isProjectileNucleus ) {
786 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
787 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
788 Pprojectile.setPz( Pzprojectile );
789 Pprojectile.setE( Eprojectile );
791 #ifdef debugPutOnMassShell
792 G4cout <<
"Proj after in CMS " << Pprojectile <<
G4endl;
795 Pprojectile.transform( toLab );
804 #ifdef debugPutOnMassShell
814 #ifdef debugPutOnMassShell
818 if ( ! isOk )
return false;
822 #ifdef debugPutOnMassShell
832 #ifdef debugPutOnMassShell
836 if ( ! isOk )
return false;
840 #ifdef debugPutOnMassShell
853 #ifdef debugBuildString
854 G4cout <<
"G4FTFModel::ExciteParticipants() " <<
G4endl;
857 G4bool Successfull(
true );
859 if ( MaxNumOfInelCollisions > 0 ) {
861 if (
G4UniformRand() < ProbMaxNumber ) MaxNumOfInelCollisions++;
864 MaxNumOfInelCollisions = 1;
867 #ifdef debugBuildString
868 G4cout <<
"MaxNumOfInelCollisions MaxNumOfInelCollisions " << MaxNumOfInelCollisions <<
G4endl;
871 G4int CurrentInteraction( 0 );
876 CurrentInteraction++;
883 #ifdef debugBuildString
884 G4cout << G4endl <<
"Interaction # Status " << CurrentInteraction <<
" "
885 << collision.
GetStatus() << G4endl <<
"Pr* Tr* " << projectile <<
" "
886 << target << G4endl <<
"projectile->GetStatus target->GetStatus "
896 #ifdef debugBuildString
901 G4bool Annihilation =
false;
903 TargetNucleon, Annihilation );
904 if ( ! Result )
continue;
911 #ifdef debugBuildString
912 G4cout <<
"Inelastic interaction" << G4endl
913 <<
"MaxNumOfInelCollisions " << MaxNumOfInelCollisions <<
G4endl;
917 G4bool Annihilation =
false;
919 TargetNucleon, Annihilation );
920 if ( ! Result )
continue;
933 #ifdef debugBuildString
946 #ifdef debugBuildString
947 G4cout <<
"FTF excitation Non Successfull -> Elastic scattering "
953 #ifdef debugBuildString
954 G4cout <<
"Elastic scat. at rejection inelastic scattering" <<
G4endl;
968 #ifdef debugBuildString
977 if ( projectile == NextProjectileNucleon || target == NextTargetNucleon ) {
988 G4bool Annihilation =
true;
990 TargetNucleon, Annihilation );
991 if ( ! Result )
continue;
995 Successfull = Successfull ||
true;
997 #ifdef debugBuildString
998 G4cout <<
"Annihilation successfull. " <<
"*AdditionalString "
999 << AdditionalString <<
G4endl;
1024 #ifdef debugBuildString
1025 G4cout <<
"----------------------------- Final properties " << G4endl
1026 <<
"projectile->GetStatus target->GetStatus " << projectile->
GetStatus()
1027 <<
" " << target->
GetStatus() << G4endl <<
"projectile->GetSoftC target->GetSoftC "
1029 << G4endl <<
"ExciteParticipants() Successfull? " << Successfull <<
G4endl;
1047 G4cout <<
"AdjustNucleons ---------------------------------------" <<
G4endl
1051 <<
"Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
1054 <<
"Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
1066 G4int interactionCase = 0;
1076 interactionCase = 1;
1078 G4cout <<
"case 1, hA prcol=0 trcol=0, AA prcol#0 trcol=0" <<
G4endl;
1097 interactionCase = 2;
1117 interactionCase = 3;
1127 <<
"ProjectileResidualMassNumber ProjectileResidualCharge "
1136 ProjectileNucleon, SelectedTargetNucleon,
1137 TargetNucleon, Annihilation, common );
1138 G4bool returnResult =
false;
1139 if ( returnCode == 0 ) {
1140 returnResult =
true;
1141 }
else if ( returnCode == 1 ) {
1144 if ( returnResult ) {
1146 SelectedTargetNucleon, common );
1150 return returnResult;
1167 G4int returnCode = 99;
1172 if ( interactionCase == 1 ) {
1178 }
else if ( interactionCase == 2 ) {
1181 }
else if ( interactionCase == 3 ) {
1198 if ( interactionCase == 1 ) {
1219 }
else if ( interactionCase == 2 ) {
1239 G4cout <<
"SelectedTN.mag() PNMass + PResidualMass "
1243 }
else if ( interactionCase == 3 ) {
1278 G4cout <<
"PNucleonMass PResidualMass TNucleonMass TResidualMass " << common.
PNucleonMass
1286 if ( ! Annihilation ) {
1293 if ( interactionCase == 1 || interactionCase == 2 ) {
1305 }
else if ( interactionCase == 3 ) {
1307 G4cout <<
"SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy "
1328 if ( Annihilation ) {
1330 G4cout <<
"SqrtS < SumMasses - TNucleonMass " << common.
SqrtS <<
" "
1340 if ( interactionCase == 2 || interactionCase == 3 ) {
1354 if ( interactionCase == 1 || interactionCase == 2 ) {
1359 }
else if ( interactionCase == 3 ) {
1385 if ( interactionCase == 1 ) {
1387 }
else if ( interactionCase == 2 ) {
1389 }
else if ( interactionCase == 3 ) {
1399 if ( interactionCase == 1 || interactionCase == 3 ) {
1401 }
else if ( interactionCase == 2 ) {
1411 if ( interactionCase == 1 || interactionCase == 3 ) {
1424 if ( interactionCase == 2 || interactionCase == 3 ) {
1426 if ( interactionCase == 2 ) {
1443 return returnCode = 0;
1447 if ( interactionCase == 1 ) {
1453 }
else if ( interactionCase == 2 ) {
1459 }
else if ( interactionCase == 3 ) {
1471 return returnCode = 1;
1490 G4bool OuterSuccess =
true;
1491 const G4int maxNumberOfLoops = 1000;
1492 const G4int maxNumberOfTries = 10000;
1493 G4int loopCounter = 0;
1494 G4int NumberOfTries = 0;
1496 OuterSuccess =
true;
1497 G4bool loopCondition =
false;
1499 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1502 DcorP *= ScaleFactor;
1503 DcorT *= ScaleFactor;
1504 AveragePt2 *= ScaleFactor;
1511 if ( interactionCase == 1 ) {
1512 }
else if ( interactionCase == 2 ) {
1530 OuterSuccess =
false;
1531 loopCondition =
true;
1534 }
else if ( interactionCase == 3 ) {
1554 OuterSuccess =
false;
1555 loopCondition =
true;
1560 G4int numberOfTimesExecuteInnerLoop = 1;
1561 if ( interactionCase == 3 ) numberOfTimesExecuteInnerLoop = 2;
1562 for (
G4int iExecute = 0; iExecute < numberOfTimesExecuteInnerLoop; iExecute++ ) {
1564 G4bool InnerSuccess =
true;
1565 G4bool isTargetToBeHandled = ( interactionCase == 1 ||
1566 ( interactionCase == 3 && iExecute == 1 ) );
1568 if ( isTargetToBeHandled ) {
1574 const G4int maxNumberOfInnerLoops = 1000;
1575 G4int innerLoopCounter = 0;
1577 InnerSuccess =
true;
1578 if ( isTargetToBeHandled ) {
1580 if ( interactionCase == 1 ) {
1586 InnerSuccess =
false;
1598 InnerSuccess =
false;
1605 if ( interactionCase == 2 ) {
1614 InnerSuccess =
false;
1619 }
while ( ( ! InnerSuccess ) &&
1620 ++innerLoopCounter < maxNumberOfInnerLoops );
1621 if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
1623 G4cout <<
"BAD situation: forced exit of the inner while loop!" <<
G4endl;
1628 if ( isTargetToBeHandled ) {
1639 if ( interactionCase == 1 ) {
1645 }
else if ( interactionCase == 2 ) {
1649 <<
"TResidualMass PtResidual XplusResidual " << common.
TResidualMass <<
" "
1657 G4cout <<
"SqrtS < Mtarget + std::sqrt(M2projectile) " << common.
SqrtS <<
" "
1662 }
else if ( interactionCase == 3 ) {
1665 <<
"XplusNucleon XplusResidual " << common.
XplusNucleon
1680 + std::sqrt( common.
M2target ) ) );
1683 }
while ( loopCondition &&
1684 ++NumberOfTries < maxNumberOfTries );
1685 if ( NumberOfTries >= maxNumberOfTries ) {
1687 G4cout <<
"BAD situation: forced exit of the intermediate while loop!" <<
G4endl;
1693 G4double Yprojectile = 0.0, YprojectileNucleon = 0.0, Ytarget = 0.0, YtargetNucleon = 0.0;
1697 if ( interactionCase == 1 ) {
1699 + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.
SqrtS;
1708 G4cout <<
"DecayMomentum2 " << DecayMomentum2 <<
G4endl
1709 <<
"WminusTarget WplusProjectile " << common.
WminusTarget
1711 <<
"Yprojectile " << Yprojectile <<
G4endl;
1725 <<
"YtN Ypr YtN-Ypr " <<
" " << YtargetNucleon <<
" " << Yprojectile
1726 <<
" " << YtargetNucleon - Yprojectile <<
G4endl;
1729 Yprojectile < YtargetNucleon ) {
1730 OuterSuccess =
false;
1733 }
else if ( interactionCase == 2 ) {
1735 + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.
SqrtS;
1742 G4cout <<
"DecayMomentum2 " << DecayMomentum2 <<
G4endl
1743 <<
"WminusTarget WplusProjectile " << common.
WminusTarget
1745 <<
"Ytarget " << Ytarget <<
G4endl;
1759 <<
"YpN Ytr YpN-Ytr " <<
" " << YprojectileNucleon <<
" " << Ytarget
1760 <<
" " << YprojectileNucleon - Ytarget <<
G4endl;
1763 Ytarget > YprojectileNucleon ) {
1764 OuterSuccess =
false;
1767 }
else if ( interactionCase == 3 ) {
1769 + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.
SqrtS;
1791 YprojectileNucleon < YtargetNucleon ) {
1792 OuterSuccess =
false;
1797 }
while ( ( ! OuterSuccess ) &&
1798 ++loopCounter < maxNumberOfLoops );
1799 if ( loopCounter >= maxNumberOfLoops ) {
1801 G4cout <<
"BAD situation: forced exit of the while loop!" <<
G4endl;
1819 if ( interactionCase == 1 ) {
1822 }
else if ( interactionCase == 2 ) {
1827 }
else if ( interactionCase == 3 ) {
1843 if ( interactionCase == 1 ) {
1848 }
else if ( interactionCase == 2 ) {
1851 }
else if ( interactionCase == 3 ) {
1867 if ( interactionCase == 1 || interactionCase == 3 ) {
1872 G4cout <<
"TargetResidualMassNumber TargetResidualCharge TargetResidualExcitationEnergy "
1878 if ( interactionCase == 1 ) {
1903 if ( interactionCase == 2 || interactionCase == 3 ) {
1904 if ( interactionCase == 2 ) {
1914 G4cout <<
"ProjectileResidualMassNumber ProjectileResidualCharge ProjectileResidualExcitationEnergy "
1920 if ( interactionCase == 2 ) {
1960 std::vector< G4VSplitableHadron* > primaries;
1966 if ( primaries.end() == std::find( primaries.begin(), primaries.end(),
1973 #ifdef debugBuildString
1975 <<
"Number of projectile strings " << primaries.size() <<
G4endl;
1978 for (
unsigned int ahadron = 0; ahadron < primaries.size(); ahadron++ ) {
1979 G4bool isProjectile(
true );
1982 FirstString = 0; SecondString = 0;
1983 if ( primaries[ahadron]->GetStatus() == 0 ) {
1986 }
else if ( primaries[ahadron]->GetStatus() == 1
1987 && primaries[ahadron]->GetSoftCollisionCount() != 0 ) {
1990 }
else if ( primaries[ahadron]->GetStatus() == 1
1991 && primaries[ahadron]->GetSoftCollisionCount() == 0 ) {
1994 primaries[ahadron]->GetTimeOfCreation(),
1995 primaries[ahadron]->GetPosition(),
1998 }
else if(primaries[ahadron]->GetStatus() == 2) {
2001 primaries[ahadron]->GetTimeOfCreation(),
2002 primaries[ahadron]->GetPosition(),
2006 G4cout <<
"Something wrong in FTF Model Build String" <<
G4endl;
2009 if ( FirstString != 0 ) strings->push_back( FirstString );
2010 if ( SecondString != 0 ) strings->push_back( SecondString );
2012 #ifdef debugBuildString
2013 G4cout <<
"FirstString & SecondString? " << FirstString <<
" " << SecondString <<
G4endl;
2024 #ifdef debugBuildString
2026 G4cout <<
"Check 1 string " << strings->operator[](0)->GetRightParton()->GetPDGcode()
2027 <<
" " << strings->operator[](0)->GetLeftParton()->GetPDGcode() <<
G4endl <<
G4endl;
2036 #ifdef debugBuildString
2037 G4cout <<
"Building of projectile-like strings" <<
G4endl;
2040 G4bool isProjectile =
true;
2043 #ifdef debugBuildString
2044 G4cout <<
"Nucleon #, status, intCount " << ahadron <<
" "
2053 #ifdef debugBuildString
2054 G4cout << G4endl <<
"ahadron aProjectile Status " << ahadron <<
" " << aProjectile
2058 FirstString = 0; SecondString = 0;
2061 #ifdef debugBuildString
2062 G4cout <<
"Case1 aProjectile->GetStatus() == 0 " <<
G4endl;
2071 #ifdef debugBuildString
2072 G4cout <<
"Case2 aProjectile->GetStatus() !=0 St==1 SoftCol!=0" <<
G4endl;
2084 #ifdef debugBuildString
2085 G4cout <<
"Case3 aProjectile->GetStatus() !=0 St==1 SoftCol==0" <<
G4endl;
2095 #ifdef debugBuildString
2096 G4cout <<
" Strings are built for nucleon marked for an interaction, but"
2097 <<
" the interaction was skipped." <<
G4endl;
2103 #ifdef debugBuildString
2104 G4cout <<
"Case4 aProjectile->GetStatus() !=0 St==2 " <<
G4endl;
2114 #ifdef debugBuildString
2115 G4cout <<
" Strings are build for involved nucleon." <<
G4endl;
2120 #ifdef debugBuildString
2127 #ifdef debugBuildString
2133 if ( FirstString != 0 ) strings->push_back( FirstString );
2134 if ( SecondString != 0 ) strings->push_back( SecondString );
2138 #ifdef debugBuildString
2142 G4bool isProjectile =
false;
2146 #ifdef debugBuildString
2147 G4cout <<
"Nucleon #, status, intCount " << aNucleon <<
" " << ahadron <<
" "
2151 FirstString = 0 ; SecondString = 0;
2157 #ifdef debugBuildString
2166 #ifdef debugBuildString
2167 G4cout <<
" 2 case A string is build, nucleon was excited." <<
G4endl;
2184 #ifdef debugBuildString
2196 #ifdef debugBuildString
2200 }
else if ( aNucleon->
GetStatus() == 2 ||
2209 #ifdef debugBuildString
2215 #ifdef debugBuildString
2221 if ( FirstString != 0 ) strings->push_back( FirstString );
2222 if ( SecondString != 0 ) strings->push_back( SecondString );
2226 #ifdef debugBuildString
2231 isProjectile =
true;
2235 FirstString = 0; SecondString = 0;
2238 if ( FirstString != 0 ) strings->push_back( FirstString );
2239 if ( SecondString != 0 ) strings->push_back( SecondString );
2258 #ifdef debugFTFmodel
2259 G4cout <<
"GetResiduals(): HighEnergyInter? GetProjectileNucleus()?"
2265 #ifdef debugFTFmodel
2277 #ifdef debugFTFmodel
2299 residualMomentum +=
tmp;
2321 const G4int maxNumberOfLoops = 1000;
2322 G4int loopCounter = 0;
2324 C = ( Chigh + Clow ) / 2.0;
2337 if ( SumMasses > Mass ) Chigh =
C;
2340 }
while ( Chigh - Clow > 0.01 &&
2341 ++loopCounter < maxNumberOfLoops );
2342 if ( loopCounter >= maxNumberOfLoops ) {
2343 #ifdef debugFTFmodel
2344 G4cout <<
"BAD situation: forced exit of the first while loop in G4FTFModel::GetResidual" << G4endl
2345 <<
"\t return immediately from the method!" <<
G4endl;
2365 #ifdef debugFTFmodel
2367 << G4endl <<
"ProjectileResidualExcitationEnergy ProjectileResidual4Momentum "
2379 #ifdef debugFTFmodel
2397 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2401 residualMomentum +=
tmp;
2412 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2423 const G4int maxNumberOfLoops = 1000;
2424 G4int loopCounter = 0;
2426 C = ( Chigh + Clow ) / 2.0;
2431 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2440 if ( SumMasses > Mass) Chigh =
C;
2443 }
while ( Chigh - Clow > 0.01 &&
2444 ++loopCounter < maxNumberOfLoops );
2445 if ( loopCounter >= maxNumberOfLoops ) {
2446 #ifdef debugFTFmodel
2447 G4cout <<
"BAD situation: forced exit of the second while loop in G4FTFModel::GetResidual" << G4endl
2448 <<
"\t return immediately from the method!" <<
G4endl;
2455 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2466 #ifdef debugFTFmodel
2472 #ifdef debugFTFmodel
2473 G4cout <<
"Low energy interaction: Target nucleus --------------" << G4endl
2474 <<
"Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
2479 G4int NumberOfTargetParticipant( 0 );
2489 if ( NumberOfTargetParticipant != 0 ) {
2502 delete targetSplitable;
2503 targetSplitable = 0;
2504 aNucleon->
Hit( targetSplitable );
2509 #ifdef debugFTFmodel
2510 G4cout <<
"NumberOfTargetParticipant " << NumberOfTargetParticipant << G4endl
2516 #ifdef debugFTFmodel
2517 G4cout <<
"Low energy interaction: Projectile nucleus --------------" << G4endl
2518 <<
"Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
2523 G4int NumberOfProjectileParticipant( 0 );
2528 NumberOfProjectileParticipant++;
2531 #ifdef debugFTFmodel
2535 DeltaExcitationE = 0.0;
2538 if ( NumberOfProjectileParticipant != 0 ) {
2540 G4double( NumberOfProjectileParticipant );
2542 G4double( NumberOfProjectileParticipant );
2554 delete projectileSplitable;
2555 projectileSplitable = 0;
2556 aNucleon->
Hit( projectileSplitable );
2561 #ifdef debugFTFmodel
2562 G4cout <<
"NumberOfProjectileParticipant " << NumberOfProjectileParticipant << G4endl
2568 #ifdef debugFTFmodel
2569 G4cout <<
"End GetResiduals -----------------" <<
G4endl;
2581 if(AveragePt2 > 0.0) {
2582 if(maxPtSquare/AveragePt2 < 1.0
e+9) {
2584 (
G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
2593 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
2603 G4double& residualExcitationEnergy,
2605 G4int& residualMassNumber,
2606 G4int& residualCharge ) {
2618 if ( ! nucleus )
return false;
2620 G4double ExcitationEnergyPerWoundedNucleon =
2643 sumMasses += 20.0*
MeV;
2647 residualExcitationEnergy += -ExcitationEnergyPerWoundedNucleon*
G4Log(
G4UniformRand() );
2649 residualMassNumber--;
2656 #ifdef debugPutOnMassShell
2657 G4cout <<
"ExcitationEnergyPerWoundedNucleon " << ExcitationEnergyPerWoundedNucleon <<
G4endl
2658 <<
"\t Residual Charge, MassNumber " << residualCharge <<
" " << residualMassNumber
2659 <<
G4endl <<
"\t Initial Momentum " << nucleusMomentum
2660 <<
G4endl <<
"\t Residual Momentum " << residualMomentum <<
G4endl;
2662 residualMomentum.
setPz( 0.0 );
2663 residualMomentum.
setE( 0.0 );
2664 if ( residualMassNumber == 0 ) {
2666 residualExcitationEnergy = 0.0;
2669 GetIonMass( residualCharge, residualMassNumber );
2670 if ( residualMassNumber == 1 ) {
2671 residualExcitationEnergy = 0.0;
2673 residualMass += residualExcitationEnergy;
2675 sumMasses += std::sqrt(
sqr( residualMass ) + residualMomentum.
perp2() );
2684 const G4int numberOfInvolvedNucleons,
2702 if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 )
return false;
2704 const G4double probDeltaIsobar = 0.05;
2706 G4int maxNumberOfDeltas =
G4int( (sqrtS - sumMasses)/(400.0*
MeV) );
2707 G4int numberOfDeltas = 0;
2709 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2712 if (
G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
2714 if ( ! involvedNucleons[i] )
continue;
2721 G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4;
2730 if ( sqrtS < sumMasses + massDelta - massNuc ) {
2734 sumMasses += ( massDelta - massNuc );
2753 const G4int residualMassNumber,
2754 const G4int numberOfInvolvedNucleons,
2768 if ( ! nucleus )
return false;
2770 if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
2779 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2780 G4Nucleon* aNucleon = involvedNucleons[i];
2781 if ( ! aNucleon )
continue;
2786 const G4int maxNumberOfLoops = 1000;
2787 G4int loopCounter = 0;
2795 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2796 G4Nucleon* aNucleon = involvedNucleons[i];
2797 if ( ! aNucleon )
continue;
2804 G4double deltaPx = ( ptSum.
x() - pResidual.
x() ) / numberOfInvolvedNucleons;
2805 G4double deltaPy = ( ptSum.
y() - pResidual.
y() ) / numberOfInvolvedNucleons;
2807 SumMasses = residualMass;
2808 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2809 G4Nucleon* aNucleon = involvedNucleons[i];
2810 if ( ! aNucleon )
continue;
2814 +
sqr( px ) +
sqr( py ) );
2823 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2824 G4Nucleon* aNucleon = involvedNucleons[i];
2825 if ( ! aNucleon )
continue;
2830 if ( x < 0.0 || x > 1.0 ) {
2842 if ( xSum < 0.0 || xSum > 1.0 ) success =
false;
2844 if ( ! success )
continue;
2849 if ( residualMassNumber == 0 ) {
2850 delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
2857 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2858 G4Nucleon* aNucleon = involvedNucleons[i];
2859 if ( ! aNucleon )
continue;
2862 if ( residualMassNumber == 0 ) {
2863 if ( x <= 0.0 || x > 1.0 ) {
2868 if ( x <= 0.0 || x > 1.0 || xSum <= 0.0 || xSum > 1.0 ) {
2887 if ( ! success )
continue;
2889 if ( success && residualMassNumber != 0 ) {
2890 mass2 += (
sqr( residualMass ) + pResidual.
perp2() ) / xSum;
2894 #ifdef debugPutOnMassShell
2898 }
while ( ( ! success ) &&
2899 ++loopCounter < maxNumberOfLoops );
2900 if ( loopCounter >= maxNumberOfLoops ) {
2916 const G4bool isProjectileNucleus,
2917 const G4int numberOfInvolvedNucleons,
2932 G4double decayMomentum2 =
sqr( sValue ) +
sqr( projectileMass2 ) +
sqr( targetMass2 )
2933 - 2.0*( sValue*( projectileMass2 + targetMass2 )
2934 + projectileMass2*targetMass2 );
2935 targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
2937 projectileWplus = sqrtS - targetMass2/targetWminus;
2938 G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
2939 G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
2940 G4double projectileY = 0.5 *
G4Log( (projectileE + projectilePz)/
2941 (projectileE - projectilePz) );
2942 G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
2943 G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
2944 G4double targetY = 0.5 *
G4Log( (targetE + targetPz)/(targetE - targetPz) );
2946 #ifdef debugPutOnMassShell
2947 G4cout <<
"decayMomentum2 " << decayMomentum2 <<
G4endl
2948 <<
"\t targetWminus projectileWplus " << targetWminus <<
" " << projectileWplus <<
G4endl
2949 <<
"\t projectileY targetY " << projectileY <<
" " << targetY <<
G4endl;
2952 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2953 G4Nucleon* aNucleon = involvedNucleons[i];
2954 if ( ! aNucleon )
continue;
2959 G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
2960 G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
2961 if ( isProjectileNucleus ) {
2962 pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*
x);
2963 e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*
x);
2967 #ifdef debugPutOnMassShell
2968 G4cout <<
"i nY pY nY-AY AY " << i <<
" " << nucleonY <<
" " << projectileY <<
G4endl;
2971 if ( std::abs( nucleonY - nucleusY ) > 2 ||
2972 ( isProjectileNucleus && targetY > nucleonY ) ||
2973 ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
2986 const G4bool isProjectileNucleus,
2989 const G4int residualMassNumber,
2990 const G4int numberOfInvolvedNucleons,
3008 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3009 G4Nucleon* aNucleon = involvedNucleons[i];
3010 if ( ! aNucleon )
continue;
3012 residual3Momentum -= tmp.
vect();
3016 G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w *
x );
3017 G4double e = w * x / 2.0 + mt2 / ( 2.0 * w *
x );
3019 if ( isProjectileNucleus ) pz *= -1.0;
3028 G4double residualMt2 =
sqr( residualMass ) +
sqr( residual3Momentum.
x() )
3029 +
sqr( residual3Momentum.
y() );
3031 #ifdef debugPutOnMassShell
3032 G4cout <<
"w residual3Momentum.z() " << w <<
" " << residual3Momentum.
z() <<
G4endl;
3037 if ( residualMassNumber != 0 ) {
3038 residualPz = -w * residual3Momentum.
z() / 2.0 +
3039 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
3040 residualE = w * residual3Momentum.
z() / 2.0 +
3041 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
3043 if ( isProjectileNucleus ) residualPz *= -1.0;
3046 residual4Momentum.
setPx( residual3Momentum.
x() );
3047 residual4Momentum.
setPy( residual3Momentum.
y() );
3048 residual4Momentum.
setPz( residualPz );
3049 residual4Momentum.
setE( residualE );
3058 desc <<
" FTF (Fritiof) Model \n"
3059 <<
"The FTF model is based on the well-known FRITIOF \n"
3060 <<
"model (B. Andersson et al., Nucl. Phys. B281, 289 \n"
3061 <<
"(1987)). Its first program implementation was given\n"
3062 <<
"by B. Nilsson-Almquist and E. Stenlund (Comp. Phys.\n"
3063 <<
"Comm. 43, 387 (1987)). The Fritiof model assumes \n"
3064 <<
"that all hadron-hadron interactions are binary \n"
3065 <<
"reactions, h_1+h_2->h_1'+h_2' where h_1' and h_2' \n"
3066 <<
"are excited states of the hadrons with continuous \n"
3067 <<
"mass spectra. The excited hadrons are considered as\n"
3068 <<
"QCD-strings, and the corresponding LUND-string \n"
3069 <<
"fragmentation model is applied for a simulation of \n"
3070 <<
"their decays. \n"
3071 <<
" The Fritiof model assumes that in the course of \n"
3072 <<
"a hadron-nucleus interaction a string originated \n"
3073 <<
"from the projectile can interact with various intra\n"
3074 <<
"nuclear nucleons and becomes into highly excited \n"
3075 <<
"states. The probability of multiple interactions is\n"
3076 <<
"calculated in the Glauber approximation. A cascading\n"
3077 <<
"of secondary particles was neglected as a rule. Due\n"
3078 <<
"to these, the original Fritiof model fails to des- \n"
3079 <<
"cribe a nuclear destruction and slow particle spectra.\n"
3080 <<
" In order to overcome the difficulties we enlarge\n"
3081 <<
"the model by the reggeon theory inspired model of \n"
3082 <<
"nuclear desctruction (Kh. Abdel-Waged and V.V. Uzhi-\n"
3083 <<
"nsky, Phys. Atom. Nucl. 60, 828 (1997); Yad. Fiz. 60, 925\n"
3084 <<
"(1997)). Momenta of the nucleons ejected from a nuc-\n"
3085 <<
"leus in the reggeon cascading are sampled according\n"
3086 <<
"to a Fermi motion algorithm presented in (EMU-01 \n"
3087 <<
"Collaboration (M.I. Adamovich et al.) Zeit. fur Phys.\n"
3088 <<
"A358, 337 (1997)). \n"
3089 <<
" New features were also added to the Fritiof model\n"
3090 <<
"implemented in Geant4: a simulation of elastic had-\n"
3091 <<
"ron-nucleon scatterings, a simulation of binary \n"
3092 <<
"reactions like NN>NN* in hadron-nucleon interactions,\n"
3093 <<
"a separate simulation of single diffractive and non-\n"
3094 <<
" diffractive events. These allowed to describe after\n"
3095 <<
"model parameter tuning a wide set of experimental \n"
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
void SetStatus(const G4int aStatus)
G4LorentzVector TResidual4Momentum
G4double GetMaxPt2ofNuclearDestruction()
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4DiffractiveExcitation * theExcitation
G4double TargetResidualExcitationEnergy
G4bool FinalizeKinematics(const G4double w, const G4bool isProjectileNucleus, const G4LorentzRotation &boostFromCmsToLab, const G4double residualMass, const G4int residualMassNumber, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4LorentzVector &residual4Momentum)
const G4ParticleDefinition * GetDefinition() const
G4Parton * GetRightParton(void) const
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
G4FTFAnnihilation * theAnnihilation
HepLorentzVector & transform(const HepRotation &)
virtual void DoLorentzContraction(const G4LorentzVector &theBoost)=0
void SetStatus(G4int aValue)
CLHEP::Hep3Vector G4ThreeVector
G4bool AdjustNucleonsAlgorithm_Sampling(G4int interactionCase, CommonVariables &common)
G4double GetCofNuclearDestruction()
void operator()(G4VSplitableHadron *aH)
G4int NumberOfInvolvedNucleonsOfTarget
G4double GetCofNuclearDestructionPr()
static G4ParticleTable * GetParticleTable()
static constexpr double MeV
void SetMomentum(G4LorentzVector &aMomentum)
G4int TResidualMassNumber
G4LorentzVector TargetResidual4Momentum
void SetBindingEnergy(G4double anEnergy)
static G4AntiNeutron * AntiNeutron()
std::vector< G4ExcitedString * > G4ExcitedStringVector
G4Nucleon * GetTargetNucleon() const
G4VSplitableHadron * GetTarget() const
G4ThreeVector PtResidualP
G4double GetProbabilityOfAnnihilation()
void SetMomentum(const G4double x, const G4double y, const G4double z)
void GetList(const G4ReactionProduct &thePrimary, G4FTFParameters *theParameters)
G4double Mt2targetNucleon
static G4AntiProton * AntiProton()
static constexpr double perCent
G4double GetPt2ofNuclearDestruction()
Hep3Vector findBoostToCM() const
const G4String & GetParticleName() const
G4double PResidualExcitationEnergy
G4bool ExciteParticipants()
G4IonTable * GetIonTable() const
G4double GetPDGCharge() const
std::vector< G4VSplitableHadron * > theAdditionalString
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
G4double YprojectileNucleus
void setVect(const Hep3Vector &)
static int FASTCALL common(PROLOG_STATE *state, int tok)
G4double PzprojectileNucleon
const G4ParticleDefinition * GetDefinition() const
void InitForInteraction(const G4ParticleDefinition *, G4int theA, G4int theZ, G4double s)
static G4Proton * Proton()
G4double condition(const G4ErrorSymMatrix &m)
G4V3DNucleus * theProjectileNucleus
G4double GetPDGMass() const
G4double GetProbabilityOfElasticScatt()
virtual void ModelDescription(std::ostream &) const
G4FTFModel(const G4String &modelName="FTF")
G4double G4Log(G4double x)
G4Nucleon * TheInvolvedNucleonsOfProjectile[250]
G4VSplitableHadron * GetSplitableHadron() const
G4V3DNucleus * GetProjectileNucleus() const
void SetParticleType(G4Proton *aProton)
virtual void SetProjectileNucleus(G4V3DNucleus *aNucleus)
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
void SetTotalEnergy(const G4double en)
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
G4VSplitableHadron * GetProjectile() const
G4InteractionContent & GetInteraction()
void Init(const G4Nucleus &aNucleus, const G4DynamicParticle &aProjectile)
G4Nucleon * TheInvolvedNucleonsOfTarget[250]
G4int ProjectileResidualMassNumber
const G4ThreeVector & GetPosition() const
CLHEP::HepLorentzRotation G4LorentzRotation
G4LorentzVector ProjectileResidual4Momentum
G4double GetTotalEnergy() const
G4int NumberOfInvolvedNucleonsOfProjectile
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
void SetTimeOfCreation(G4double aTime)
G4ReactionProduct theProjectile
G4double GetBindingEnergy() const
HepLorentzRotation & rotateZ(double delta)
const G4ThreeVector & GetPosition() const
virtual void Init(G4int theZ, G4int theA)
G4int GetBaryonNumber() const
void AdjustNucleonsAlgorithm_afterSampling(G4int interactionCase, G4VSplitableHadron *SelectedAntiBaryon, G4VSplitableHadron *SelectedTargetNucleon, CommonVariables &common)
G4double EprojectileNucleon
static constexpr double twopi
G4double GetTotalMomentum() const
HepLorentzRotation inverse() const
G4ElasticHNScattering * theElastic
G4Parton * GetLeftParton(void) const
G4Nucleon * GetProjectileNucleon() const
G4LorentzVector Pprojectile
virtual void DoLorentzBoost(const G4LorentzVector &theBoost)=0
G4ThreeVector GetMomentum() const
G4double GetMaxNumberOfCollisions()
G4int AdjustNucleonsAlgorithm_beforeSampling(G4int interactionCase, G4VSplitableHadron *SelectedAntiBaryon, G4Nucleon *ProjectileNucleon, G4VSplitableHadron *SelectedTargetNucleon, G4Nucleon *TargetNucleon, G4bool Annihilation, CommonVariables &common)
G4int ProjectileResidualCharge
void SetThisPointer(G4VPartonStringModel *aPointer)
G4ExcitedStringVector * GetStrings()
virtual G4int GetMassNumber()=0
virtual G4Nucleon * GetNextNucleon()=0
G4double Mt2projectileNucleon
G4int PResidualMassNumber
virtual const G4LorentzVector & Get4Momentum() const
G4int GetSoftCollisionCount()
virtual const G4ParticleDefinition * GetDefinition() const
virtual void InitProjectileNucleus(G4int theZ, G4int theA)
G4bool SamplingNucleonKinematics(G4double averagePt2, const G4double maxPt2, G4double dCor, G4V3DNucleus *nucleus, const G4LorentzVector &pResidual, const G4double residualMass, const G4int residualMassNumber, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4double &mass2)
G4int TargetResidualMassNumber
void Hit(G4VSplitableHadron *aHit)
G4int TargetResidualCharge
static G4Neutron * Neutron()
void SetDefinition(const G4ParticleDefinition *aDefinition)
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double AbsoluteLevel)
G4FTFParticipants theParticipants
virtual G4bool Annihilate(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4VSplitableHadron *&AdditionalString, G4FTFParameters *theParameters) const
G4int GetPDGEncoding() const
G4GLOB_DLL std::ostream G4cout
G4ThreeVector PtResidualT
Hep3Vector boostVector() const
const G4LorentzVector & Get4Momentum() const
CLHEP::HepLorentzVector G4LorentzVector
HepLorentzRotation & rotateY(double delta)
void StoreInvolvedNucleon()
G4V3DNucleus * GetTargetNucleus() const
G4double GetExcitationEnergyPerWoundedNucleon()
G4double ProjectileResidualExcitationEnergy
G4double GetR2ofNuclearDestruction()
G4bool AdjustNucleons(G4VSplitableHadron *SelectedAntiBaryon, G4Nucleon *ProjectileNucleon, G4VSplitableHadron *SelectedTargetNucleon, G4Nucleon *TargetNucleon, G4bool Annihilation)
G4bool CheckKinematics(const G4double sValue, const G4double sqrtS, const G4double projectileMass2, const G4double targetMass2, const G4double nucleusY, const G4bool isProjectileNucleus, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4double &targetWminus, G4double &projectileWplus, G4bool &success)
virtual G4bool StartLoop()=0
static constexpr double GeV
G4bool GenerateDeltaIsobar(const G4double sqrtS, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4double &sumMasses)
G4double GetTimeOfCreation()
void Set4Momentum(const G4LorentzVector &a4Momentum)
HepLorentzVector & boost(double, double, double)
G4ExcitedStringVector * BuildStrings()
G4double TResidualExcitationEnergy
G4LorentzVector PResidual4Momentum
G4bool ComputeNucleusProperties(G4V3DNucleus *nucleus, G4LorentzVector &nucleusMomentum, G4LorentzVector &residualMomentum, G4double &sumMasses, G4double &residualExcitationEnergy, G4double &residualMass, G4int &residualMassNumber, G4int &residualCharge)
G4FTFParameters * theParameters
G4double GetDofNuclearDestruction()