82 for (
G4int i = 0; i < 250; i++ ) {
134 if ( aNucleon )
delete aNucleon;
142 if ( aNucleon )
delete aNucleon;
162 <<
"FTF init Target A Z " << aNucleus.
GetA_asInt()
294 G4cout <<
"FTF PutOnMassShell Success? " << Success <<
G4endl;
306 G4cout <<
"FTF ExciteParticipants Success? " << Success <<
G4endl;
312 G4cout <<
"FTF BuildStrings ";
318 G4cout <<
"FTF BuildStrings " << theStrings <<
" OK" << G4endl
319 <<
"FTF GetResiduals of Nuclei " <<
G4endl;
330 std::vector< G4VSplitableHadron* > primaries;
335 if ( primaries.end() ==
336 std::find( primaries.begin(), primaries.end(), interaction.
GetProjectile() ) ) {
350 if ( aNucleon )
delete aNucleon;
352 NumberOfInvolvedNucleonsOfProjectile = 0;
357 if ( aNucleon )
delete aNucleon;
359 NumberOfInvolvedNucleonsOfTarget = 0;
362 G4cout <<
"End of FTF. Go to fragmentation" << G4endl
363 <<
"To continue - enter 1, to stop - ^C" <<
G4endl;
392 G4cout <<
"G4FTFModel::StoreInvolvedNucleon -------------" <<
G4endl;
407 while ( ( aProjectileNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
430 #ifdef debugReggeonCascade
431 G4cout <<
"G4FTFModel::ReggeonCascade -----------" <<
G4endl
440 for (
G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
468 Neighbour->
Hit( targetSplitable );
476 #ifdef debugReggeonCascade
477 G4cout <<
"Final NumberOfInvolvedNucleonsOfTarget "
487 for (
G4int InvPN = 0; InvPN < InitNINp; InvPN++ ) {
500 while ( ( Neighbour = theProjectileNucleus->
GetNextNucleon() ) ) {
515 Neighbour->
Hit( projectileSplitable );
523 #ifdef debugReggeonCascade
524 G4cout <<
"NumberOfInvolvedNucleonsOfProjectile "
534 G4bool isProjectileNucleus =
false;
536 isProjectileNucleus =
true;
539 #ifdef debugPutOnMassShell
541 if ( isProjectileNucleus ) {
542 G4cout <<
"PutOnMassShell for Nucleus_Nucleus " <<
G4endl;
547 if ( Pprojectile.z() < 0.0 ) {
559 #ifdef debugPutOnMassShell
565 if ( ! isOk )
return false;
574 if ( ! isProjectileNucleus ) {
575 Mprojectile = Pprojectile.mag();
576 M2projectile = Pprojectile.mag2();
577 SumMasses += Mprojectile + 20.0*
MeV;
579 #ifdef debugPutOnMassShell
580 G4cout <<
"Projectile : ";
585 if ( ! isOk )
return false;
592 #ifdef debugPutOnMassShell
593 G4cout <<
"Psum " << Psum/
GeV <<
" GeV" << G4endl <<
"SqrtS " << SqrtS/
GeV <<
" GeV" << G4endl
594 <<
"SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/
GeV <<
" "
595 << PrResidualMass/
GeV <<
" " << TargetResidualMass/
GeV <<
" GeV" <<
G4endl;
598 if ( SqrtS < SumMasses ) {
604 G4double savedSumMasses = SumMasses;
605 if ( isProjectileNucleus ) {
606 SumMasses -= std::sqrt(
sqr( PrResidualMass ) + PprojResidual.perp2() );
608 + PprojResidual.perp2() );
610 SumMasses -= std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.perp2() );
612 + PtargetResidual.perp2() );
614 if ( SqrtS < SumMasses ) {
615 SumMasses = savedSumMasses;
616 if ( isProjectileNucleus ) {
623 if ( isProjectileNucleus ) {
627 #ifdef debugPutOnMassShell
628 if ( isProjectileNucleus ) {
629 G4cout <<
"PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/
GeV <<
" "
632 G4cout <<
"TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/
GeV <<
" "
634 <<
"Sum masses " << SumMasses/
GeV <<
G4endl;
638 if ( isProjectileNucleus && thePrNucleus->
GetMassNumber() != 1 ) {
647 if ( ! isOk )
return false;
658 if ( Ptmp.pz() <= 0.0 ) {
665 if ( isProjectileNucleus ) {
667 YprojectileNucleus = Ptmp.rapidity();
669 Ptmp = toCms*Ptarget;
670 G4double YtargetNucleus = Ptmp.rapidity();
674 if ( isProjectileNucleus ) {
681 #ifdef debugPutOnMassShell
682 if ( isProjectileNucleus ) {
683 G4cout <<
"Y projectileNucleus " << YprojectileNucleus <<
G4endl;
685 G4cout <<
"Y targetNucleus " << YtargetNucleus << G4endl
687 <<
" DcorP DcorT " << DcorP <<
" " << DcorT <<
" AveragePt2 " << AveragePt2 <<
G4endl;
694 G4int NumberOfTries = 0;
696 G4bool OuterSuccess =
true;
698 const G4int maxNumberOfLoops = 1000;
699 G4int loopCounter = 0;
702 const G4int maxNumberOfInnerLoops = 10000;
705 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
711 DcorP *= ScaleFactor;
712 DcorT *= ScaleFactor;
713 AveragePt2 *= ScaleFactor;
715 if ( isProjectileNucleus ) {
718 thePrNucleus, PprojResidual,
726 theTargetNucleus, PtargetResidual,
731 #ifdef debugPutOnMassShell
732 G4cout <<
"SqrtS, Mp+Mt, Mp, Mt " << SqrtS/
GeV <<
" "
733 << ( std::sqrt( M2proj ) + std::sqrt( M2target) )/
GeV <<
" "
734 << std::sqrt( M2proj )/
GeV <<
" " << std::sqrt( M2target )/
GeV <<
G4endl;
737 if ( ! isOk )
return false;
738 }
while ( ( SqrtS < std::sqrt( M2proj ) + std::sqrt( M2target ) ) &&
739 NumberOfTries < maxNumberOfInnerLoops );
740 if ( NumberOfTries >= maxNumberOfInnerLoops ) {
741 #ifdef debugPutOnMassShell
742 G4cout <<
"BAD situation: forced exit of the inner while loop!" <<
G4endl;
746 if ( isProjectileNucleus ) {
747 isOk =
CheckKinematics( S, SqrtS, M2proj, M2target, YprojectileNucleus,
true,
750 WminusTarget, WplusProjectile, OuterSuccess );
755 WminusTarget, WplusProjectile, OuterSuccess );
756 if ( ! isOk )
return false;
757 }
while ( ( ! OuterSuccess ) &&
758 ++loopCounter < maxNumberOfLoops );
759 if ( loopCounter >= maxNumberOfLoops ) {
760 #ifdef debugPutOnMassShell
761 G4cout <<
"BAD situation: forced exit of the while loop!" <<
G4endl;
772 if ( ! isProjectileNucleus ) {
774 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
775 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
776 Pprojectile.setPz( Pzprojectile );
777 Pprojectile.setE( Eprojectile );
779 #ifdef debugPutOnMassShell
780 G4cout <<
"Proj after in CMS " << Pprojectile <<
G4endl;
783 Pprojectile.transform( toLab );
792 #ifdef debugPutOnMassShell
802 #ifdef debugPutOnMassShell
806 if ( ! isOk )
return false;
810 #ifdef debugPutOnMassShell
820 #ifdef debugPutOnMassShell
824 if ( ! isOk )
return false;
828 #ifdef debugPutOnMassShell
841 #ifdef debugBuildString
842 G4cout <<
"G4FTFModel::ExciteParticipants() " <<
G4endl;
845 G4bool Successfull(
true );
847 if ( MaxNumOfInelCollisions > 0 ) {
849 if (
G4UniformRand() < ProbMaxNumber ) MaxNumOfInelCollisions++;
852 MaxNumOfInelCollisions = 1;
855 #ifdef debugBuildString
856 G4cout <<
"MaxNumOfInelCollisions MaxNumOfInelCollisions " << MaxNumOfInelCollisions <<
G4endl;
859 G4int CurrentInteraction( 0 );
864 CurrentInteraction++;
871 #ifdef debugBuildString
872 G4cout << G4endl <<
"Interaction # Status " << CurrentInteraction <<
" "
873 << collision.
GetStatus() << G4endl <<
"Pr* Tr* " << projectile <<
" "
874 << target << G4endl <<
"projectile->GetStatus target->GetStatus "
884 #ifdef debugBuildString
889 G4bool Annihilation =
false;
891 TargetNucleon, Annihilation );
892 if ( ! Result )
continue;
899 #ifdef debugBuildString
900 G4cout <<
"Inelastic interaction" << G4endl
901 <<
"MaxNumOfInelCollisions " << MaxNumOfInelCollisions <<
G4endl;
905 G4bool Annihilation =
false;
907 TargetNucleon, Annihilation );
908 if ( ! Result )
continue;
920 #ifdef debugBuildString
933 #ifdef debugBuildString
934 G4cout <<
"FTF excitation Non Successfull -> Elastic scattering "
940 #ifdef debugBuildString
941 G4cout <<
"Elastic scat. at rejection inelastic scattering" <<
G4endl;
955 #ifdef debugBuildString
964 if ( projectile == NextProjectileNucleon || target == NextTargetNucleon ) {
975 G4bool Annihilation =
true;
977 TargetNucleon, Annihilation );
978 if ( ! Result )
continue;
982 Successfull = Successfull ||
true;
984 #ifdef debugBuildString
985 G4cout <<
"Annihilation successfull. " <<
"*AdditionalString "
986 << AdditionalString <<
G4endl;
996 #ifdef debugBuildString
997 G4cout <<
"----------------------------- Final properties " << G4endl
998 <<
"projectile->GetStatus target->GetStatus " << projectile->
GetStatus()
999 <<
" " << target->
GetStatus() << G4endl <<
"projectile->GetSoftC target->GetSoftC "
1001 << G4endl <<
"ExciteParticipants() Successfull? " << Successfull <<
G4endl;
1019 G4cout <<
"AdjustNucleons ---------------------------------------" <<
G4endl
1023 <<
"Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
1026 <<
"Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
1041 G4double DcorP( 0.0 ), DcorT( 0.0 );
1047 G4double ExcitationEnergyPerWoundedNucleon =
1061 G4cout <<
"case 1, hA prcol=0 trcol=0, AA prcol#0 trcol=0" <<
G4endl;
1091 toCms.rotateZ( -1*Ptmp.phi() );
1092 toCms.rotateY( -1*Ptmp.theta() );
1093 Pprojectile.transform( toCms );
1108 if ( TResidualMassNumber <= 1 ) {
1109 TResidualExcitationEnergy = 0.0;
1113 if ( TResidualMassNumber != 0 ) {
1115 ->
GetIonMass( TResidualCharge, TResidualMassNumber );
1127 if ( ! Annihilation ) {
1128 if ( SqrtS < SumMasses ) {
1131 if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1134 G4cout <<
"TResidualExcitationEnergy " << TResidualExcitationEnergy <<
G4endl;
1137 TResidualExcitationEnergy = SqrtS - SumMasses;
1140 G4cout <<
"TResidualExcitationEnergy " << TResidualExcitationEnergy <<
G4endl;
1148 if ( Annihilation ) {
1151 G4cout <<
"SqrtS < SumMasses - TNucleonMass " << SqrtS <<
" "
1152 << SumMasses - TNucleonMass <<
G4endl;
1155 if ( SqrtS < SumMasses - TNucleonMass ) {
1160 G4cout <<
"SqrtS < SumMasses " << SqrtS <<
" " << SumMasses <<
G4endl;
1163 if ( SqrtS < SumMasses ) {
1164 TNucleonMass = SqrtS - (SumMasses - TNucleonMass) - TResidualExcitationEnergy;
1170 SumMasses = SqrtS - TResidualExcitationEnergy;
1176 G4cout <<
"SqrtS < SumMasses " << SqrtS <<
" " << SumMasses <<
G4endl;
1179 if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1180 TResidualExcitationEnergy = SqrtS - SumMasses;
1192 Ptmp.setPx( 0.0 ); Ptmp.setPy( 0.0 ); Ptmp.setPz( 0.0 );
1199 Pprojectile = Ptmp; Pprojectile.transform( toLab );
1203 Ptmp.setE( TNucleonMass );
1209 Ptarget = Ptmp; Ptarget.transform( toLab );
1223 Ptmp.transform( toLab );
1233 G4double Mprojectile = Pprojectile.mag();
1234 G4double M2projectile = Pprojectile.mag2();
1238 G4double YtargetNucleus = TResidual4Momentum.rapidity();
1240 TResidualMass += TResidualExcitationEnergy;
1249 G4int NumberOfTries( 0 );
1251 G4bool OuterSuccess(
true );
1253 const G4int maxNumberOfLoops = 1000;
1254 G4int loopCounter = 0;
1256 OuterSuccess =
true;
1258 const G4int maxNumberOfTries = 10000;
1263 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1266 DcorT *= ScaleFactor;
1267 AveragePt2 *= ScaleFactor;
1277 G4bool InerSuccess =
true;
1279 const G4int maxNumberOfInnerLoops = 1000;
1280 G4int innerLoopCounter = 0;
1284 PtNucleon =
GaussianPt( AveragePt2, maxPtSquare );
1285 PtResidual = -PtNucleon;
1287 G4double Mtarget = std::sqrt(
sqr( TNucleonMass ) + PtNucleon.mag2() ) +
1288 std::sqrt(
sqr( TResidualMass ) + PtResidual.mag2() );
1289 if ( SqrtS < Mprojectile + Mtarget ) {
1290 InerSuccess =
false;
1295 G4double Xcenter = std::sqrt(
sqr( TNucleonMass ) + PtNucleon.mag2() ) / Mtarget;
1296 XminusNucleon = Xcenter + tmpX.x();
1297 if ( XminusNucleon <= 0.0 || XminusNucleon >= 1.0 ) {
1298 InerSuccess =
false;
1302 XminusResidual = 1.0 - XminusNucleon;
1303 }
while ( ( ! InerSuccess ) &&
1304 ++innerLoopCounter < maxNumberOfInnerLoops );
1305 if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
1307 G4cout <<
"BAD situation: forced exit of the inner while loop!" <<
G4endl;
1313 XminusNucleon = 1.0;
1314 XminusResidual = 1.0;
1318 M2target = (
sqr( TNucleonMass ) + PtNucleon.mag2() ) / XminusNucleon +
1319 (
sqr( TResidualMass ) + PtResidual.mag2() ) / XminusResidual;
1321 }
while ( ( SqrtS < Mprojectile + std::sqrt( M2target) ) &&
1322 ++NumberOfTries < maxNumberOfTries );
1323 if ( NumberOfTries >= maxNumberOfTries ) {
1325 G4cout <<
"BAD situation: forced exit of the intermediate while loop!" <<
G4endl;
1331 - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
1333 WminusTarget = ( S - M2projectile + M2target + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS;
1334 WplusProjectile = SqrtS - M2target / WminusTarget;
1336 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
1337 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
1338 G4double Yprojectile = 0.5 *
G4Log( (Eprojectile + Pzprojectile) /
1339 (Eprojectile - Pzprojectile) );
1342 G4cout <<
"DecayMomentum2 " << DecayMomentum2 << G4endl
1343 <<
"WminusTarget WplusProjectile " << WminusTarget <<
" " << WplusProjectile
1344 << G4endl <<
"Yprojectile " << Yprojectile <<
G4endl;
1347 G4double Mt2 =
sqr( TNucleonMass ) + PtNucleon.mag2();
1348 G4double Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1349 G4double E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1353 G4cout <<
"YtN Ytr YtN-Ytr " <<
" " << YtargetNucleon <<
" " << YtargetNucleus <<
" "
1354 << YtargetNucleon - YtargetNucleus << G4endl
1355 <<
"YtN Ypr YtN-Ypr " <<
" " << YtargetNucleon <<
" " << Yprojectile
1356 <<
" " << YtargetNucleon - Yprojectile <<
G4endl;
1359 if ( std::abs( YtargetNucleon - YtargetNucleus ) > 2 || Yprojectile < YtargetNucleon ) {
1360 OuterSuccess =
false;
1364 }
while ( ( ! OuterSuccess ) &&
1365 ++loopCounter < maxNumberOfLoops );
1366 if ( loopCounter >= maxNumberOfLoops ) {
1368 G4cout <<
"BAD situation: forced exit of the while loop!" <<
G4endl;
1373 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
1374 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
1375 Pprojectile.setPz( Pzprojectile ); Pprojectile.setE( Eprojectile );
1378 G4cout <<
"Proj after in CMS " << Pprojectile <<
G4endl;
1381 Pprojectile.transform( toLab );
1389 G4double Mt2 =
sqr( TNucleonMass ) + PtNucleon.mag2();
1390 G4double Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1391 G4double E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1393 Ptarget.setPx( PtNucleon.x() ); Ptarget.setPy( PtNucleon.y() );
1394 Ptarget.setPz( Pz ); Ptarget.setE( E );
1395 Ptarget.transform( toLab );
1408 G4cout <<
"TargetResidualMassNumber TargetResidualCharge TargetResidualExcitationEnergy "
1414 Mt2 =
sqr( TResidualMass ) + PtResidual.mag2();
1415 Pz = -WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
1416 E = WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
1418 TargetResidual4Momentum.setPx( PtResidual.x() );
1419 TargetResidual4Momentum.setPy( PtResidual.y() );
1420 TargetResidual4Momentum.setPz( Pz );
1421 TargetResidual4Momentum.setE( E );
1424 G4cout <<
"New Residu " << TargetResidual4Momentum <<
" CMS" <<
G4endl;
1427 TargetResidual4Momentum.transform( toLab );
1430 G4cout <<
"New Residu " << TargetResidual4Momentum <<
" Lab" <<
G4endl;
1468 toCms.rotateZ( -1*Ptmp.phi() );
1469 toCms.rotateY( -1*Ptmp.theta() );
1472 Pprojectile.transform( toCms );
1484 if ( TResidualMassNumber <= 1 ) {
1485 TResidualExcitationEnergy = 0.0;
1489 if ( TResidualMassNumber != 0 ) {
1491 ->
GetIonMass( TResidualCharge , TResidualMassNumber );
1497 TNucleonMass + TResidualMass;
1500 G4cout <<
"SelectedTN.mag() PNMass + PResidualMass "
1502 << TNucleonMass <<
" " << TResidualMass <<
G4endl;
1507 if ( ! Annihilation ) {
1508 if ( SqrtS < SumMasses ) {
1511 if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1512 TResidualExcitationEnergy = SqrtS - SumMasses;
1518 if ( Annihilation ) {
1519 if ( SqrtS < SumMasses - TNucleonMass ) {
1522 if ( SqrtS < SumMasses ) {
1523 TNucleonMass = SqrtS - (SumMasses - TNucleonMass);
1525 TResidualExcitationEnergy = 0.0;
1529 if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1530 TResidualExcitationEnergy = SqrtS - SumMasses;
1542 Ptmp.setPx( 0.0 ); Ptmp.setPy( 0.0 ); Ptmp.setPz( 0.0 );
1543 Ptmp.setE( SelectedTargetNucleon->
Get4Momentum().mag() );
1544 Ptarget = Ptmp; Ptarget.transform( toLab );
1548 Ptmp.setE( TNucleonMass );
1549 Pprojectile = Ptmp; Pprojectile.transform( toLab );
1558 Ptmp.transform( toLab );
1565 G4double M2target = Ptarget.mag2();
1568 G4double YprojectileNucleus = TResidual4Momentum.rapidity();
1570 TResidualMass += TResidualExcitationEnergy;
1579 G4int NumberOfTries( 0 );
1581 G4bool OuterSuccess(
true );
1583 const G4int maxNumberOfLoops = 1000;
1584 G4int loopCounter = 0;
1587 OuterSuccess =
true;
1588 const G4int maxNumberOfTries = 10000;
1593 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1596 DcorP *= ScaleFactor;
1597 AveragePt2 *= ScaleFactor;
1605 PtNucleon =
GaussianPt( AveragePt2, maxPtSquare );
1609 PtResidual = -PtNucleon;
1611 G4double Mprojectile = std::sqrt(
sqr( TNucleonMass ) + PtNucleon.mag2() ) +
1612 std::sqrt(
sqr( TResidualMass ) + PtResidual.mag2() );
1615 G4cout <<
"SqrtS < Mtarget + Mprojectile " << SqrtS <<
" " << Mtarget
1616 <<
" " << Mprojectile <<
" " << Mtarget + Mprojectile <<
G4endl;
1619 M2projectile =
sqr( Mprojectile );
1620 if ( SqrtS < Mtarget + Mprojectile ) {
1621 OuterSuccess =
false;
1625 G4double Xcenter = std::sqrt(
sqr( TNucleonMass ) + PtNucleon.mag2() ) / Mprojectile;
1627 G4bool InerSuccess =
true;
1629 const G4int maxNumberOfInnerLoops = 1000;
1630 G4int innerLoopCounter = 0;
1634 XplusNucleon = Xcenter + tmpX.x();
1635 if ( XplusNucleon <= 0.0 || XplusNucleon >= 1.0 ) {
1636 InerSuccess =
false;
1639 XplusResidual = 1.0 - XplusNucleon;
1640 }
while ( ( ! InerSuccess ) &&
1641 ++innerLoopCounter < maxNumberOfInnerLoops );
1642 if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
1644 G4cout <<
"BAD situation: forced exit of the inner while loop!" <<
G4endl;
1651 XplusResidual = 1.0;
1656 G4cout <<
"TNucleonMass PtNucleon XplusNucleon " << TNucleonMass <<
" " << PtNucleon
1657 <<
" " << XplusNucleon << G4endl
1658 <<
"TResidualMass PtResidual XplusResidual " << TResidualMass <<
" " << PtResidual
1659 <<
" " << XplusResidual <<
G4endl;
1662 M2projectile = (
sqr( TNucleonMass ) + PtNucleon.mag2() ) / XplusNucleon +
1663 (
sqr( TResidualMass ) + PtResidual.mag2() ) / XplusResidual;
1666 G4cout <<
"SqrtS < Mtarget + std::sqrt(M2projectile) " << SqrtS <<
" " << Mtarget
1667 <<
" " << std::sqrt( M2projectile ) <<
" " << Mtarget + std::sqrt( M2projectile )
1671 }
while ( ( SqrtS < Mtarget + std::sqrt( M2projectile ) ) &&
1672 ++NumberOfTries < maxNumberOfTries );
1673 if ( NumberOfTries >= maxNumberOfTries ) {
1675 G4cout <<
"BAD situation: forced exit of the intermediate while loop!" <<
G4endl;
1681 - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
1683 WplusProjectile = ( S + M2projectile - M2target + std::sqrt( DecayMomentum2 ) )/2.0/SqrtS;
1684 WminusTarget = SqrtS - M2projectile/WplusProjectile;
1686 G4double Pztarget = -WminusTarget/2.0 + M2target/2.0/WminusTarget;
1687 G4double Etarget = WminusTarget/2.0 + M2target/2.0/WminusTarget;
1688 G4double Ytarget = 0.5 *
G4Log( (Etarget + Pztarget)/(Etarget - Pztarget) );
1691 G4cout <<
"DecayMomentum2 " << DecayMomentum2 << G4endl
1692 <<
"WminusTarget WplusProjectile " << WminusTarget <<
" " << WplusProjectile
1693 << G4endl <<
"YtargetNucleon " << Ytarget <<
G4endl;
1696 G4double Mt2 =
sqr( TNucleonMass ) + PtNucleon.mag2();
1697 G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
1698 G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
1699 G4double YprojectileNucleon = 0.5 *
G4Log( (E + Pz)/(E - Pz) );
1702 G4cout <<
"YpN Ypr YpN-Ypr " <<
" " << YprojectileNucleon <<
" " << YprojectileNucleus
1703 <<
" " << YprojectileNucleon - YprojectileNucleus << G4endl
1704 <<
"YpN Ytr YpN-Ytr " <<
" " << YprojectileNucleon <<
" " << Ytarget
1705 <<
" " << YprojectileNucleon - Ytarget <<
G4endl;
1708 if ( std::abs( YprojectileNucleon - YprojectileNucleus ) > 2 ||
1709 Ytarget > YprojectileNucleon ) {
1710 OuterSuccess =
false;
1714 }
while ( ( ! OuterSuccess ) &&
1715 ++loopCounter < maxNumberOfLoops );
1716 if ( loopCounter >= maxNumberOfLoops ) {
1718 G4cout <<
"BAD situation: forced exit of the while loop!" <<
G4endl;
1724 G4double Pztarget = -WminusTarget/2.0 + M2target/2.0/WminusTarget;
1725 G4double Etarget = WminusTarget/2.0 + M2target/2.0/WminusTarget;
1726 Ptarget.setPz( Pztarget ); Ptarget.setE( Etarget );
1727 Ptarget.transform( toLab );
1735 G4double Mt2 =
sqr( TNucleonMass ) + PtNucleon.mag2();
1736 G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
1737 G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
1738 Pprojectile.setPx( PtNucleon.x() ); Pprojectile.setPy( PtNucleon.y() );
1739 Pprojectile.setPz( Pz ); Pprojectile.setE( E );
1740 Pprojectile.transform( toLab );
1744 G4cout <<
"Proj after in Lab " << Pprojectile <<
G4endl;
1753 Mt2 =
sqr( TResidualMass ) + PtResidual.mag2();
1754 Pz = WplusProjectile*XplusResidual/2.0 - Mt2/(2.0*WplusProjectile*XplusResidual);
1755 E = WplusProjectile*XplusResidual/2.0 + Mt2/(2.0*WplusProjectile*XplusResidual);
1756 ProjectileResidual4Momentum.setPx( PtResidual.x() );
1757 ProjectileResidual4Momentum.setPy( PtResidual.y() );
1758 ProjectileResidual4Momentum.setPz( Pz );
1759 ProjectileResidual4Momentum.setE( E );
1760 ProjectileResidual4Momentum.transform( toLab );
1780 <<
"ProjectileResidualMassNumber ProjectileResidualCharge "
1792 toCms.rotateZ( -1*Ptmp.phi() );
1793 toCms.rotateY( -1*Ptmp.theta() );
1795 Pprojectile.transform( toCms );
1808 if ( PResidualMassNumber <= 1 ) {
1809 PResidualExcitationEnergy = 0.0;
1813 if ( PResidualMassNumber != 0 ) {
1815 ->
GetIonMass( PResidualCharge, PResidualMassNumber );
1827 if ( TResidualMassNumber <= 1 ) {
1828 TResidualExcitationEnergy = 0.0;
1831 if ( TResidualMassNumber != 0 ) {
1833 ->
GetIonMass( TResidualCharge, TResidualMassNumber );
1838 G4double SumMasses = PNucleonMass + PResidualMass + TNucleonMass + TResidualMass;
1841 G4cout <<
"PNucleonMass PResidualMass TNucleonMass TResidualMass " << PNucleonMass
1842 <<
" " << PResidualMass <<
" " << TNucleonMass <<
" " << TResidualMass << G4endl
1843 <<
"PResidualExcitationEnergy " << PResidualExcitationEnergy << G4endl
1844 <<
"TResidualExcitationEnergy " << TResidualExcitationEnergy <<
G4endl;
1849 if ( ! Annihilation ) {
1852 G4cout <<
"SqrtS < SumMasses " << SqrtS <<
" " << SumMasses <<
G4endl;
1855 if ( SqrtS < SumMasses ) {
1860 G4cout <<
"SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy "
1861 << SqrtS <<
" " << SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy
1865 if ( SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy ) {
1868 if ( PResidualExcitationEnergy <= 0.0 ) {
1869 TResidualExcitationEnergy = SqrtS - SumMasses;
1870 }
else if ( TResidualExcitationEnergy <= 0.0 ) {
1871 PResidualExcitationEnergy = SqrtS - SumMasses;
1873 G4double Fraction = (SqrtS - SumMasses) /
1874 (PResidualExcitationEnergy + TResidualExcitationEnergy);
1875 PResidualExcitationEnergy *= Fraction;
1876 TResidualExcitationEnergy *= Fraction;
1885 if ( Annihilation ) {
1886 if ( SqrtS < SumMasses - TNucleonMass ) {
1889 if ( SqrtS < SumMasses ) {
1891 TNucleonMass = SqrtS - (SumMasses - TNucleonMass);
1893 TResidualExcitationEnergy = 0.0;
1895 if ( SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy ) {
1897 if ( PResidualExcitationEnergy <= 0.0 ) {
1898 TResidualExcitationEnergy = SqrtS - SumMasses;
1899 }
else if ( TResidualExcitationEnergy <= 0.0 ) {
1900 PResidualExcitationEnergy = SqrtS - SumMasses;
1902 G4double Fraction = (SqrtS - SumMasses) /
1903 (PResidualExcitationEnergy + TResidualExcitationEnergy);
1904 PResidualExcitationEnergy *= Fraction;
1905 TResidualExcitationEnergy *= Fraction;
1913 Ptmp.setPx( 0.0 ); Ptmp.setPy( 0.0 ); Ptmp.setPz( 0.0 );
1914 Ptmp.setE( PNucleonMass );
1915 Pprojectile = Ptmp; Pprojectile.transform( toLab );
1924 Ptmp.transform( toLab );
1928 Ptmp.setPx( 0.0 ); Ptmp.setPy( 0.0 ); Ptmp.setPz( 0.0 );
1929 Ptmp.setE( TNucleonMass );
1930 Ptarget = Ptmp; Ptarget.transform( toLab );
1939 Ptmp.transform( toLab );
1940 TargetResidual4Momentum = Ptmp;
1946 G4double YprojectileNucleus = PResidual4Momentum.rapidity();
1949 G4cout <<
"YprojectileNucleus XcenterP " << YprojectileNucleus <<
G4endl;
1953 G4double YtargetNucleus = TResidual4Momentum.rapidity();
1955 PResidualMass += PResidualExcitationEnergy;
1956 TResidualMass += TResidualExcitationEnergy;
1973 G4int NumberOfTries( 0 );
1975 G4bool OuterSuccess(
true );
1977 const G4int maxNumberOfLoops = 1000;
1978 G4int loopCounter = 0;
1981 OuterSuccess =
true;
1982 const G4int maxNumberOfTries = 10000;
1987 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1990 DcorP *= ScaleFactor;
1991 DcorT *= ScaleFactor;
1992 AveragePt2 *= ScaleFactor;
2000 PtNucleonP =
GaussianPt( AveragePt2, maxPtSquare );
2004 PtResidualP = -PtNucleonP;
2007 PtNucleonT =
GaussianPt( AveragePt2, maxPtSquare );
2011 PtResidualT = -PtNucleonT;
2013 G4double Mprojectile = std::sqrt(
sqr( PNucleonMass ) + PtNucleonP.mag2() ) +
2014 std::sqrt(
sqr( PResidualMass ) + PtResidualP.mag2() );
2015 M2projectile =
sqr( Mprojectile );
2017 G4double Mtarget = std::sqrt(
sqr( TNucleonMass ) + PtNucleonT.mag2() ) +
2018 std::sqrt(
sqr( TResidualMass ) + PtResidualT.mag2() );
2019 M2target =
sqr( Mtarget );
2021 if ( SqrtS < Mprojectile + Mtarget ) {
2022 OuterSuccess =
false;
2026 G4bool InerSuccess =
true;
2029 const G4int maxNumberOfInnerLoops = 1000;
2030 G4int innerLoopCounter = 0;
2034 G4double XcenterP = std::sqrt(
sqr( PNucleonMass ) + PtNucleonP.mag2() ) / Mprojectile;
2035 XplusNucleon = XcenterP + tmpX.x();
2042 if ( XplusNucleon <= 0.0 || XplusNucleon >= 1.0 ) {
2043 InerSuccess =
false;
2046 XplusResidual = 1.0 - XplusNucleon;
2047 }
while ( ( ! InerSuccess ) &&
2048 ++innerLoopCounter < maxNumberOfInnerLoops );
2049 if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
2051 G4cout <<
"BAD situation: forced exit of the first inner while loop!" <<
G4endl;
2064 XplusResidual = 1.0;
2069 const G4int maxNumberOfInnerLoops = 1000;
2070 G4int innerLoopCounter = 0;
2075 G4double XcenterT = std::sqrt(
sqr( TNucleonMass ) + PtNucleonT.mag2() ) / Mtarget;
2076 XminusNucleon = XcenterT + tmpX.x();
2077 if ( XminusNucleon <= 0.0 || XminusNucleon >= 1.0 ) {
2078 InerSuccess =
false;
2081 XminusResidual = 1.0 - XminusNucleon;
2082 }
while ( ( ! InerSuccess ) &&
2083 ++innerLoopCounter < maxNumberOfInnerLoops );
2084 if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
2086 G4cout <<
"BAD situation: forced exit of the second inner while loop!" <<
G4endl;
2091 XminusNucleon = 1.0;
2092 XminusResidual = 1.0;
2096 G4cout <<
"PtNucleonP " << PtNucleonP <<
" " << PtResidualP << G4endl
2097 <<
"XplusNucleon XplusResidual " << XplusNucleon <<
" " << XplusResidual << G4endl
2098 <<
"PtNucleonT " << PtNucleonT <<
" " << PtResidualT << G4endl
2099 <<
"XminusNucleon XminusResidual " << XminusNucleon <<
" " << XminusResidual
2103 M2projectile = (
sqr( PNucleonMass ) + PtNucleonP.mag2() ) / XplusNucleon +
2104 (
sqr( PResidualMass) + PtResidualP.mag2() ) / XplusResidual;
2105 M2target = (
sqr( TNucleonMass ) + PtNucleonT.mag2() ) / XminusNucleon +
2106 (
sqr( TResidualMass ) + PtResidualT.mag2() ) / XminusResidual;
2108 }
while ( ( SqrtS < std::sqrt( M2projectile ) + std::sqrt( M2target ) ) &&
2109 ++NumberOfTries < maxNumberOfTries );
2110 if ( NumberOfTries >= maxNumberOfTries ) {
2112 G4cout <<
"BAD situation: forced exit of the intermediate while loop!" <<
G4endl;
2118 - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
2120 WplusProjectile = ( S + M2projectile - M2target + std::sqrt( DecayMomentum2 ) )/2.0/SqrtS;
2121 WminusTarget = SqrtS - M2projectile/WplusProjectile;
2123 G4double Mt2 =
sqr( PNucleonMass ) + PtNucleonP.mag2();
2124 G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
2125 G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
2126 G4double YprojectileNucleon = 0.5 *
G4Log( (E + Pz)/(E - Pz) );
2128 Mt2 =
sqr( TNucleonMass ) + PtNucleonT.mag2();
2129 Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2130 E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2133 if ( std::abs( YtargetNucleon - YtargetNucleus ) > 2 ||
2134 std::abs( YprojectileNucleon - YprojectileNucleus ) > 2 ||
2135 YprojectileNucleon < YtargetNucleon ) {
2136 OuterSuccess =
false;
2140 }
while ( ( ! OuterSuccess ) &&
2141 ++loopCounter < maxNumberOfLoops );
2142 if ( loopCounter >= maxNumberOfLoops ) {
2144 G4cout <<
"BAD situation: forced exit of the while loop!" <<
G4endl;
2153 G4double Mt2 =
sqr( PNucleonMass ) + PtNucleonP.mag2();
2154 G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
2155 G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
2157 Pprojectile.setPx( PtNucleonP.x() ); Pprojectile.setPy( PtNucleonP.y() );
2158 Pprojectile.setPz( Pz ); Pprojectile.setE( E );
2159 Pprojectile.transform( toLab );
2168 G4cout <<
"PResidualMass PtResidualP " << PResidualMass <<
" " << PtResidualP <<
G4endl;
2172 Mt2 =
sqr( PResidualMass ) + PtResidualP.mag2();
2173 Pz = WplusProjectile*XplusResidual/2.0 - Mt2/(2.0*WplusProjectile*XplusResidual);
2174 E = WplusProjectile*XplusResidual/2.0 + Mt2/(2.0*WplusProjectile*XplusResidual);
2175 ProjectileResidual4Momentum.setPx( PtResidualP.x() );
2176 ProjectileResidual4Momentum.setPy( PtResidualP.y() );
2177 ProjectileResidual4Momentum.setPz( Pz );
2178 ProjectileResidual4Momentum.setE( E );
2179 ProjectileResidual4Momentum.transform( toLab );
2185 G4cout <<
"Pr N R " << Pprojectile << G4endl <<
" "
2186 << ProjectileResidual4Momentum <<
G4endl;
2189 Mt2 =
sqr( TNucleonMass ) + PtNucleonT.mag2();
2190 Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2191 E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2193 Ptarget.setPx( PtNucleonT.x() ); Ptarget.setPy( PtNucleonT.y() );
2194 Ptarget.setPz( Pz ); Ptarget.setE( E );
2195 Ptarget.transform( toLab );
2204 Mt2 =
sqr( TResidualMass ) + PtResidualT.mag2();
2205 Pz = -WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
2206 E = WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
2208 TargetResidual4Momentum.setPx( PtResidualT.x() );
2209 TargetResidual4Momentum.setPy( PtResidualT.y() );
2210 TargetResidual4Momentum.setPz( Pz );
2211 TargetResidual4Momentum.setE( E) ;
2212 TargetResidual4Momentum.transform( toLab );
2218 G4cout <<
"Tr N R " << Ptarget << G4endl <<
" " << TargetResidual4Momentum <<
G4endl;
2240 std::vector< G4VSplitableHadron* > primaries;
2246 if ( primaries.end() == std::find( primaries.begin(), primaries.end(),
2253 #ifdef debugBuildString
2255 <<
"Number of projectile strings " << primaries.size() <<
G4endl;
2258 for (
unsigned int ahadron = 0; ahadron < primaries.size(); ahadron++ ) {
2259 G4bool isProjectile(
true );
2262 FirstString = 0; SecondString = 0;
2263 if ( primaries[ahadron]->GetStatus() <= 1 )
2268 else if(primaries[ahadron]->GetStatus() == 2)
2272 primaries[ahadron]->GetDefinition(),
2273 primaries[ahadron]->GetTimeOfCreation(),
2274 primaries[ahadron]->GetPosition(),
2278 else {
G4cout<<
"Something wrong in FTF Model Build String" <<
G4endl;}
2280 if ( FirstString != 0 ) strings->push_back( FirstString );
2281 if ( SecondString != 0 ) strings->push_back( SecondString );
2283 #ifdef debugBuildString
2284 G4cout <<
"FirstString & SecondString? " << FirstString <<
" " << SecondString <<
G4endl;
2294 #ifdef debugBuildString
2297 G4cout <<
"Check 1 string " << strings->operator[](0)->GetRightParton()->GetPDGcode()
2298 <<
" " << strings->operator[](0)->GetLeftParton()->GetPDGcode() <<
G4endl <<
G4endl;
2307 #ifdef debugBuildString
2308 G4cout <<
"Building of projectile-like strings" <<
G4endl;
2311 G4bool isProjectile =
true;
2314 #ifdef debugBuildString
2315 G4cout <<
"Nucleon #, status, intCount " << ahadron <<
" "
2324 #ifdef debugBuildString
2325 G4cout << G4endl <<
"ahadron aProjectile Status " << ahadron <<
" " << aProjectile
2331 #ifdef debugBuildString
2332 G4cout <<
"Case1 aProjectile->GetStatus() == 0 " <<
G4endl;
2335 FirstString = 0; SecondString = 0;
2339 if ( FirstString != 0 ) strings->push_back( FirstString );
2340 if ( SecondString != 0 ) strings->push_back( SecondString );
2344 #ifdef debugBuildString
2345 G4cout <<
"Case2 aProjectile->GetStatus() !=0 St==1 SoftCol!=0" <<
G4endl;
2348 FirstString = 0; SecondString = 0;
2352 if ( FirstString != 0 ) strings->push_back( FirstString );
2353 if ( SecondString != 0 ) strings->push_back( SecondString );
2360 #ifdef debugBuildString
2361 G4cout <<
"Case3 aProjectile->GetStatus() !=0 St==1 SoftCol==0" <<
G4endl;
2364 FirstString = 0; SecondString = 0;
2368 if ( FirstString != 0 ) strings->push_back( FirstString );
2369 if ( SecondString != 0 ) strings->push_back( SecondString );
2371 #ifdef debugBuildString
2372 G4cout <<
" Strings are built for nucleon marked for an interaction, but"
2373 <<
" the interaction was skipped." <<
G4endl;
2379 #ifdef debugBuildString
2380 G4cout <<
"Case4 aProjectile->GetStatus() !=0 St==2 " <<
G4endl;
2383 FirstString = 0; SecondString = 0;
2387 if ( FirstString != 0 ) strings->push_back( FirstString );
2388 if ( SecondString != 0 ) strings->push_back( SecondString );
2390 #ifdef debugBuildString
2391 G4cout <<
" Strings are build for involved nucleon." <<
G4endl;
2396 #ifdef debugBuildString
2403 #ifdef debugBuildString
2411 #ifdef debugBuildString
2415 G4bool isProjectile =
false;
2419 #ifdef debugBuildString
2420 G4cout <<
"Nucleon #, status, intCount " << aNucleon <<
" " << ahadron <<
" "
2425 FirstString = 0 ; SecondString = 0;
2428 if ( FirstString != 0 ) strings->push_back( FirstString );
2429 if ( SecondString != 0 ) strings->push_back( SecondString );
2431 #ifdef debugBuildString
2437 FirstString = 0; SecondString = 0;
2440 if ( FirstString != 0 ) strings->push_back( FirstString );
2441 if ( SecondString != 0 ) strings->push_back( SecondString );
2443 #ifdef debugBuildString
2444 G4cout <<
" 2 case A string is build, nucleon was excited." <<
G4endl;
2452 FirstString = 0; SecondString = 0;
2456 if(SecondString == 0)
2468 if ( FirstString != 0 ) strings->push_back( FirstString );
2469 if ( SecondString != 0 ) strings->push_back( SecondString );
2471 #ifdef debugBuildString
2483 #ifdef debugBuildString
2487 }
else if(( aNucleon->
GetStatus() == 2 )||
2489 FirstString = 0; SecondString = 0;
2493 if(SecondString == 0)
2505 if ( FirstString != 0 ) strings->push_back( FirstString );
2506 if ( SecondString != 0 ) strings->push_back( SecondString );
2508 #ifdef debugBuildString
2514 #ifdef debugBuildString
2521 #ifdef debugBuildString
2526 isProjectile =
true;
2530 FirstString = 0; SecondString = 0;
2533 if ( FirstString != 0 ) strings->push_back( FirstString );
2534 if ( SecondString != 0 ) strings->push_back( SecondString );
2553 #ifdef debugFTFmodel
2554 G4cout <<
"GetResiduals(): HighEnergyInter? GetProjectileNucleus()?"
2560 #ifdef debugFTFmodel
2572 #ifdef debugFTFmodel
2596 residualMomentum +=tmp;
2610 G4double E=std::sqrt(tmp.vect().mag2()+
2618 const G4int maxNumberOfLoops = 1000;
2619 G4int loopCounter = 0;
2636 if(SumMasses > Mass) {Chigh=
C;}
2639 }
while( (Chigh-Clow > 0.01) &&
2640 ++loopCounter < maxNumberOfLoops );
2641 if ( loopCounter >= maxNumberOfLoops ) {
2642 #ifdef debugFTFmodel
2643 G4cout <<
"BAD situation: forced exit of the first while loop in G4FTFModel::GetResidual" << G4endl
2644 <<
"\t return immediately from the method!" <<
G4endl;
2654 G4double E=std::sqrt(tmp.vect().mag2()+
2656 tmp.setE(E); tmp.boost(-bstToCM);
2665 #ifdef debugFTFmodel
2667 << G4endl <<
"ProjectileResidualExcitationEnergy ProjectileResidual4Momentum "
2679 #ifdef debugFTFmodel
2699 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2703 residualMomentum +=tmp;
2714 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2717 G4double E=std::sqrt(tmp.vect().mag2()+
2725 const G4int maxNumberOfLoops = 1000;
2726 G4int loopCounter = 0;
2734 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2743 if(SumMasses > Mass) {Chigh=
C;}
2746 }
while( (Chigh-Clow > 0.01) &&
2747 ++loopCounter < maxNumberOfLoops );
2748 if ( loopCounter >= maxNumberOfLoops ) {
2749 #ifdef debugFTFmodel
2750 G4cout <<
"BAD situation: forced exit of the second while loop in G4FTFModel::GetResidual" << G4endl
2751 <<
"\t return immediately from the method!" <<
G4endl;
2758 while ( ( aNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
2761 G4double E=std::sqrt(tmp.vect().mag2()+
2763 tmp.setE(E); tmp.boost(-bstToCM);
2769 #ifdef debugFTFmodel
2775 #ifdef debugFTFmodel
2776 G4cout <<
"Low energy interaction: Target nucleus --------------" << G4endl
2777 <<
"Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
2782 G4int NumberOfTargetParticipant( 0 );
2792 if ( NumberOfTargetParticipant != 0 ) {
2805 delete targetSplitable;
2806 targetSplitable = 0;
2807 aNucleon->
Hit( targetSplitable );
2812 #ifdef debugFTFmodel
2813 G4cout <<
"NumberOfTargetParticipant " << NumberOfTargetParticipant << G4endl
2819 #ifdef debugFTFmodel
2820 G4cout <<
"Low energy interaction: Projectile nucleus --------------" << G4endl
2821 <<
"Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
2826 G4int NumberOfProjectileParticipant( 0 );
2831 NumberOfProjectileParticipant++;
2834 #ifdef debugFTFmodel
2838 DeltaExcitationE = 0.0;
2841 if ( NumberOfProjectileParticipant != 0 ) {
2843 G4double( NumberOfProjectileParticipant );
2845 G4double( NumberOfProjectileParticipant );
2857 delete projectileSplitable;
2858 projectileSplitable = 0;
2859 aNucleon->
Hit( projectileSplitable );
2864 #ifdef debugFTFmodel
2865 G4cout <<
"NumberOfProjectileParticipant " << NumberOfProjectileParticipant << G4endl
2871 #ifdef debugFTFmodel
2872 G4cout <<
"End GetResiduals -----------------" <<
G4endl;
2884 if ( AveragePt2 <= 0.0 ) {
2888 (
G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
2893 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
2904 G4double& residualExcitationEnergy,
2906 G4int& residualMassNumber,
2907 G4int& residualCharge ) {
2919 if ( ! nucleus )
return false;
2921 G4double ExcitationEnergyPerWoundedNucleon =
2944 sumMasses += 20.0*
MeV;
2947 residualExcitationEnergy += -ExcitationEnergyPerWoundedNucleon*
2949 residualMassNumber--;
2956 #ifdef debugPutOnMassShell
2957 G4cout <<
"ExcitationEnergyPerWoundedNucleon " << ExcitationEnergyPerWoundedNucleon <<
G4endl
2958 <<
"\t Residual Charge, MassNumber " << residualCharge <<
" " << residualMassNumber
2959 <<
G4endl <<
"\t Initial Momentum " << nucleusMomentum
2960 <<
G4endl <<
"\t Residual Momentum " << residualMomentum <<
G4endl;
2962 residualMomentum.setPz( 0.0 );
2963 residualMomentum.setE( 0.0 );
2964 if ( residualMassNumber == 0 ) {
2966 residualExcitationEnergy = 0.0;
2969 GetIonMass( residualCharge, residualMassNumber );
2970 if ( residualMassNumber == 1 ) {
2971 residualExcitationEnergy = 0.0;
2974 sumMasses += std::sqrt(
sqr( residualMass ) + residualMomentum.perp2() );
2983 const G4int numberOfInvolvedNucleons,
3001 if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 )
return false;
3005 const G4double probDeltaIsobar = 0.10;
3007 G4int maxNumberOfDeltas =
G4int( (sqrtS - sumMasses)/(400.0*
MeV) );
3008 G4int numberOfDeltas = 0;
3010 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3013 if (
G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
3015 if ( ! involvedNucleons[i] )
continue;
3022 G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4;
3031 if ( sqrtS < sumMasses + massDelta - massNuc ) {
3035 sumMasses += ( massDelta - massNuc );
3054 const G4int residualMassNumber,
3055 const G4int numberOfInvolvedNucleons,
3069 if ( ! nucleus )
return false;
3071 if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
3079 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3080 G4Nucleon* aNucleon = involvedNucleons[i];
3081 if ( ! aNucleon )
continue;
3085 const G4int maxNumberOfLoops = 1000;
3086 G4int loopCounter = 0;
3093 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3094 G4Nucleon* aNucleon = involvedNucleons[i];
3095 if ( ! aNucleon )
continue;
3101 if ( x < 0.0 || x > 1.0 ) {
3111 if ( xSum < 0.0 || xSum > 1.0 ) success =
false;
3113 if ( ! success )
continue;
3115 G4double deltaPx = ( ptSum.x() - pResidual.x() ) / numberOfInvolvedNucleons;
3116 G4double deltaPy = ( ptSum.y() - pResidual.y() ) / numberOfInvolvedNucleons;
3118 if ( residualMassNumber == 0 ) {
3119 delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
3126 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3127 G4Nucleon* aNucleon = involvedNucleons[i];
3128 if ( ! aNucleon )
continue;
3131 if ( residualMassNumber == 0 ) {
3132 if ( x <= 0.0 || x > 1.0 ) {
3137 if ( x <= 0.0 || x > 1.0 || xSum <= 0.0 || xSum > 1.0 ) {
3145 +
sqr( px ) +
sqr( py ) ) / x;
3150 if ( success && residualMassNumber != 0 ) {
3151 mass2 += (
sqr( residualMass ) + pResidual.perp2() ) / xSum;
3154 #ifdef debugPutOnMassShell
3158 }
while ( ( ! success ) &&
3159 ++loopCounter < maxNumberOfLoops );
3160 if ( loopCounter >= maxNumberOfLoops ) {
3176 const G4bool isProjectileNucleus,
3177 const G4int numberOfInvolvedNucleons,
3192 G4double decayMomentum2 =
sqr( sValue ) +
sqr( projectileMass2 ) +
sqr( targetMass2 )
3193 - 2.0*sValue*projectileMass2 - 2.0*sValue*targetMass2
3194 - 2.0*projectileMass2*targetMass2;
3195 targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
3197 projectileWplus = sqrtS - targetMass2/targetWminus;
3198 G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
3199 G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
3200 G4double projectileY = 0.5 *
G4Log( (projectileE + projectilePz)/
3201 (projectileE - projectilePz) );
3202 G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
3203 G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
3204 G4double targetY = 0.5 *
G4Log( (targetE + targetPz)/(targetE - targetPz) );
3206 #ifdef debugPutOnMassShell
3207 G4cout <<
"decayMomentum2 " << decayMomentum2 <<
G4endl
3208 <<
"\t targetWminus projectileWplus " << targetWminus <<
" " << projectileWplus <<
G4endl
3209 <<
"\t projectileY targetY " << projectileY <<
" " << targetY <<
G4endl;
3212 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3213 G4Nucleon* aNucleon = involvedNucleons[i];
3214 if ( ! aNucleon )
continue;
3219 G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
3220 G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
3221 if ( isProjectileNucleus ) {
3222 pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*
x);
3223 e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*
x);
3227 #ifdef debugPutOnMassShell
3228 G4cout <<
"i nY pY nY-AY AY " << i <<
" " << nucleonY <<
" " << projectileY <<
G4endl;
3231 if ( std::abs( nucleonY - nucleusY ) > 2 ||
3232 ( isProjectileNucleus && targetY > nucleonY ) ||
3233 ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
3246 const G4bool isProjectileNucleus,
3249 const G4int residualMassNumber,
3250 const G4int numberOfInvolvedNucleons,
3268 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3269 G4Nucleon* aNucleon = involvedNucleons[i];
3270 if ( ! aNucleon )
continue;
3272 residual3Momentum -= tmp.vect();
3276 G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w *
x );
3277 G4double e = w * x / 2.0 + mt2 / ( 2.0 * w *
x );
3279 if ( isProjectileNucleus ) pz *= -1.0;
3282 tmp.transform( boostFromCmsToLab );
3288 G4double residualMt2 =
sqr( residualMass ) +
sqr( residual3Momentum.x() )
3289 +
sqr( residual3Momentum.y() );
3291 #ifdef debugPutOnMassShell
3292 G4cout <<
"w residual3Momentum.z() " << w <<
" " << residual3Momentum.z() <<
G4endl;
3297 if ( residualMassNumber != 0 ) {
3298 residualPz = -w * residual3Momentum.z() / 2.0 +
3299 residualMt2 / ( 2.0 * w * residual3Momentum.z() );
3300 residualE = w * residual3Momentum.z() / 2.0 +
3301 residualMt2 / ( 2.0 * w * residual3Momentum.z() );
3303 if ( isProjectileNucleus ) residualPz *= -1.0;
3306 residual4Momentum.setPx( residual3Momentum.x() );
3307 residual4Momentum.setPy( residual3Momentum.y() );
3308 residual4Momentum.setPz( residualPz );
3309 residual4Momentum.setE( residualE );
3318 desc <<
" FTF (Fritiof) Model \n"
3319 <<
"The FTF model is based on the well-known FRITIOF \n"
3320 <<
"model (B. Andersson et al., Nucl. Phys. B281, 289 \n"
3321 <<
"(1987)). Its first program implementation was given\n"
3322 <<
"by B. Nilsson-Almquist and E. Stenlund (Comp. Phys.\n"
3323 <<
"Comm. 43, 387 (1987)). The Fritiof model assumes \n"
3324 <<
"that all hadron-hadron interactions are binary \n"
3325 <<
"reactions, h_1+h_2->h_1'+h_2' where h_1' and h_2' \n"
3326 <<
"are excited states of the hadrons with continuous \n"
3327 <<
"mass spectra. The excited hadrons are considered as\n"
3328 <<
"QCD-strings, and the corresponding LUND-string \n"
3329 <<
"fragmentation model is applied for a simulation of \n"
3330 <<
"their decays. \n"
3331 <<
" The Fritiof model assumes that in the course of \n"
3332 <<
"a hadron-nucleus interaction a string originated \n"
3333 <<
"from the projectile can interact with various intra\n"
3334 <<
"nuclear nucleons and becomes into highly excited \n"
3335 <<
"states. The probability of multiple interactions is\n"
3336 <<
"calculated in the Glauber approximation. A cascading\n"
3337 <<
"of secondary particles was neglected as a rule. Due\n"
3338 <<
"to these, the original Fritiof model fails to des- \n"
3339 <<
"cribe a nuclear destruction and slow particle spectra.\n"
3340 <<
" In order to overcome the difficulties we enlarge\n"
3341 <<
"the model by the reggeon theory inspired model of \n"
3342 <<
"nuclear desctruction (Kh. Abdel-Waged and V.V. Uzhi-\n"
3343 <<
"nsky, Phys. Atom. Nucl. 60, 828 (1997); Yad. Fiz. 60, 925\n"
3344 <<
"(1997)). Momenta of the nucleons ejected from a nuc-\n"
3345 <<
"leus in the reggeon cascading are sampled according\n"
3346 <<
"to a Fermi motion algorithm presented in (EMU-01 \n"
3347 <<
"Collaboration (M.I. Adamovich et al.) Zeit. fur Phys.\n"
3348 <<
"A358, 337 (1997)). \n"
3349 <<
" New features were also added to the Fritiof model\n"
3350 <<
"implemented in Geant4: a simulation of elastic had-\n"
3351 <<
"ron-nucleon scatterings, a simulation of binary \n"
3352 <<
"reactions like NN>NN* in hadron-nucleon interactions,\n"
3353 <<
"a separate simulation of single diffractive and non-\n"
3354 <<
" diffractive events. These allowed to describe after\n"
3355 <<
"model parameter tuning a wide set of experimental \n"
G4ElasticHNScattering * theElastic
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
void GetList(const G4ReactionProduct &thePrimary, G4FTFParameters *theParameters)
std::vector< G4ExcitedString * > G4ExcitedStringVector
G4Nucleon * GetTargetNucleon() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4Nucleon * GetProjectileNucleon() const
void SetParticleType(G4Proton *aProton)
G4double GetProbabilityOfElasticScatt()
G4double GetTotalMomentum() const
G4FTFModel(const G4String &modelName="FTF")
void SetMomentum(G4LorentzVector &aMomentum)
CLHEP::Hep3Vector G4ThreeVector
void operator()(G4VSplitableHadron *aH)
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
virtual G4bool StartLoop()=0
G4int GetSoftCollisionCount()
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4bool AdjustNucleons(G4VSplitableHadron *SelectedAntiBaryon, G4Nucleon *ProjectileNucleon, G4VSplitableHadron *SelectedTargetNucleon, G4Nucleon *TargetNucleon, G4bool Annihilation)
virtual void DoLorentzContraction(const G4LorentzVector &theBoost)=0
virtual G4int GetMassNumber()=0
G4Parton * GetLeftParton(void) const
const G4double w[NPOINTSGL]
void SetTimeOfCreation(G4double aTime)
G4int GetPDGEncoding() const
virtual const G4LorentzVector & Get4Momentum() const
G4ReactionProduct theProjectile
G4V3DNucleus * GetTargetNucleus() const
virtual const G4ThreeVector & GetPosition() const
G4int TargetResidualMassNumber
void SetDefinition(const G4ParticleDefinition *aDefinition)
virtual void SetProjectileNucleus(G4V3DNucleus *aNucleus)
void SetStatus(const G4int aStatus)
G4FTFParameters * theParameters
G4FTFAnnihilation * theAnnihilation
G4bool ExciteParticipants()
const G4String & GetParticleName() const
virtual const G4ParticleDefinition * GetDefinition() const
void SetThisPointer(G4VPartonStringModel *aPointer)
G4VSplitableHadron * GetSplitableHadron() const
G4double GetCofNuclearDestructionPr()
const G4ParticleDefinition * GetDefinition() const
G4Nucleon * TheInvolvedNucleonsOfTarget[250]
const G4ThreeVector & GetPosition() const
std::vector< G4VSplitableHadron * > theAdditionalString
G4InteractionContent & GetInteraction()
const G4ParticleDefinition * GetDefinition() const
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
G4double GetMaxNumberOfCollisions()
virtual void Init(G4int theZ, G4int theA)
static G4AntiProton * AntiProton()
G4FTFParticipants theParticipants
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double AbsoluteLevel)
G4V3DNucleus * GetProjectileNucleus() const
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
G4int ProjectileResidualMassNumber
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
G4DiffractiveExcitation * theExcitation
void SetTotalEnergy(const G4double en)
static const double twopi
G4double GetTimeOfCreation()
G4ExcitedStringVector * BuildStrings()
G4V3DNucleus * theProjectileNucleus
static G4Proton * Proton()
G4double GetMaxPt2ofNuclearDestruction()
static const double perCent
G4bool ComputeNucleusProperties(G4V3DNucleus *nucleus, G4LorentzVector &nucleusMomentum, G4LorentzVector &residualMomentum, G4double &sumMasses, G4double &residualExcitationEnergy, G4double &residualMass, G4int &residualMassNumber, G4int &residualCharge)
G4Nucleon * TheInvolvedNucleonsOfProjectile[250]
static G4Neutron * Neutron()
void SetBindingEnergy(G4double anEnergy)
void SetProbabilityOfElasticScatt(const G4double Xtotal, const G4double Xelastic)
void Init(const G4Nucleus &aNucleus, const G4DynamicParticle &aProjectile)
const G4LorentzVector & Get4Momentum() const
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
virtual void ModelDescription(std::ostream &) const
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
G4int ProjectileResidualCharge
G4double GetTotalEnergy() const
G4Parton * GetRightParton(void) const
virtual void InitProjectileNucleus(G4int theZ, G4int theA)
G4double GetPDGMass() const
virtual void DoLorentzBoost(const G4LorentzVector &theBoost)=0
void SetStatus(G4int aValue)
static G4ParticleTable * GetParticleTable()
G4double GetPt2ofNuclearDestruction()
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)
G4double GetR2ofNuclearDestruction()
G4int NumberOfInvolvedNucleonsOfProjectile
void StoreInvolvedNucleon()
const G4double x[NPOINTSGL]
G4double GetExcitationEnergyPerWoundedNucleon()
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)
G4bool GenerateDeltaIsobar(const G4double sqrtS, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4double &sumMasses)
G4VSplitableHadron * GetTarget() const
G4VSplitableHadron * GetProjectile() const
G4ThreeVector GetMomentum() const
const G4ThreeVector & GetPosition() const
virtual G4bool Annihilate(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4VSplitableHadron *&AdditionalString, G4FTFParameters *theParameters) const
G4bool FinalizeKinematics(const G4double w, const G4bool isProjectileNucleus, const G4LorentzRotation &boostFromCmsToLab, const G4double residualMass, const G4int residualMassNumber, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4LorentzVector &residual4Momentum)
G4int TargetResidualCharge
G4double GetCofNuclearDestruction()
G4double GetProbabilityOfAnnihilation()
G4LorentzVector TargetResidual4Momentum
virtual G4Nucleon * GetNextNucleon()=0
void Set4Momentum(const G4LorentzVector &a4Momentum)
void Hit(G4VSplitableHadron *aHit)
G4double GetPDGCharge() const
G4LorentzVector ProjectileResidual4Momentum
G4double GetBindingEnergy() const
G4double ProjectileResidualExcitationEnergy
G4ExcitedStringVector * GetStrings()
static G4AntiNeutron * AntiNeutron()
G4int GetBaryonNumber() const
G4int NumberOfInvolvedNucleonsOfTarget
G4double TargetResidualExcitationEnergy
CLHEP::HepLorentzVector G4LorentzVector
G4double GetDofNuclearDestruction()