82 for (
G4int i = 0; i < 250; i++ ) {
138 if ( aNucleon )
delete aNucleon;
146 if ( aNucleon )
delete aNucleon;
166 <<
"FTF init Target A Z " << aNucleus.
GetA_asInt()
313 G4cout <<
"FTF PutOnMassShell Success? " << Success <<
G4endl;
325 G4cout <<
"FTF ExciteParticipants Success? " << Success <<
G4endl;
331 G4cout <<
"FTF BuildStrings ";
337 G4cout <<
"FTF BuildStrings " << theStrings <<
" OK" << G4endl
338 <<
"FTF GetResiduals of Nuclei " <<
G4endl;
351 std::vector< G4VSplitableHadron* > primaries;
356 if ( primaries.end() ==
357 std::find( primaries.begin(), primaries.end(), interaction.
GetProjectile() ) ) {
371 if ( aNucleon )
delete aNucleon;
373 NumberOfInvolvedNucleonsOfProjectile = 0;
378 if ( aNucleon )
delete aNucleon;
380 NumberOfInvolvedNucleonsOfTarget = 0;
383 G4cout <<
"End of FTF. Go to fragmentation" << G4endl
384 <<
"To continue - enter 1, to stop - ^C" <<
G4endl;
413 G4cout <<
"G4FTFModel::StoreInvolvedNucleon -------------" <<
G4endl;
428 while ( ( aProjectileNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
449 #ifdef debugReggeonCascade
450 G4cout <<
"G4FTFModel::ReggeonCascade -----------" <<
G4endl
459 for (
G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
486 Neighbour->
Hit( targetSplitable );
494 #ifdef debugReggeonCascade
495 G4cout <<
"Final NumberOfInvolvedNucleonsOfTarget "
505 for (
G4int InvPN = 0; InvPN < InitNINp; InvPN++ ) {
517 while ( ( Neighbour = theProjectileNucleus->
GetNextNucleon() ) ) {
532 Neighbour->
Hit( projectileSplitable );
540 #ifdef debugReggeonCascade
541 G4cout <<
"NumberOfInvolvedNucleonsOfProjectile "
551 G4bool isProjectileNucleus =
false;
553 isProjectileNucleus =
true;
556 #ifdef debugPutOnMassShell
558 if ( isProjectileNucleus ) {
559 G4cout <<
"PutOnMassShell for Nucleus_Nucleus " <<
G4endl;
564 if ( Pprojectile.z() < 0.0 ) {
576 #ifdef debugPutOnMassShell
582 if ( ! isOk )
return false;
591 if ( ! isProjectileNucleus ) {
592 Mprojectile = Pprojectile.mag();
593 M2projectile = Pprojectile.mag2();
594 SumMasses += Mprojectile + 20.0*
MeV;
596 #ifdef debugPutOnMassShell
597 G4cout <<
"Projectile : ";
602 if ( ! isOk )
return false;
609 #ifdef debugPutOnMassShell
610 G4cout <<
"Psum " << Psum/
GeV <<
" GeV" << G4endl <<
"SqrtS " << SqrtS/
GeV <<
" GeV" << G4endl
611 <<
"SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/
GeV <<
" "
612 << PrResidualMass/
GeV <<
" " << TargetResidualMass/
GeV <<
" GeV" <<
G4endl;
615 if ( SqrtS < SumMasses ) {
621 G4double savedSumMasses = SumMasses;
622 if ( isProjectileNucleus ) {
623 SumMasses -= std::sqrt(
sqr( PrResidualMass ) + PprojResidual.
perp2() );
625 + PprojResidual.
perp2() );
627 SumMasses -= std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.
perp2() );
629 + PtargetResidual.
perp2() );
631 if ( SqrtS < SumMasses ) {
632 SumMasses = savedSumMasses;
633 if ( isProjectileNucleus ) {
640 if ( isProjectileNucleus ) {
644 #ifdef debugPutOnMassShell
645 if ( isProjectileNucleus ) {
646 G4cout <<
"PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/
GeV <<
" "
649 G4cout <<
"TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/
GeV <<
" "
651 <<
"Sum masses " << SumMasses/
GeV <<
G4endl;
655 if ( isProjectileNucleus && thePrNucleus->
GetMassNumber() != 1 ) {
664 if ( ! isOk )
return false;
675 if ( Ptmp.pz() <= 0.0 ) {
682 if ( isProjectileNucleus ) {
684 YprojectileNucleus = Ptmp.
rapidity();
686 Ptmp = toCms*Ptarget;
687 G4double YtargetNucleus = Ptmp.rapidity();
691 if ( isProjectileNucleus ) {
698 #ifdef debugPutOnMassShell
699 if ( isProjectileNucleus ) {
700 G4cout <<
"Y projectileNucleus " << YprojectileNucleus <<
G4endl;
702 G4cout <<
"Y targetNucleus " << YtargetNucleus << G4endl
704 <<
" DcorP DcorT " << DcorP <<
" " << DcorT <<
" AveragePt2 " << AveragePt2 <<
G4endl;
711 G4int NumberOfTries = 0;
713 G4bool OuterSuccess =
true;
715 const G4int maxNumberOfLoops = 1000;
716 G4int loopCounter = 0;
719 const G4int maxNumberOfInnerLoops = 10000;
722 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
728 DcorP *= ScaleFactor;
729 DcorT *= ScaleFactor;
730 AveragePt2 *= ScaleFactor;
732 if ( isProjectileNucleus ) {
735 thePrNucleus, PprojResidual,
743 theTargetNucleus, PtargetResidual,
748 #ifdef debugPutOnMassShell
749 G4cout <<
"SqrtS, Mp+Mt, Mp, Mt " << SqrtS/
GeV <<
" "
750 << ( std::sqrt( M2proj ) + std::sqrt( M2target) )/
GeV <<
" "
751 << std::sqrt( M2proj )/
GeV <<
" " << std::sqrt( M2target )/
GeV <<
G4endl;
754 if ( ! isOk )
return false;
755 }
while ( ( SqrtS < std::sqrt( M2proj ) + std::sqrt( M2target ) ) &&
756 NumberOfTries < maxNumberOfInnerLoops );
757 if ( NumberOfTries >= maxNumberOfInnerLoops ) {
758 #ifdef debugPutOnMassShell
759 G4cout <<
"BAD situation: forced exit of the inner while loop!" <<
G4endl;
763 if ( isProjectileNucleus ) {
764 isOk =
CheckKinematics( S, SqrtS, M2proj, M2target, YprojectileNucleus,
true,
767 WminusTarget, WplusProjectile, OuterSuccess );
772 WminusTarget, WplusProjectile, OuterSuccess );
773 if ( ! isOk )
return false;
774 }
while ( ( ! OuterSuccess ) &&
775 ++loopCounter < maxNumberOfLoops );
776 if ( loopCounter >= maxNumberOfLoops ) {
777 #ifdef debugPutOnMassShell
778 G4cout <<
"BAD situation: forced exit of the while loop!" <<
G4endl;
789 if ( ! isProjectileNucleus ) {
791 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
792 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
793 Pprojectile.setPz( Pzprojectile );
794 Pprojectile.setE( Eprojectile );
796 #ifdef debugPutOnMassShell
797 G4cout <<
"Proj after in CMS " << Pprojectile <<
G4endl;
800 Pprojectile.transform( toLab );
809 #ifdef debugPutOnMassShell
819 #ifdef debugPutOnMassShell
823 if ( ! isOk )
return false;
827 #ifdef debugPutOnMassShell
837 #ifdef debugPutOnMassShell
841 if ( ! isOk )
return false;
845 #ifdef debugPutOnMassShell
858 #ifdef debugBuildString
859 G4cout <<
"G4FTFModel::ExciteParticipants() " <<
G4endl;
864 if ( MaxNumOfInelCollisions > 0 ) {
866 if (
G4UniformRand() < ProbMaxNumber ) MaxNumOfInelCollisions++;
869 MaxNumOfInelCollisions = 1;
872 #ifdef debugBuildString
873 G4cout <<
"MaxNumOfInelCollisions per hadron/nucleon " << MaxNumOfInelCollisions <<
G4endl;
876 G4int CurrentInteraction( 0 );
879 G4bool InnerSuccess(
true );
881 CurrentInteraction++;
888 #ifdef debugBuildString
889 G4cout << G4endl <<
"Interaction # Status " << CurrentInteraction <<
" "
890 << collision.
GetStatus() << G4endl <<
"Pr* Tr* " << projectile <<
" "
891 << target << G4endl <<
"projectile->GetStatus target->GetStatus "
901 #ifdef debugBuildString
906 G4bool Annihilation =
false;
908 TargetNucleon, Annihilation );
909 if ( ! Result )
continue;
915 #ifdef debugBuildString
916 G4cout <<
"Inelastic interaction" << G4endl
917 <<
"MaxNumOfInelCollisions per hadron/nucleon " << MaxNumOfInelCollisions <<
G4endl;
921 G4bool Annihilation =
false;
923 TargetNucleon, Annihilation );
924 if ( ! Result )
continue;
937 #ifdef debugBuildString
946 #ifdef debugBuildString
947 G4cout <<
"FTF excitation Non InnerSuccess of Elastic scattering "
948 << InnerSuccess <<
G4endl;
952 #ifdef debugBuildString
953 G4cout <<
"Elastic scat. at rejection inelastic scattering" <<
G4endl;
965 #ifdef debugBuildString
974 if ( projectile == NextProjectileNucleon || target == NextTargetNucleon ) {
985 G4bool Annihilation =
true;
987 TargetNucleon, Annihilation );
988 if ( ! Result )
continue;
993 #ifdef debugBuildString
994 G4cout <<
"Annihilation successfull. " <<
"*AdditionalString "
995 << AdditionalString <<
G4endl;
1020 if( InnerSuccess ) Success =
true;
1022 #ifdef debugBuildString
1023 G4cout <<
"----------------------------- Final properties " << G4endl
1024 <<
"projectile->GetStatus target->GetStatus " << projectile->
GetStatus()
1025 <<
" " << target->
GetStatus() << G4endl <<
"projectile->GetSoftC target->GetSoftC "
1027 << G4endl <<
"ExciteParticipants() Success? " << Success <<
G4endl;
1045 G4cout <<
"AdjustNucleons ---------------------------------------" <<
G4endl
1049 <<
"Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
1052 <<
"Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
1064 G4int interactionCase = 0;
1074 interactionCase = 1;
1076 G4cout <<
"case 1, hA prcol=0 trcol=0, AA prcol#0 trcol=0" <<
G4endl;
1095 interactionCase = 2;
1115 interactionCase = 3;
1125 <<
"ProjectileResidualMassNumber ProjectileResidualCharge "
1134 ProjectileNucleon, SelectedTargetNucleon,
1135 TargetNucleon, Annihilation, common );
1136 G4bool returnResult =
false;
1137 if ( returnCode == 0 ) {
1138 returnResult =
true;
1139 }
else if ( returnCode == 1 ) {
1142 if ( returnResult ) {
1144 SelectedTargetNucleon, common );
1148 return returnResult;
1165 G4int returnCode = 99;
1170 if ( interactionCase == 1 ) {
1176 }
else if ( interactionCase == 2 ) {
1179 }
else if ( interactionCase == 3 ) {
1196 if ( interactionCase == 1 ) {
1217 }
else if ( interactionCase == 2 ) {
1237 G4cout <<
"SelectedTN.mag() PNMass + PResidualMass "
1241 }
else if ( interactionCase == 3 ) {
1276 G4cout <<
"PNucleonMass PResidualMass TNucleonMass TResidualMass " << common.
PNucleonMass
1284 if ( ! Annihilation ) {
1291 if ( interactionCase == 1 || interactionCase == 2 ) {
1303 }
else if ( interactionCase == 3 ) {
1305 G4cout <<
"SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy "
1326 if ( Annihilation ) {
1328 G4cout <<
"SqrtS < SumMasses - TNucleonMass " << common.
SqrtS <<
" "
1338 if ( interactionCase == 2 || interactionCase == 3 ) {
1352 if ( interactionCase == 1 || interactionCase == 2 ) {
1357 }
else if ( interactionCase == 3 ) {
1383 if ( interactionCase == 1 ) {
1385 }
else if ( interactionCase == 2 ) {
1387 }
else if ( interactionCase == 3 ) {
1402 if ( interactionCase == 1 || interactionCase == 3 ) {
1404 }
else if ( interactionCase == 2 ) {
1419 if ( interactionCase == 1 || interactionCase == 3 ) {
1428 common.
Ptmp.
setPx( -saveSelectedTargetNucleon4Momentum.
x() );
1429 common.
Ptmp.
setPy( -saveSelectedTargetNucleon4Momentum.
y() );
1430 common.
Ptmp.
setPz( -saveSelectedTargetNucleon4Momentum.
z() );
1440 if ( interactionCase == 2 || interactionCase == 3 ) {
1442 if ( interactionCase == 2 ) {
1455 common.
Ptmp.
setPx( -saveSelectedAntiBaryon4Momentum.
x() );
1456 common.
Ptmp.
setPy( -saveSelectedAntiBaryon4Momentum.
y() );
1457 common.
Ptmp.
setPz( -saveSelectedAntiBaryon4Momentum.
z() );
1467 return returnCode = 0;
1471 if ( interactionCase == 1 ) {
1477 }
else if ( interactionCase == 2 ) {
1483 }
else if ( interactionCase == 3 ) {
1495 return returnCode = 1;
1515 G4bool OuterSuccess =
true;
1516 const G4int maxNumberOfLoops = 1000;
1517 const G4int maxNumberOfTries = 10000;
1518 G4int loopCounter = 0;
1519 G4int NumberOfTries = 0;
1521 OuterSuccess =
true;
1522 G4bool loopCondition =
false;
1524 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1527 DcorP *= ScaleFactor;
1528 DcorT *= ScaleFactor;
1529 AveragePt2 *= ScaleFactor;
1536 if ( interactionCase == 1 ) {
1537 }
else if ( interactionCase == 2 ) {
1555 OuterSuccess =
false;
1556 loopCondition =
true;
1559 }
else if ( interactionCase == 3 ) {
1579 OuterSuccess =
false;
1580 loopCondition =
true;
1585 G4int numberOfTimesExecuteInnerLoop = 1;
1586 if ( interactionCase == 3 ) numberOfTimesExecuteInnerLoop = 2;
1587 for (
G4int iExecute = 0; iExecute < numberOfTimesExecuteInnerLoop; iExecute++ ) {
1589 G4bool InnerSuccess =
true;
1590 G4bool isTargetToBeHandled = ( interactionCase == 1 ||
1591 ( interactionCase == 3 && iExecute == 1 ) );
1593 if ( isTargetToBeHandled ) {
1599 const G4int maxNumberOfInnerLoops = 1000;
1600 G4int innerLoopCounter = 0;
1602 InnerSuccess =
true;
1603 if ( isTargetToBeHandled ) {
1605 if ( interactionCase == 1 ) {
1611 InnerSuccess =
false;
1623 InnerSuccess =
false;
1630 if ( interactionCase == 2 ) {
1639 InnerSuccess =
false;
1644 }
while ( ( ! InnerSuccess ) &&
1645 ++innerLoopCounter < maxNumberOfInnerLoops );
1646 if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
1648 G4cout <<
"BAD situation: forced exit of the inner while loop!" <<
G4endl;
1653 if ( isTargetToBeHandled ) {
1664 if ( interactionCase == 1 ) {
1670 }
else if ( interactionCase == 2 ) {
1674 <<
"TResidualMass PtResidual XplusResidual " << common.
TResidualMass <<
" "
1682 G4cout <<
"SqrtS < Mtarget + std::sqrt(M2projectile) " << common.
SqrtS <<
" "
1687 }
else if ( interactionCase == 3 ) {
1690 <<
"XplusNucleon XplusResidual " << common.
XplusNucleon
1705 + std::sqrt( common.
M2target ) ) );
1708 }
while ( loopCondition &&
1709 ++NumberOfTries < maxNumberOfTries );
1710 if ( NumberOfTries >= maxNumberOfTries ) {
1712 G4cout <<
"BAD situation: forced exit of the intermediate while loop!" <<
G4endl;
1718 G4double Yprojectile = 0.0, YprojectileNucleon = 0.0, Ytarget = 0.0, YtargetNucleon = 0.0;
1722 if ( interactionCase == 1 ) {
1724 + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.
SqrtS;
1733 G4cout <<
"DecayMomentum2 " << DecayMomentum2 <<
G4endl
1734 <<
"WminusTarget WplusProjectile " << common.
WminusTarget
1736 <<
"Yprojectile " << Yprojectile <<
G4endl;
1750 <<
"YtN Ypr YtN-Ypr " <<
" " << YtargetNucleon <<
" " << Yprojectile
1751 <<
" " << YtargetNucleon - Yprojectile <<
G4endl;
1754 Yprojectile < YtargetNucleon ) {
1755 OuterSuccess =
false;
1758 }
else if ( interactionCase == 2 ) {
1760 + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.
SqrtS;
1767 G4cout <<
"DecayMomentum2 " << DecayMomentum2 <<
G4endl
1768 <<
"WminusTarget WplusProjectile " << common.
WminusTarget
1770 <<
"Ytarget " << Ytarget <<
G4endl;
1784 <<
"YpN Ytr YpN-Ytr " <<
" " << YprojectileNucleon <<
" " << Ytarget
1785 <<
" " << YprojectileNucleon - Ytarget <<
G4endl;
1788 Ytarget > YprojectileNucleon ) {
1789 OuterSuccess =
false;
1792 }
else if ( interactionCase == 3 ) {
1794 + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.
SqrtS;
1816 YprojectileNucleon < YtargetNucleon ) {
1817 OuterSuccess =
false;
1822 }
while ( ( ! OuterSuccess ) &&
1823 ++loopCounter < maxNumberOfLoops );
1824 if ( loopCounter >= maxNumberOfLoops ) {
1826 G4cout <<
"BAD situation: forced exit of the while loop!" <<
G4endl;
1844 if ( interactionCase == 1 ) {
1847 }
else if ( interactionCase == 2 ) {
1852 }
else if ( interactionCase == 3 ) {
1868 if ( interactionCase == 1 ) {
1873 }
else if ( interactionCase == 2 ) {
1876 }
else if ( interactionCase == 3 ) {
1892 if ( interactionCase == 1 || interactionCase == 3 ) {
1897 G4cout <<
"TargetResidualMassNumber TargetResidualCharge TargetResidualExcitationEnergy "
1903 if ( interactionCase == 1 ) {
1928 if ( interactionCase == 2 || interactionCase == 3 ) {
1929 if ( interactionCase == 2 ) {
1939 G4cout <<
"ProjectileResidualMassNumber ProjectileResidualCharge ProjectileResidualExcitationEnergy "
1945 if ( interactionCase == 2 ) {
1984 std::vector< G4VSplitableHadron* > primaries;
1990 if ( primaries.end() == std::find( primaries.begin(), primaries.end(),
1997 #ifdef debugBuildString
1999 <<
"Number of projectile strings " << primaries.size() <<
G4endl;
2002 for (
unsigned int ahadron = 0; ahadron < primaries.size(); ahadron++ ) {
2003 G4bool isProjectile(
true );
2006 FirstString = 0; SecondString = 0;
2007 if ( primaries[ahadron]->GetStatus() == 0 ) {
2010 }
else if ( primaries[ahadron]->GetStatus() == 1
2011 && primaries[ahadron]->GetSoftCollisionCount() != 0 ) {
2014 }
else if ( primaries[ahadron]->GetStatus() == 1
2015 && primaries[ahadron]->GetSoftCollisionCount() == 0 ) {
2018 primaries[ahadron]->GetTimeOfCreation(),
2019 primaries[ahadron]->GetPosition(),
2022 }
else if (primaries[ahadron]->GetStatus() == 2) {
2025 primaries[ahadron]->GetTimeOfCreation(),
2026 primaries[ahadron]->GetPosition(),
2030 G4cout <<
"Something wrong in FTF Model Build String" <<
G4endl;
2033 if ( FirstString != 0 ) strings->push_back( FirstString );
2034 if ( SecondString != 0 ) strings->push_back( SecondString );
2036 #ifdef debugBuildString
2037 G4cout <<
"FirstString & SecondString? " << FirstString <<
" " << SecondString <<
G4endl;
2048 #ifdef debugBuildString
2050 G4cout <<
"Check 1 string " << strings->operator[](0)->GetRightParton()->GetPDGcode()
2051 <<
" " << strings->operator[](0)->GetLeftParton()->GetPDGcode() <<
G4endl <<
G4endl;
2060 #ifdef debugBuildString
2061 G4cout <<
"Building of projectile-like strings" <<
G4endl;
2064 G4bool isProjectile =
true;
2067 #ifdef debugBuildString
2068 G4cout <<
"Nucleon #, status, intCount " << ahadron <<
" "
2077 #ifdef debugBuildString
2078 G4cout << G4endl <<
"ahadron aProjectile Status " << ahadron <<
" " << aProjectile
2082 FirstString = 0; SecondString = 0;
2085 #ifdef debugBuildString
2086 G4cout <<
"Case1 aProjectile->GetStatus() == 0 " <<
G4endl;
2095 #ifdef debugBuildString
2096 G4cout <<
"Case2 aProjectile->GetStatus() !=0 St==1 SoftCol!=0" <<
G4endl;
2108 #ifdef debugBuildString
2109 G4cout <<
"Case3 aProjectile->GetStatus() !=0 St==1 SoftCol==0" <<
G4endl;
2119 #ifdef debugBuildString
2120 G4cout <<
" Strings are built for nucleon marked for an interaction, but"
2121 <<
" the interaction was skipped." <<
G4endl;
2127 #ifdef debugBuildString
2128 G4cout <<
"Case4 aProjectile->GetStatus() !=0 St==2 " <<
G4endl;
2138 #ifdef debugBuildString
2139 G4cout <<
" Strings are build for involved nucleon." <<
G4endl;
2144 #ifdef debugBuildString
2151 #ifdef debugBuildString
2157 if ( FirstString != 0 ) strings->push_back( FirstString );
2158 if ( SecondString != 0 ) strings->push_back( SecondString );
2162 #ifdef debugBuildString
2166 G4bool isProjectile =
false;
2170 #ifdef debugBuildString
2171 G4cout <<
"Nucleon #, status, intCount " << aNucleon <<
" " << ahadron <<
" "
2175 FirstString = 0 ; SecondString = 0;
2181 #ifdef debugBuildString
2190 #ifdef debugBuildString
2191 G4cout <<
" 2 case A string is build, nucleon was excited." <<
G4endl;
2208 #ifdef debugBuildString
2220 #ifdef debugBuildString
2224 }
else if ( aNucleon->
GetStatus() == 2 ||
2233 #ifdef debugBuildString
2239 #ifdef debugBuildString
2245 if ( FirstString != 0 ) strings->push_back( FirstString );
2246 if ( SecondString != 0 ) strings->push_back( SecondString );
2250 #ifdef debugBuildString
2255 isProjectile =
true;
2259 FirstString = 0; SecondString = 0;
2262 if ( FirstString != 0 ) strings->push_back( FirstString );
2263 if ( SecondString != 0 ) strings->push_back( SecondString );
2282 #ifdef debugFTFmodel
2283 G4cout <<
"GetResiduals(): HighEnergyInter? GetProjectileNucleus()?"
2289 #ifdef debugFTFmodel
2301 #ifdef debugFTFmodel
2323 residualMomentum +=
tmp;
2345 const G4int maxNumberOfLoops = 1000;
2346 G4int loopCounter = 0;
2348 C = ( Chigh + Clow ) / 2.0;
2361 if ( SumMasses > Mass ) Chigh =
C;
2364 }
while ( Chigh - Clow > 0.01 &&
2365 ++loopCounter < maxNumberOfLoops );
2366 if ( loopCounter >= maxNumberOfLoops ) {
2367 #ifdef debugFTFmodel
2368 G4cout <<
"BAD situation: forced exit of the first while loop in G4FTFModel::GetResidual" << G4endl
2369 <<
"\t return immediately from the method!" <<
G4endl;
2389 #ifdef debugFTFmodel
2391 << G4endl <<
"ProjectileResidualExcitationEnergy ProjectileResidual4Momentum "
2403 #ifdef debugFTFmodel
2421 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2425 residualMomentum +=
tmp;
2436 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2447 const G4int maxNumberOfLoops = 1000;
2448 G4int loopCounter = 0;
2450 C = ( Chigh + Clow ) / 2.0;
2455 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2464 if ( SumMasses > Mass) Chigh =
C;
2467 }
while ( Chigh - Clow > 0.01 &&
2468 ++loopCounter < maxNumberOfLoops );
2469 if ( loopCounter >= maxNumberOfLoops ) {
2470 #ifdef debugFTFmodel
2471 G4cout <<
"BAD situation: forced exit of the second while loop in G4FTFModel::GetResidual" << G4endl
2472 <<
"\t return immediately from the method!" <<
G4endl;
2479 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2490 #ifdef debugFTFmodel
2496 #ifdef debugFTFmodel
2497 G4cout <<
"Low energy interaction: Target nucleus --------------" << G4endl
2498 <<
"Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
2503 G4int NumberOfTargetParticipant( 0 );
2513 if ( NumberOfTargetParticipant != 0 ) {
2526 delete targetSplitable;
2527 targetSplitable = 0;
2528 aNucleon->
Hit( targetSplitable );
2533 #ifdef debugFTFmodel
2534 G4cout <<
"NumberOfTargetParticipant " << NumberOfTargetParticipant << G4endl
2540 #ifdef debugFTFmodel
2541 G4cout <<
"Low energy interaction: Projectile nucleus --------------" << G4endl
2542 <<
"Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
2547 G4int NumberOfProjectileParticipant( 0 );
2554 #ifdef debugFTFmodel
2558 DeltaExcitationE = 0.0;
2561 if ( NumberOfProjectileParticipant != 0 ) {
2575 delete projectileSplitable;
2576 projectileSplitable = 0;
2577 aNucleon->
Hit( projectileSplitable );
2582 #ifdef debugFTFmodel
2583 G4cout <<
"NumberOfProjectileParticipant " << NumberOfProjectileParticipant << G4endl
2589 #ifdef debugFTFmodel
2590 G4cout <<
"End GetResiduals -----------------" <<
G4endl;
2602 if (AveragePt2 > 0.0) {
2603 if (maxPtSquare/AveragePt2 < 1.0
e+9) {
2605 (
G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
2614 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
2624 G4double& residualExcitationEnergy,
2626 G4int& residualMassNumber,
2627 G4int& residualCharge ) {
2639 if ( ! nucleus )
return false;
2641 G4double ExcitationEnergyPerWoundedNucleon =
2664 sumMasses += 20.0*
MeV;
2667 residualExcitationEnergy += -ExcitationEnergyPerWoundedNucleon*
G4Log(
G4UniformRand() );
2669 residualMassNumber--;
2676 #ifdef debugPutOnMassShell
2677 G4cout <<
"ExcitationEnergyPerWoundedNucleon " << ExcitationEnergyPerWoundedNucleon <<
G4endl
2678 <<
"\t Residual Charge, MassNumber " << residualCharge <<
" " << residualMassNumber
2679 <<
G4endl <<
"\t Initial Momentum " << nucleusMomentum
2680 <<
G4endl <<
"\t Residual Momentum " << residualMomentum <<
G4endl;
2682 residualMomentum.
setPz( 0.0 );
2683 residualMomentum.
setE( 0.0 );
2684 if ( residualMassNumber == 0 ) {
2686 residualExcitationEnergy = 0.0;
2689 GetIonMass( residualCharge, residualMassNumber );
2690 if ( residualMassNumber == 1 ) {
2691 residualExcitationEnergy = 0.0;
2693 residualMass += residualExcitationEnergy;
2695 sumMasses += std::sqrt(
sqr( residualMass ) + residualMomentum.
perp2() );
2704 const G4int numberOfInvolvedNucleons,
2722 if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 )
return false;
2724 const G4double probDeltaIsobar = 0.05;
2726 G4int maxNumberOfDeltas =
G4int( (sqrtS - sumMasses)/(400.0*
MeV) );
2727 G4int numberOfDeltas = 0;
2729 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2731 if (
G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
2733 if ( ! involvedNucleons[i] )
continue;
2740 G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4;
2749 if ( sqrtS < sumMasses + massDelta - massNuc ) {
2753 sumMasses += ( massDelta - massNuc );
2771 const G4int residualMassNumber,
2772 const G4int numberOfInvolvedNucleons,
2786 if ( ! nucleus )
return false;
2788 if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
2797 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2798 G4Nucleon* aNucleon = involvedNucleons[i];
2799 if ( ! aNucleon )
continue;
2804 const G4int maxNumberOfLoops = 1000;
2805 G4int loopCounter = 0;
2813 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2814 G4Nucleon* aNucleon = involvedNucleons[i];
2815 if ( ! aNucleon )
continue;
2822 G4double deltaPx = ( ptSum.
x() - pResidual.
x() ) / numberOfInvolvedNucleons;
2823 G4double deltaPy = ( ptSum.
y() - pResidual.
y() ) / numberOfInvolvedNucleons;
2825 SumMasses = residualMass;
2826 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2827 G4Nucleon* aNucleon = involvedNucleons[i];
2828 if ( ! aNucleon )
continue;
2832 +
sqr( px ) +
sqr( py ) );
2841 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2842 G4Nucleon* aNucleon = involvedNucleons[i];
2843 if ( ! aNucleon )
continue;
2848 if ( x < 0.0 || x > 1.0 ) {
2860 if ( xSum < 0.0 || xSum > 1.0 ) success =
false;
2862 if ( ! success )
continue;
2867 if ( residualMassNumber == 0 ) {
2868 delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
2875 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2876 G4Nucleon* aNucleon = involvedNucleons[i];
2877 if ( ! aNucleon )
continue;
2880 if ( residualMassNumber == 0 ) {
2881 if ( x <= 0.0 || x > 1.0 ) {
2886 if ( x <= 0.0 || x > 1.0 || xSum <= 0.0 || xSum > 1.0 ) {
2905 if ( ! success )
continue;
2907 if ( success && residualMassNumber != 0 ) {
2908 mass2 += (
sqr( residualMass ) + pResidual.
perp2() ) / xSum;
2912 #ifdef debugPutOnMassShell
2916 }
while ( ( ! success ) &&
2917 ++loopCounter < maxNumberOfLoops );
2918 if ( loopCounter >= maxNumberOfLoops ) {
2934 const G4bool isProjectileNucleus,
2935 const G4int numberOfInvolvedNucleons,
2950 G4double decayMomentum2 =
sqr( sValue ) +
sqr( projectileMass2 ) +
sqr( targetMass2 )
2951 - 2.0*( sValue*( projectileMass2 + targetMass2 )
2952 + projectileMass2*targetMass2 );
2953 targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
2955 projectileWplus = sqrtS - targetMass2/targetWminus;
2956 G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
2957 G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
2958 G4double projectileY = 0.5 *
G4Log( (projectileE + projectilePz)/
2959 (projectileE - projectilePz) );
2960 G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
2961 G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
2962 G4double targetY = 0.5 *
G4Log( (targetE + targetPz)/(targetE - targetPz) );
2964 #ifdef debugPutOnMassShell
2965 G4cout <<
"decayMomentum2 " << decayMomentum2 <<
G4endl
2966 <<
"\t targetWminus projectileWplus " << targetWminus <<
" " << projectileWplus <<
G4endl
2967 <<
"\t projectileY targetY " << projectileY <<
" " << targetY <<
G4endl;
2970 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
2971 G4Nucleon* aNucleon = involvedNucleons[i];
2972 if ( ! aNucleon )
continue;
2977 G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
2978 G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
2979 if ( isProjectileNucleus ) {
2980 pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*
x);
2981 e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*
x);
2985 #ifdef debugPutOnMassShell
2986 G4cout <<
"i nY pY nY-AY AY " << i <<
" " << nucleonY <<
" " << projectileY <<
G4endl;
2989 if (
std::abs( nucleonY - nucleusY ) > 2 ||
2990 ( isProjectileNucleus && targetY > nucleonY ) ||
2991 ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
3004 const G4bool isProjectileNucleus,
3007 const G4int residualMassNumber,
3008 const G4int numberOfInvolvedNucleons,
3026 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3027 G4Nucleon* aNucleon = involvedNucleons[i];
3028 if ( ! aNucleon )
continue;
3030 residual3Momentum -= tmp.
vect();
3034 G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w *
x );
3035 G4double e = w * x / 2.0 + mt2 / ( 2.0 * w *
x );
3037 if ( isProjectileNucleus ) pz *= -1.0;
3046 G4double residualMt2 =
sqr( residualMass ) +
sqr( residual3Momentum.
x() )
3047 +
sqr( residual3Momentum.
y() );
3049 #ifdef debugPutOnMassShell
3050 G4cout <<
"w residual3Momentum.z() " << w <<
" " << residual3Momentum.
z() <<
G4endl;
3055 if ( residualMassNumber != 0 ) {
3056 residualPz = -w * residual3Momentum.
z() / 2.0 +
3057 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
3058 residualE = w * residual3Momentum.
z() / 2.0 +
3059 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
3061 if ( isProjectileNucleus ) residualPz *= -1.0;
3064 residual4Momentum.
setPx( residual3Momentum.
x() );
3065 residual4Momentum.
setPy( residual3Momentum.
y() );
3066 residual4Momentum.
setPz( residualPz );
3067 residual4Momentum.
setE( residualE );
3076 desc <<
" FTF (Fritiof) Model \n"
3077 <<
"The FTF model is based on the well-known FRITIOF \n"
3078 <<
"model (B. Andersson et al., Nucl. Phys. B281, 289 \n"
3079 <<
"(1987)). Its first program implementation was given\n"
3080 <<
"by B. Nilsson-Almquist and E. Stenlund (Comp. Phys.\n"
3081 <<
"Comm. 43, 387 (1987)). The Fritiof model assumes \n"
3082 <<
"that all hadron-hadron interactions are binary \n"
3083 <<
"reactions, h_1+h_2->h_1'+h_2' where h_1' and h_2' \n"
3084 <<
"are excited states of the hadrons with continuous \n"
3085 <<
"mass spectra. The excited hadrons are considered as\n"
3086 <<
"QCD-strings, and the corresponding LUND-string \n"
3087 <<
"fragmentation model is applied for a simulation of \n"
3088 <<
"their decays. \n"
3089 <<
" The Fritiof model assumes that in the course of \n"
3090 <<
"a hadron-nucleus interaction a string originated \n"
3091 <<
"from the projectile can interact with various intra\n"
3092 <<
"nuclear nucleons and becomes into highly excited \n"
3093 <<
"states. The probability of multiple interactions is\n"
3094 <<
"calculated in the Glauber approximation. A cascading\n"
3095 <<
"of secondary particles was neglected as a rule. Due\n"
3096 <<
"to these, the original Fritiof model fails to des- \n"
3097 <<
"cribe a nuclear destruction and slow particle spectra.\n"
3098 <<
" In order to overcome the difficulties we enlarge\n"
3099 <<
"the model by the reggeon theory inspired model of \n"
3100 <<
"nuclear desctruction (Kh. Abdel-Waged and V.V. Uzhi-\n"
3101 <<
"nsky, Phys. Atom. Nucl. 60, 828 (1997); Yad. Fiz. 60, 925\n"
3102 <<
"(1997)). Momenta of the nucleons ejected from a nuc-\n"
3103 <<
"leus in the reggeon cascading are sampled according\n"
3104 <<
"to a Fermi motion algorithm presented in (EMU-01 \n"
3105 <<
"Collaboration (M.I. Adamovich et al.) Zeit. fur Phys.\n"
3106 <<
"A358, 337 (1997)). \n"
3107 <<
" New features were also added to the Fritiof model\n"
3108 <<
"implemented in Geant4: a simulation of elastic had-\n"
3109 <<
"ron-nucleon scatterings, a simulation of binary \n"
3110 <<
"reactions like NN>NN* in hadron-nucleon interactions,\n"
3111 <<
"a separate simulation of single diffractive and non-\n"
3112 <<
" diffractive events. These allowed to describe after\n"
3113 <<
"model parameter tuning a wide set of experimental \n"