77 NumberOfInvolvedNucleonsOfTarget = 0;
78 NumberOfInvolvedNucleonsOfProjectile= 0;
80 LowEnergyLimit = 5000.0*
MeV;
81 HighEnergyInter =
true;
84 ProjectileResidual4Momentum =
tmp;
85 ProjectileResidualMassNumber = 0;
86 ProjectileResidualCharge = 0;
87 ProjectileResidualExcitationEnergy = 0.0;
89 TargetResidual4Momentum =
tmp;
90 TargetResidualMassNumber = 0;
91 TargetResidualCharge = 0;
92 TargetResidualExcitationEnergy = 0.0;
109 if ( theParameters != 0 )
delete theParameters;
110 if ( theExcitation != 0 )
delete theExcitation;
111 if ( theElastic != 0 )
delete theElastic;
112 if ( theAnnihilation != 0 )
delete theAnnihilation;
115 if ( theAdditionalString.size() != 0 ) {
116 std::for_each( theAdditionalString.begin(), theAdditionalString.end(),
119 theAdditionalString.clear();
122 if ( NumberOfInvolvedNucleonsOfTarget != 0 ) {
123 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
125 if ( aNucleon )
delete aNucleon;
130 if ( NumberOfInvolvedNucleonsOfProjectile != 0 ) {
131 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
133 if ( aNucleon )
delete aNucleon;
143 theProjectile = aProjectile;
149 <<
"FTF init Proj Mass " << theProjectile.
GetMass()
153 <<
"FTF init Target A Z " << aNucleus.
GetA_asInt()
160 ProjectileResidualMassNumber = 0;
161 ProjectileResidualCharge = 0;
162 ProjectileResidualExcitationEnergy = 0.0;
163 ProjectileResidual4Momentum =
tmp;
165 TargetResidualMassNumber = aNucleus.
GetA_asInt();
167 TargetResidualExcitationEnergy = 0.0;
168 TargetResidual4Momentum =
tmp;
170 ->
GetIonMass( TargetResidualCharge, TargetResidualMassNumber );
172 TargetResidual4Momentum.
setE( TargetResidualMass );
179 ProjectileResidualExcitationEnergy = 0.0;
183 if ( PlabPerParticle < LowEnergyLimit ) {
184 HighEnergyInter =
false;
186 HighEnergyInter =
true;
197 if ( PlabPerParticle < LowEnergyLimit ) {
198 HighEnergyInter =
false;
200 HighEnergyInter =
true;
220 if ( PlabPerParticle < LowEnergyLimit ) {
221 HighEnergyInter =
false;
223 HighEnergyInter =
true;
229 ProjectileResidualExcitationEnergy = 0.0;
238 if ( theParameters != 0 )
delete theParameters;
242 if ( theAdditionalString.size() != 0 ) {
243 std::for_each( theAdditionalString.begin(), theAdditionalString.end(),
246 theAdditionalString.clear();
264 theParticipants.
GetList( theProjectile, theParameters );
265 StoreInvolvedNucleon();
269 if ( HighEnergyInter ) {
276 Success = PutOnMassShell();
279 G4cout <<
"FTF PutOnMassShell Success? " << Success <<
G4endl;
288 if ( Success ) Success = ExciteParticipants();
291 G4cout <<
"FTF ExciteParticipants Success? " << Success <<
G4endl;
297 G4cout <<
"FTF BuildStrings ";
300 theStrings = BuildStrings();
303 G4cout <<
"FTF BuildStrings " << theStrings <<
" OK" << G4endl
304 <<
"FTF GetResiduals of Nuclei " <<
G4endl;
309 if ( theParameters != 0 ) {
310 delete theParameters;
315 std::vector< G4VSplitableHadron* > primaries;
317 while ( theParticipants.
Next() ) {
320 if ( primaries.end() ==
321 std::find( primaries.begin(), primaries.end(), interaction.
GetProjectile() ) ) {
333 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
335 if ( aNucleon )
delete aNucleon;
337 NumberOfInvolvedNucleonsOfProjectile = 0;
340 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
342 if ( aNucleon )
delete aNucleon;
344 NumberOfInvolvedNucleonsOfTarget = 0;
347 G4cout <<
"End of FTF. Go to fragmentation" << G4endl
348 <<
"To continue - enter 1, to stop - ^C" <<
G4endl;
358 void G4FTFModel::StoreInvolvedNucleon() {
361 NumberOfInvolvedNucleonsOfTarget = 0;
369 TheInvolvedNucleonsOfTarget[NumberOfInvolvedNucleonsOfTarget] = aNucleon;
370 NumberOfInvolvedNucleonsOfTarget++;
375 G4cout <<
"G4FTFModel::StoreInvolvedNucleon -------------" <<
G4endl;
376 G4cout <<
"NumberOfInvolvedNucleonsOfTarget " << NumberOfInvolvedNucleonsOfTarget
384 NumberOfInvolvedNucleonsOfProjectile = 0;
390 while ( ( aProjectileNucleon = theProjectileNucleus->
GetNextNucleon() ) ) {
393 TheInvolvedNucleonsOfProjectile[NumberOfInvolvedNucleonsOfProjectile] = aProjectileNucleon;
394 NumberOfInvolvedNucleonsOfProjectile++;
399 G4cout <<
"NumberOfInvolvedNucleonsOfProjectile " << NumberOfInvolvedNucleonsOfProjectile
409 void G4FTFModel::ReggeonCascade() {
414 #ifdef debugReggeonCascade
415 G4cout <<
"G4FTFModel::ReggeonCascade -----------" << G4endl
416 <<
"theProjectile.GetTotalMomentum() " << theProjectile.
GetTotalMomentum() << G4endl
417 <<
"theProjectile.GetTotalEnergy() " << theProjectile.
GetTotalEnergy() << G4endl
418 <<
"ExcitationE/WN " << ExcitationE <<
G4endl;
421 G4int InitNINt = NumberOfInvolvedNucleonsOfTarget;
424 for (
G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
425 G4Nucleon* aTargetNucleon = TheInvolvedNucleonsOfTarget[ InvTN ];
438 if ( ! Neighbour->AreYouHit() ) {
439 G4double impact2 =
sqr( XofWoundedNucleon - Neighbour->GetPosition().x() ) +
440 sqr( YofWoundedNucleon - Neighbour->GetPosition().y() );
446 TheInvolvedNucleonsOfTarget[ NumberOfInvolvedNucleonsOfTarget ] = Neighbour;
447 NumberOfInvolvedNucleonsOfTarget++;
452 Neighbour->Hit( targetSplitable );
460 #ifdef debugReggeonCascade
461 G4cout <<
"Final NumberOfInvolvedNucleonsOfTarget "
462 << NumberOfInvolvedNucleonsOfTarget << G4endl <<
G4endl;
468 for (
G4int InvPN = 0; InvPN < NumberOfInvolvedNucleonsOfProjectile; InvPN++ ) {
469 G4Nucleon* aProjectileNucleon = TheInvolvedNucleonsOfProjectile[ InvPN ];
481 while ( ( Neighbour = theProjectileNucleus->
GetNextNucleon() ) ) {
482 if ( ! Neighbour->AreYouHit() ) {
483 G4double impact2=
sqr( XofWoundedNucleon - Neighbour->GetPosition().x() ) +
484 sqr( YofWoundedNucleon - Neighbour->GetPosition().y() );
490 TheInvolvedNucleonsOfProjectile[ NumberOfInvolvedNucleonsOfProjectile ] = Neighbour;
491 NumberOfInvolvedNucleonsOfProjectile++;
496 Neighbour->Hit( projectileSplitable );
504 #ifdef debugReggeonCascade
505 G4cout <<
"NumberOfInvolvedNucleonsOfProjectile "
506 << NumberOfInvolvedNucleonsOfProjectile << G4endl <<
G4endl;
513 G4bool G4FTFModel::PutOnMassShell() {
515 #ifdef debugPutOnMassShell
522 if ( Pprojectile.z() < 0.0 ) {
526 G4double M2projectile = Pprojectile.mag2();
534 G4double ExcitationEnergyPerWoundedNucleon =
537 #ifdef debugPutOnMassShell
538 G4cout <<
"ExcitationEnergyPerWoundedNucleon " << ExcitationEnergyPerWoundedNucleon <<
G4endl;
549 SumMasses += 20.0*
MeV;
550 TargetResidualExcitationEnergy += ExcitationEnergyPerWoundedNucleon;
551 TargetResidualMassNumber--;
558 #ifdef debugPutOnMassShell
559 G4cout <<
"Target residual: Charge, MassNumber " << TargetResidualCharge <<
" "
560 << TargetResidualMassNumber << G4endl <<
"Target Initial Momentum " << Ptarget
561 << G4endl <<
"Target Residual Momentum " << PtargetResidual <<
G4endl;
565 PtargetResidual.
setPz( 0.0 ); PtargetResidual.setE( 0.0 );
567 if ( TargetResidualMassNumber == 0 ) {
568 TargetResidualMass = 0.0;
569 TargetResidualExcitationEnergy = 0.0;
572 GetIonMass( TargetResidualCharge, TargetResidualMassNumber );
573 if ( TargetResidualMassNumber == 1 ) {
574 TargetResidualExcitationEnergy = 0.0;
577 SumMasses += std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.perp2() );
582 #ifdef debugPutOnMassShell
583 G4cout <<
"Psum " << Psum/
GeV <<
" GeV" << G4endl <<
"SqrtS " << SqrtS/
GeV <<
" GeV" << G4endl
584 <<
"SumMasses And TargetResidualMass " << SumMasses/
GeV <<
" "
585 << TargetResidualMass/
GeV <<
" GeV" <<
G4endl;
588 if ( SqrtS < SumMasses ) {
592 SumMasses -= std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.perp2() );
593 SumMasses += std::sqrt(
sqr( TargetResidualMass + TargetResidualExcitationEnergy )
594 + PtargetResidual.perp2() );
596 if ( SqrtS < SumMasses ) {
597 SumMasses -= std::sqrt(
sqr( TargetResidualMass + TargetResidualExcitationEnergy )
598 + PtargetResidual.perp2() );
599 SumMasses += std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.perp2() );
600 TargetResidualExcitationEnergy = 0.0;
603 TargetResidualMass +=TargetResidualExcitationEnergy;
605 #ifdef debugPutOnMassShell
606 G4cout <<
"TargetResidualMass SumMasses TargetResidualExcitationEnergy "
607 << TargetResidualMass/
GeV <<
" " << SumMasses/
GeV <<
" "
608 << TargetResidualExcitationEnergy <<
" MeV" <<
G4endl;
617 G4int MaxNumberOfDeltas =
G4int( (SqrtS - SumMasses)/(400.0*
MeV) );
618 G4int NumberOfDeltas( 0 );
624 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
627 if (
G4UniformRand() < ProbDeltaIsobar && NumberOfDeltas < MaxNumberOfDeltas ) {
635 G4int newPDGcode = PDGcode/10; newPDGcode = newPDGcode*10 + 4;
643 if ( SqrtS < SumMasses + MassDel - MassNuc ) {
646 ProbDeltaIsobar = 0.0;
648 SumMasses += (MassDel - MassNuc);
658 if ( Ptmp.pz() <= 0.0 ) {
664 Ptmp = toCms*Ptarget;
665 G4double YtargetNucleus = Ptmp.rapidity();
672 #ifdef debugPutOnMassShell
673 G4cout <<
"Y targetNucleus " << YtargetNucleus << G4endl
675 <<
" AveragePt2 " << AveragePt2 <<
G4endl;
682 G4int NumberOfTries( 0 );
684 G4bool OuterSuccess(
true );
694 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
698 AveragePt2 *= ScaleFactor;
704 G4bool InerSuccess =
true;
712 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
713 aNucleon = TheInvolvedNucleonsOfTarget[i];
723 G4double DeltaX = ( PtSum.x() - PtargetResidual.x() )/NumberOfInvolvedNucleonsOfTarget;
724 G4double DeltaY = ( PtSum.y() - PtargetResidual.y() )/NumberOfInvolvedNucleonsOfTarget;
726 if ( TargetResidualMassNumber == 0 ) {
727 DeltaXminus = (XminusSum - 1.0) / NumberOfInvolvedNucleonsOfTarget;
735 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
736 aNucleon = TheInvolvedNucleonsOfTarget[i];
739 if ( TargetResidualMassNumber == 0 ) {
740 if ( Xminus <= 0.0 || Xminus > 1.0 ) {
745 if ( Xminus <= 0.0 || Xminus > 1.0 || XminusSum <= 0.0 || XminusSum > 1.0 ) {
753 +
sqr( Px ) +
sqr( Py ) ) / Xminus;
758 if ( InerSuccess && TargetResidualMassNumber != 0 ) {
759 M2target += (
sqr( TargetResidualMass ) + PtargetResidual.perp2() ) / XminusSum;
762 #ifdef debugPutOnMassShell
763 G4cout <<
"InerSuccess " << InerSuccess << G4endl <<
"SqrtS, Mp+Mt, Mt " << SqrtS/
GeV
764 <<
" " << ( Mprojectile + std::sqrt( M2target ) )/
GeV <<
" "
765 << std::sqrt( M2target )/
GeV << G4endl
766 <<
"To continue - enter 1, to stop - ^C" <<
G4endl;
770 }
while ( ! InerSuccess );
772 }
while ( SqrtS < Mprojectile + std::sqrt( M2target ) );
775 - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
776 WminusTarget = ( S - M2projectile + M2target + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS;
777 WplusProjectile = SqrtS - M2target / WminusTarget;
778 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
779 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
780 G4double Yprojectile = 0.5 * std::log( (Eprojectile + Pzprojectile)/
781 (Eprojectile - Pzprojectile) );
783 #ifdef debugPutOnMassShell
784 G4cout <<
"DecayMomentum2 " << DecayMomentum2 << G4endl
785 <<
"WminusTarget WplusProjectile " << WminusTarget <<
" " << WplusProjectile
786 << G4endl <<
"Yprojectile " << Yprojectile <<
G4endl;
790 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
791 aNucleon = TheInvolvedNucleonsOfTarget[i];
796 G4double Pz = -WminusTarget*Xminus/2.0 + Mt2/(2.0*WminusTarget*Xminus);
797 G4double E = WminusTarget*Xminus/2.0 + Mt2/(2.0*WminusTarget*Xminus);
798 G4double YtargetNucleon = 0.5 * std::log( (E + Pz)/(E - Pz) );
800 #ifdef debugPutOnMassShell
801 G4cout <<
"i YtN Ypr YtN-YtA " << i <<
" " << YtargetNucleon <<
" " << YtargetNucleus
802 <<
" " << YtargetNucleon - YtargetNucleus <<
G4endl;
805 if ( std::abs( YtargetNucleon - YtargetNucleus ) > 2 || Yprojectile < YtargetNucleon ) {
806 OuterSuccess =
false;
811 }
while ( ! OuterSuccess );
813 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
814 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
815 Pprojectile.setPz( Pzprojectile ); Pprojectile.setE( Eprojectile );
817 #ifdef debugPutOnMassShell
818 G4cout <<
"Proj after in CMS " << Pprojectile <<
G4endl;
822 Pprojectile.transform( toLab );
827 theParticipants.
Next();
831 #ifdef debugPutOnMassShell
837 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
838 aNucleon = TheInvolvedNucleonsOfTarget[i];
840 TargetResidual3Momentum -= tmp.
vect();
844 G4double Pz = -WminusTarget*Xminus/2.0 + Mt2/(2.0*WminusTarget*Xminus);
845 G4double E = WminusTarget*Xminus/2.0 + Mt2/(2.0*WminusTarget*Xminus);
854 G4double Mt2Residual =
sqr( TargetResidualMass ) +
sqr( TargetResidual3Momentum.x() )
855 +
sqr( TargetResidual3Momentum.y() );
857 #ifdef debugPutOnMassShell
858 G4cout <<
"WminusTarget TargetResidual3Momentum.z() " << WminusTarget <<
" "
859 << TargetResidual3Momentum.z() <<
G4endl;
864 if ( TargetResidualMassNumber != 0 ) {
865 PzResidual = -WminusTarget*TargetResidual3Momentum.z()/2.0 +
866 Mt2Residual/(2.0*WminusTarget*TargetResidual3Momentum.z());
867 EResidual = WminusTarget*TargetResidual3Momentum.z()/2.0 +
868 Mt2Residual/(2.0*WminusTarget*TargetResidual3Momentum.z());
871 TargetResidual4Momentum.
setPx( TargetResidual3Momentum.x() );
872 TargetResidual4Momentum.
setPy( TargetResidual3Momentum.y() );
873 TargetResidual4Momentum.
setPz( PzResidual );
874 TargetResidual4Momentum.
setE( EResidual );
876 #ifdef debugPutOnMassShell
877 G4cout <<
"Target Residual4Momentum in CMS" << TargetResidual4Momentum <<
G4endl;
880 TargetResidual4Momentum.
transform( toLab );
882 #ifdef debugPutOnMassShell
883 G4cout <<
"Target Residual4Momentum in Lab " << TargetResidual4Momentum << G4endl
884 <<
"To continue enter - 1, to break - ^C" <<
G4endl;
894 #ifdef debugPutOnMassShell
895 G4cout <<
"PutOnMassShell for Nucleus_Nucleus " <<
G4endl;
900 if ( Pprojectile.z() < 0.0 ) {
906 #ifdef debugPutOnMassShell
907 G4cout <<
"ExcitationEnergyPerWoundedNucleon " << ExcitationEnergyPerWoundedNucleon <<
G4endl;
927 SumMasses += 20.0*
MeV;
928 ProjectileResidualExcitationEnergy += ExcitationEnergyPerWoundedNucleon;
929 ProjectileResidualMassNumber--;
936 #ifdef debugPutOnMassShell
937 G4cout <<
"Projectile residual: Charge, MassNumber " << ProjectileResidualCharge <<
" "
938 << ProjectileResidualMassNumber << G4endl <<
"Initial Momentum " << Pproj << G4endl
939 <<
"Residual Momentum " << PprojResidual <<
G4endl;
956 SumMasses += 20.0*
MeV;
957 TargetResidualExcitationEnergy += ExcitationEnergyPerWoundedNucleon;
958 TargetResidualMassNumber--;
965 #ifdef debugPutOnMassShell
966 G4cout <<
"Target residual: Charge, MassNumber " << TargetResidualCharge <<
" "
967 << TargetResidualMassNumber << G4endl <<
"Initial Momentum " << Ptarget << G4endl
968 <<
"Residual Momentum " << PtargetResidual <<
G4endl;
973 PprojResidual.
setPz( 0.0 ); PprojResidual.setE( 0.0 );
974 PtargetResidual.setPz( 0.0 ); PtargetResidual.setE( 0.0 );
977 if ( ProjectileResidualMassNumber == 0 ) {
978 PrResidualMass = 0.0;
979 ProjectileResidualExcitationEnergy = 0.0;
982 ->
GetIonMass( ProjectileResidualCharge, ProjectileResidualMassNumber );
983 if ( ProjectileResidualMassNumber == 1 ) {
984 ProjectileResidualExcitationEnergy = 0.0;
988 SumMasses += std::sqrt(
sqr( PrResidualMass ) + PprojResidual.vect().perp2() );
991 if ( TargetResidualMassNumber == 0 ) {
992 TargetResidualMass = 0.0;
993 TargetResidualExcitationEnergy = 0.0;
996 ->
GetIonMass( TargetResidualCharge, TargetResidualMassNumber );
997 if ( TargetResidualMassNumber == 1 ) {
998 TargetResidualExcitationEnergy = 0.0;
1002 SumMasses += std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.perp2() );
1007 #ifdef debugPutOnMassShell
1008 G4cout <<
"Psum " << Psum/
GeV <<
" GeV" << G4endl <<
"SqrtS " << SqrtS/
GeV <<
" GeV" << G4endl
1009 <<
"SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/
GeV <<
" "
1010 << PrResidualMass/
GeV <<
" " << TargetResidualMass/
GeV <<
" GeV" <<
G4endl;
1013 if ( SqrtS < SumMasses ) {
1017 SumMasses -= std::sqrt(
sqr( PrResidualMass ) + PprojResidual.perp2() );
1018 SumMasses += std::sqrt(
sqr( PrResidualMass + ProjectileResidualExcitationEnergy ) +
1019 PprojResidual.perp2() );
1020 SumMasses -= std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.perp2() );
1021 SumMasses += std::sqrt(
sqr( TargetResidualMass + TargetResidualExcitationEnergy ) +
1022 PtargetResidual.perp2() );
1024 if ( SqrtS < SumMasses ) {
1025 SumMasses -= std::sqrt(
sqr( PrResidualMass + ProjectileResidualExcitationEnergy ) +
1026 PprojResidual.perp2() );
1027 SumMasses += std::sqrt(
sqr( PrResidualMass ) + PprojResidual.perp2() );
1028 ProjectileResidualExcitationEnergy = 0.0;
1029 SumMasses -= std::sqrt(
sqr( TargetResidualMass + TargetResidualExcitationEnergy ) +
1030 PtargetResidual.perp2() );
1031 SumMasses += std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.perp2() );
1032 TargetResidualExcitationEnergy = 0.0;
1035 PrResidualMass += ProjectileResidualExcitationEnergy;
1036 TargetResidualMass += TargetResidualExcitationEnergy;
1038 #ifdef debugPutOnMassShell
1039 G4cout <<
"PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/
GeV <<
" "
1040 << ProjectileResidualExcitationEnergy <<
" MeV" << G4endl
1041 <<
"TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/
GeV <<
" "
1042 << TargetResidualExcitationEnergy <<
" MeV" << G4endl
1043 <<
"Sum masses " << SumMasses/
GeV <<
G4endl;
1048 G4int MaxNumberOfDeltas =
G4int( (SqrtS - SumMasses)/(400.0*
MeV) );
1049 G4int NumberOfDeltas( 0 );
1055 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
1056 if (
G4UniformRand() < ProbDeltaIsobar && NumberOfDeltas < MaxNumberOfDeltas ) {
1064 G4int newPDGcode = PDGcode/10; newPDGcode = newPDGcode*10 + 4;
1070 if ( SqrtS < SumMasses + MassDel - MassNuc ) {
1072 ProbDeltaIsobar = 0.0;
1074 SumMasses += (MassDel - MassNuc);
1083 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
1084 if (
G4UniformRand() < ProbDeltaIsobar && NumberOfDeltas < MaxNumberOfDeltas ) {
1092 G4int newPDGcode = PDGcode/10; newPDGcode=newPDGcode*10 + 4;
1097 if ( SqrtS < SumMasses + MassDel - MassNuc ) {
1099 ProbDeltaIsobar = 0.0;
1101 SumMasses += (MassDel - MassNuc);
1109 if ( Ptmp.pz() <= 0.0 ) {
1117 G4double YprojectileNucleus = Ptmp.rapidity();
1118 Ptmp = toCms*Ptarget;
1119 G4double YtargetNucleus = Ptmp.rapidity();
1127 #ifdef debugPutOnMassShell
1128 G4cout <<
"Y projectileNucleus " << YprojectileNucleus << G4endl <<
"Y targetNucleus "
1130 <<
" DcorP DcorT " << DcorP <<
" " << DcorT <<
" AveragePt2 " << AveragePt2 <<
G4endl;
1137 G4int NumberOfTries( 0 );
1139 G4bool OuterSuccess(
true );
1143 OuterSuccess =
true;
1149 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1152 DcorP *= ScaleFactor;
1153 DcorT *= ScaleFactor;
1154 AveragePt2 *= ScaleFactor;
1161 G4bool InerSuccess =
true;
1170 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
1171 aNucleon = TheInvolvedNucleonsOfProjectile[i];
1172 G4ThreeVector tmpPt = GaussianPt( AveragePt2, maxPtSquare );
1181 G4double DeltaX = (PtSum.x() - PprojResidual.x()) / NumberOfInvolvedNucleonsOfProjectile;
1182 G4double DeltaY = (PtSum.y() - PprojResidual.y()) / NumberOfInvolvedNucleonsOfProjectile;
1184 if ( ProjectileResidualMassNumber == 0 ) {
1185 DeltaXplus = (XplusSum - 1.0) / NumberOfInvolvedNucleonsOfProjectile;
1192 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
1193 aNucleon = TheInvolvedNucleonsOfProjectile[i];
1196 if ( ProjectileResidualMassNumber == 0 ) {
1197 if ( Xplus <= 0.0 || Xplus > 1.0 ) {
1198 InerSuccess =
false;
1202 if ( Xplus <= 0.0 || Xplus > 1.0 || XplusSum <= 0.0 || XplusSum > 1.0 ) {
1203 InerSuccess =
false;
1210 sqr( Px ) +
sqr( Py ) ) / Xplus;
1215 if ( InerSuccess && ProjectileResidualMassNumber != 0 ) {
1216 M2proj += (
sqr( PrResidualMass ) + PprojResidual.perp2() ) / XplusSum;
1219 }
while ( ! InerSuccess );
1233 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
1234 aNucleon = TheInvolvedNucleonsOfTarget[i];
1235 G4ThreeVector tmpPt = GaussianPt( AveragePt2, maxPtSquare );
1239 XminusSum += Xminus;
1244 G4double DeltaX = (PtSum.x() - PtargetResidual.x()) / NumberOfInvolvedNucleonsOfTarget;
1245 G4double DeltaY = (PtSum.y() - PtargetResidual.y()) / NumberOfInvolvedNucleonsOfTarget;
1247 if ( TargetResidualMassNumber == 0 ) {
1248 DeltaXminus = (XminusSum - 1.0) / NumberOfInvolvedNucleonsOfTarget;
1256 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
1257 aNucleon = TheInvolvedNucleonsOfTarget[i];
1259 XminusSum -= Xminus;
1260 if ( TargetResidualMassNumber == 0 ) {
1261 if ( Xminus <= 0.0 || Xminus > 1.0 ) {
1262 InerSuccess =
false;
1266 if ( Xminus <= 0.0 || Xminus > 1.0 || XminusSum <= 0.0 || XminusSum > 1.0 ) {
1267 InerSuccess =
false;
1274 sqr( Px ) +
sqr( Py ) ) / Xminus;
1279 if ( InerSuccess && TargetResidualMassNumber != 0 ) {
1280 M2target += (
sqr( TargetResidualMass ) + PtargetResidual.perp2() ) / XminusSum;
1283 }
while ( ! InerSuccess );
1285 #ifdef debugPutOnMassShell
1286 G4cout <<
"SqrtS, Mp+Mt, Mp, Mt " << SqrtS/
GeV <<
" "
1287 << ( std::sqrt( M2proj ) + std::sqrt( M2target) )/
GeV <<
" "
1288 << std::sqrt( M2proj )/
GeV <<
" " << std::sqrt( M2target )/
GeV << G4endl
1289 <<
"To continue - enter 1, to stop - ^C" <<
G4endl;
1293 }
while ( SqrtS < std::sqrt( M2proj ) + std::sqrt( M2target ) );
1296 - 2.0*S*M2proj - 2.0*S*M2target - 2.0*M2proj*M2target;
1297 WminusTarget = ( S - M2proj + M2target + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS;
1298 WplusProjectile = SqrtS - M2target/WminusTarget;
1299 G4double Pzprojectile = WplusProjectile/2.0 - M2proj/2.0/WplusProjectile;
1300 G4double Eprojectile = WplusProjectile/2.0 + M2proj/2.0/WplusProjectile;
1301 G4double Yprojectile = 0.5 * std::log( (Eprojectile + Pzprojectile)/
1302 (Eprojectile - Pzprojectile) );
1303 G4double Pztarget = -WminusTarget/2.0 + M2target/2.0/WminusTarget;
1304 G4double Etarget = WminusTarget/2.0 + M2target/2.0/WminusTarget;
1305 G4double Ytarget = 0.5 * std::log( (Etarget + Pztarget)/(Etarget - Pztarget) );
1307 #ifdef debugPutOnMassShell
1308 G4cout <<
"DecayMomentum2 " << DecayMomentum2 << G4endl <<
"WminusTarget WplusProjectile "
1309 << WminusTarget <<
" " << WplusProjectile << G4endl
1310 <<
"Yprojectile " << Yprojectile <<
G4endl;
1315 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
1316 aNucleon = TheInvolvedNucleonsOfProjectile[i];
1321 G4double Pz = WplusProjectile*Xplus/2.0 - Mt2/(2.0*WplusProjectile*Xplus);
1322 G4double E = WplusProjectile*Xplus/2.0 + Mt2/(2.0*WplusProjectile*Xplus);
1323 G4double YprojNucleon = 0.5 * std::log( (E + Pz)/(E - Pz) );
1325 #ifdef debugPutOnMassShell
1326 G4cout <<
"i YpN Ypr YpN-YtA Ypr0 " << i <<
" " << YprojNucleon <<
" " << Yprojectile
1327 <<
" " << YprojNucleon - Yprojectile <<
" " << YprojectileNucleus <<
G4endl;
1330 if ( std::abs( YprojNucleon - YprojectileNucleus ) > 2 || Ytarget > YprojNucleon ) {
1331 OuterSuccess =
false;
1337 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
1338 aNucleon = TheInvolvedNucleonsOfTarget[i];
1343 G4double Pz = -WminusTarget*Xminus/2.0 + Mt2/(2.0*WminusTarget*Xminus);
1344 G4double E = WminusTarget*Xminus/2.0 + Mt2/(2.0*WminusTarget*Xminus);
1345 G4double YtargetNucleon = 0.5 * std::log( (E + Pz)/(E - Pz) );
1347 #ifdef debugPutOnMassShell
1348 G4cout <<
"i YtN Ypr YtN-YtA " << i <<
" " << YtargetNucleon <<
" " << Yprojectile
1349 <<
" " << YtargetNucleon - YtargetNucleus <<
G4endl;
1352 if ( std::abs( YtargetNucleon - YtargetNucleus ) > 2 || Yprojectile < YtargetNucleon ) {
1353 OuterSuccess =
false;
1358 }
while ( ! OuterSuccess );
1362 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
1363 aNucleon = TheInvolvedNucleonsOfProjectile[i];
1365 ProjectileResidual3Momentum -= tmp.
vect();
1369 G4double Pz = WplusProjectile*Xplus/2.0 - Mt2/(2.0*WplusProjectile*Xplus);
1370 G4double E = WplusProjectile*Xplus/2.0 + Mt2/(2.0*WplusProjectile*Xplus);
1379 G4double Mt2PrResidual =
sqr( PrResidualMass ) +
sqr( ProjectileResidual3Momentum.x() ) +
1380 sqr( ProjectileResidual3Momentum.y() );
1382 #ifdef debugPutOnMassShell
1383 G4cout <<
"WplusProjectile ProjectileResidual3Momentum.z() " << WplusProjectile <<
" "
1384 << ProjectileResidual3Momentum.z() <<
G4endl;
1389 if ( ProjectileResidualMassNumber != 0 ) {
1390 PzPrResidual = WplusProjectile*ProjectileResidual3Momentum.z()/2.0 -
1391 Mt2PrResidual/( 2.0*WplusProjectile*ProjectileResidual3Momentum.z() );
1392 EPrResidual = WplusProjectile*ProjectileResidual3Momentum.z()/2.0 +
1393 Mt2PrResidual/( 2.0*WplusProjectile*ProjectileResidual3Momentum.z() );
1395 ProjectileResidual4Momentum.
setPx( ProjectileResidual3Momentum.x());
1396 ProjectileResidual4Momentum.
setPy( ProjectileResidual3Momentum.y());
1397 ProjectileResidual4Momentum.
setPz( PzPrResidual );
1398 ProjectileResidual4Momentum.
setE( EPrResidual );
1400 #ifdef debugPutOnMassShell
1401 G4cout <<
"Projectile Residual4Momentum in CMS" << ProjectileResidual4Momentum <<
G4endl;
1404 ProjectileResidual4Momentum.
transform( toLab );
1406 #ifdef debugPutOnMassShell
1407 G4cout <<
"Projectile Residual4Momentum in Lab " << ProjectileResidual4Momentum <<
G4endl;
1411 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
1412 aNucleon = TheInvolvedNucleonsOfTarget[i];
1414 TargetResidual3Momentum -= tmp.
vect();
1418 G4double Pz = -WminusTarget*Xminus/2.0 + Mt2/(2.0*WminusTarget*Xminus);
1419 G4double E = WminusTarget*Xminus/2.0 + Mt2/(2.0*WminusTarget*Xminus);
1428 G4double Mt2TrResidual =
sqr( TargetResidualMass ) +
sqr( TargetResidual3Momentum.x() ) +
1429 sqr( TargetResidual3Momentum.y() );
1431 #ifdef debugPutOnMassShell
1432 G4cout <<
"WminusTarget TargetResidual3Momentum.z() " << WminusTarget
1433 <<
" " << TargetResidual3Momentum.z() <<
G4endl;
1438 if ( TargetResidualMassNumber != 0 ) {
1439 PzTrResidual = -WminusTarget*TargetResidual3Momentum.z()/2.0 +
1440 Mt2TrResidual/( 2.0*WminusTarget*TargetResidual3Momentum.z() );
1441 ETrResidual = WminusTarget*TargetResidual3Momentum.z()/2.0 +
1442 Mt2TrResidual/( 2.0*WminusTarget*TargetResidual3Momentum.z() );
1445 TargetResidual4Momentum.
setPx( TargetResidual3Momentum.x() );
1446 TargetResidual4Momentum.
setPy( TargetResidual3Momentum.y() );
1447 TargetResidual4Momentum.
setPz( PzTrResidual );
1448 TargetResidual4Momentum.
setE( ETrResidual );
1450 #ifdef debugPutOnMassShell
1451 G4cout <<
"Target Residual4Momentum in CMS" << TargetResidual4Momentum <<
G4endl;
1454 TargetResidual4Momentum.
transform( toLab );
1456 #ifdef debugPutOnMassShell
1457 G4cout <<
"Target Residual4Momentum in Lab " << TargetResidual4Momentum << G4endl
1458 <<
"To continue enter - 1, to break - ^C" <<
G4endl;
1468 G4bool G4FTFModel::ExciteParticipants() {
1470 #ifdef debugBuildString
1471 G4cout <<
"G4FTFModel::ExciteParticipants() " <<
G4endl;
1474 G4bool Successfull(
true );
1476 if ( MaxNumOfInelCollisions > 0 ) {
1478 if (
G4UniformRand() < ProbMaxNumber ) MaxNumOfInelCollisions++;
1481 MaxNumOfInelCollisions = 1;
1484 #ifdef debugBuildString
1485 G4cout <<
"MaxNumOfInelCollisions MaxNumOfInelCollisions " << MaxNumOfInelCollisions <<
G4endl;
1488 G4int CurrentInteraction( 0 );
1491 while ( theParticipants.
Next() ) {
1493 CurrentInteraction++;
1500 #ifdef debugBuildString
1501 G4cout << G4endl <<
"Interaction # Status " << CurrentInteraction <<
" "
1502 << collision.
GetStatus() << G4endl <<
"Pr* Tr* " << projectile <<
" "
1503 << target << G4endl <<
"projectile->GetStatus target->GetStatus "
1513 #ifdef debugBuildString
1517 if ( ! HighEnergyInter ) {
1518 G4bool Annihilation =
false;
1519 G4bool Result = AdjustNucleons( projectile, ProjectileNucleon, target,
1520 TargetNucleon, Annihilation );
1521 if ( ! Result )
continue;
1523 Successfull = theElastic->
ElasticScattering( projectile, target, theParameters )
1528 #ifdef debugBuildString
1529 G4cout <<
"Inelastic interaction" << G4endl
1530 <<
"MaxNumOfInelCollisions " << MaxNumOfInelCollisions <<
G4endl;
1533 if ( ! HighEnergyInter ) {
1534 G4bool Annihilation =
false;
1535 G4bool Result = AdjustNucleons( projectile, ProjectileNucleon, target,
1536 TargetNucleon, Annihilation );
1537 if ( ! Result )
continue;
1547 if (theExcitation->
ExciteParticipants( projectile, target, theParameters, theElastic )){
1549 #ifdef debugBuildString
1559 #ifdef debugBuildString
1560 G4cout <<
"FTF excitation Non Successfull -> Elastic scattering "
1561 << Successfull <<
G4endl;
1564 Successfull = theElastic->
ElasticScattering( projectile, target, theParameters )
1569 #ifdef debugBuildString
1570 G4cout <<
"Elastic scat. at rejection inelastic scattering" <<
G4endl;
1579 Successfull = theElastic->
ElasticScattering( projectile, target, theParameters )
1584 #ifdef debugBuildString
1589 while ( theParticipants.
Next() ) {
1593 if ( projectile == NextProjectileNucleon || target == NextTargetNucleon ) {
1600 for (
G4int I = 0; I < CurrentInteraction; I++ ) theParticipants.
Next();
1603 if ( ! HighEnergyInter ) {
1604 G4bool Annihilation =
true;
1605 G4bool Result = AdjustNucleons( projectile, ProjectileNucleon, target,
1606 TargetNucleon, Annihilation );
1607 if ( ! Result )
continue;
1610 if ( theAnnihilation->
Annihilate( projectile, target, AdditionalString, theParameters ) ){
1611 Successfull = Successfull ||
true;
1613 #ifdef debugBuildString
1614 G4cout <<
"Annihilation successfull. " <<
"*AdditionalString "
1615 << AdditionalString <<
G4endl;
1620 if ( AdditionalString != 0 ) theAdditionalString.push_back( AdditionalString );
1625 #ifdef debugBuildString
1626 G4cout <<
"----------------------------- Final properties " << G4endl
1627 <<
"projectile->GetStatus target->GetStatus " << projectile->
GetStatus()
1628 <<
" " << target->
GetStatus() << G4endl <<
"projectile->GetSoftC target->GetSoftC "
1630 << G4endl <<
"ExciteParticipants() Successfull? " << Successfull <<
G4endl;
1648 G4cout <<
"AdjustNucleons ---------------------------------------" << G4endl
1650 <<
"Proj 4mom " << SelectedAntiBaryon->
Get4Momentum() << G4endl
1651 <<
"Targ 4mom " << SelectedTargetNucleon->
Get4Momentum() << G4endl
1652 <<
"Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
1653 << ProjectileResidualMassNumber <<
" " << ProjectileResidualCharge <<
" "
1654 << ProjectileResidualExcitationEnergy << G4endl
1655 <<
"Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
1656 << TargetResidualMassNumber <<
" " << TargetResidualCharge <<
" "
1657 << TargetResidualExcitationEnergy << G4endl
1670 G4double DcorP( 0.0 ), DcorT( 0.0 );
1671 if ( ProjectileResidualMassNumber != 0 ) DcorP = Dcor /
G4double(ProjectileResidualMassNumber);
1672 if ( TargetResidualMassNumber != 0 ) DcorT = Dcor /
G4double(TargetResidualMassNumber);
1676 G4double ExcitationEnergyPerWoundedNucleon =
1690 G4cout <<
"case 1, hA prcol=0 trcol=0, AA prcol#0 trcol=0" <<
G4endl;
1693 if ( TargetResidualMassNumber < 1 ) {
1701 if ( TargetResidualMassNumber == 1 ) {
1702 TargetResidualMassNumber = 0;
1703 TargetResidualCharge = 0;
1704 TargetResidualExcitationEnergy = 0.0;
1705 SelectedTargetNucleon->
Set4Momentum( TargetResidual4Momentum );
1713 G4cout <<
"Targ res Init " << TargetResidual4Momentum <<
G4endl;
1730 G4int TResidualMassNumber = TargetResidualMassNumber - 1;
1731 G4int TResidualCharge = TargetResidualCharge -
1733 G4double TResidualExcitationEnergy = TargetResidualExcitationEnergy +
1734 ExcitationEnergyPerWoundedNucleon;
1735 if ( TResidualMassNumber <= 1 ) {
1736 TResidualExcitationEnergy = 0.0;
1740 if ( TResidualMassNumber != 0 ) {
1742 ->
GetIonMass( TResidualCharge, TResidualMassNumber );
1754 if ( ! Annihilation ) {
1755 if ( SqrtS < SumMasses ) {
1758 if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1761 G4cout <<
"TResidualExcitationEnergy " << TResidualExcitationEnergy <<
G4endl;
1764 TResidualExcitationEnergy = SqrtS - SumMasses;
1767 G4cout <<
"TResidualExcitationEnergy " << TResidualExcitationEnergy <<
G4endl;
1775 if ( Annihilation ) {
1778 G4cout <<
"SqrtS < SumMasses - TNucleonMass " << SqrtS <<
" "
1779 << SumMasses - TNucleonMass <<
G4endl;
1782 if ( SqrtS < SumMasses - TNucleonMass ) {
1787 G4cout <<
"SqrtS < SumMasses " << SqrtS <<
" " << SumMasses <<
G4endl;
1790 if ( SqrtS < SumMasses ) {
1791 TNucleonMass = SqrtS - (SumMasses - TNucleonMass) - TResidualExcitationEnergy;
1797 SumMasses = SqrtS - TResidualExcitationEnergy;
1803 G4cout <<
"SqrtS < SumMasses " << SqrtS <<
" " << SumMasses <<
G4endl;
1806 if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1807 TResidualExcitationEnergy = SqrtS - SumMasses;
1826 Pprojectile = Ptmp; Pprojectile.
transform( toLab );
1830 Ptmp.
setE( TNucleonMass );
1836 Ptarget = Ptmp; Ptarget.
transform( toLab );
1840 TargetResidualMassNumber = TResidualMassNumber;
1841 TargetResidualCharge = TResidualCharge;
1842 TargetResidualExcitationEnergy = TResidualExcitationEnergy;
1844 Ptmp.
setE( TResidualMass + TargetResidualExcitationEnergy );
1851 TargetResidual4Momentum = Ptmp;
1854 G4cout << Pprojectile << G4endl << Ptarget << G4endl << TargetResidual4Momentum <<
G4endl;
1860 G4double Mprojectile = Pprojectile.mag();
1861 G4double M2projectile = Pprojectile.mag2();
1867 TResidualMass += TResidualExcitationEnergy;
1876 G4int NumberOfTries( 0 );
1878 G4bool OuterSuccess(
true );
1881 OuterSuccess =
true;
1887 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1890 DcorT *= ScaleFactor;
1891 AveragePt2 *= ScaleFactor;
1901 G4bool InerSuccess =
true;
1902 if ( TargetResidualMassNumber > 1 ) {
1906 PtNucleon = GaussianPt( AveragePt2, maxPtSquare );
1907 PtResidual = -PtNucleon;
1909 G4double Mtarget = std::sqrt(
sqr( TNucleonMass ) + PtNucleon.mag2() ) +
1910 std::sqrt(
sqr( TResidualMass ) + PtResidual.mag2() );
1911 if ( SqrtS < Mprojectile + Mtarget ) {
1912 InerSuccess =
false;
1917 G4double Xcenter = std::sqrt(
sqr( TNucleonMass ) + PtNucleon.mag2() ) / Mtarget;
1918 XminusNucleon = Xcenter + tmpX.
x();
1919 if ( XminusNucleon <= 0.0 || XminusNucleon >= 1.0 ) {
1920 InerSuccess =
false;
1924 XminusResidual = 1.0 - XminusNucleon;
1925 }
while ( ! InerSuccess );
1927 XminusNucleon = 1.0;
1928 XminusResidual = 1.0;
1932 M2target = (
sqr( TNucleonMass ) + PtNucleon.mag2() ) / XminusNucleon +
1933 (
sqr( TResidualMass ) + PtResidual.mag2() ) / XminusResidual;
1935 }
while ( SqrtS < Mprojectile + std::sqrt( M2target) );
1938 - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
1940 WminusTarget = ( S - M2projectile + M2target + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS;
1941 WplusProjectile = SqrtS - M2target / WminusTarget;
1943 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
1944 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
1945 G4double Yprojectile = 0.5 * std::log( (Eprojectile + Pzprojectile) /
1946 (Eprojectile - Pzprojectile) );
1949 G4cout <<
"DecayMomentum2 " << DecayMomentum2 << G4endl
1950 <<
"WminusTarget WplusProjectile " << WminusTarget <<
" " << WplusProjectile
1951 << G4endl <<
"Yprojectile " << Yprojectile <<
G4endl;
1954 G4double Mt2 =
sqr( TNucleonMass ) + PtNucleon.mag2();
1955 G4double Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1956 G4double E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1957 G4double YtargetNucleon = 0.5 * std::log( (E + Pz)/(E - Pz) );
1960 G4cout <<
"YtN Ytr YtN-Ytr " <<
" " << YtargetNucleon <<
" " << YtargetNucleus <<
" "
1961 << YtargetNucleon - YtargetNucleus << G4endl
1962 <<
"YtN Ypr YtN-Ypr " <<
" " << YtargetNucleon <<
" " << Yprojectile
1963 <<
" " << YtargetNucleon - Yprojectile <<
G4endl;
1966 if ( std::abs( YtargetNucleon - YtargetNucleus ) > 2 || Yprojectile < YtargetNucleon ) {
1967 OuterSuccess =
false;
1971 }
while ( ! OuterSuccess );
1973 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
1974 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
1975 Pprojectile.setPz( Pzprojectile ); Pprojectile.setE( Eprojectile );
1978 G4cout <<
"Proj after in CMS " << Pprojectile <<
G4endl;
1981 Pprojectile.transform( toLab );
1989 G4double Mt2 =
sqr( TNucleonMass ) + PtNucleon.mag2();
1990 G4double Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1991 G4double E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1993 Ptarget.setPx( PtNucleon.x() ); Ptarget.setPy( PtNucleon.y() );
1994 Ptarget.setPz( Pz ); Ptarget.setE( E );
1995 Ptarget.transform( toLab );
2003 TargetResidualMassNumber = TResidualMassNumber;
2004 TargetResidualCharge = TResidualCharge;
2005 TargetResidualExcitationEnergy = TResidualExcitationEnergy;
2008 G4cout <<
"TargetResidualMassNumber TargetResidualCharge TargetResidualExcitationEnergy "
2009 << TargetResidualMassNumber <<
" " << TargetResidualCharge <<
" "
2010 << TargetResidualExcitationEnergy <<
G4endl;
2013 if ( TargetResidualMassNumber != 0 ) {
2014 Mt2 =
sqr( TResidualMass ) + PtResidual.mag2();
2015 Pz = -WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
2016 E = WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
2018 TargetResidual4Momentum.setPx( PtResidual.x() );
2019 TargetResidual4Momentum.setPy( PtResidual.y() );
2020 TargetResidual4Momentum.setPz( Pz );
2021 TargetResidual4Momentum.setE( E );
2024 G4cout <<
"New Residu " << TargetResidual4Momentum <<
" CMS" <<
G4endl;
2027 TargetResidual4Momentum.transform( toLab );
2030 G4cout <<
"New Residu " << TargetResidual4Momentum <<
" Lab" <<
G4endl;
2046 if ( ProjectileResidualMassNumber < 1 )
return false;
2048 if ( ProjectileResidual4Momentum.
rapidity() <=
2053 if ( ProjectileResidualMassNumber == 1 ) {
2054 ProjectileResidualMassNumber = 0;
2055 ProjectileResidualCharge = 0;
2056 ProjectileResidualExcitationEnergy = 0.0;
2057 SelectedAntiBaryon->
Set4Momentum( ProjectileResidual4Momentum );
2077 G4int TResidualMassNumber = ProjectileResidualMassNumber - 1;
2078 G4int TResidualCharge = ProjectileResidualCharge
2080 G4double TResidualExcitationEnergy = ProjectileResidualExcitationEnergy +
2081 ExcitationEnergyPerWoundedNucleon;
2082 if ( TResidualMassNumber <= 1 ) {
2083 TResidualExcitationEnergy = 0.0;
2087 if ( TResidualMassNumber != 0 ) {
2089 ->
GetIonMass( TResidualCharge , TResidualMassNumber );
2095 TNucleonMass + TResidualMass;
2098 G4cout <<
"SelectedTN.mag() PNMass + PResidualMass "
2100 << TNucleonMass <<
" " << TResidualMass <<
G4endl;
2105 if ( ! Annihilation ) {
2106 if ( SqrtS < SumMasses ) {
2109 if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
2110 TResidualExcitationEnergy = SqrtS - SumMasses;
2116 if ( Annihilation ) {
2117 if ( SqrtS < SumMasses - TNucleonMass ) {
2120 if ( SqrtS < SumMasses ) {
2121 TNucleonMass = SqrtS - (SumMasses - TNucleonMass);
2123 TResidualExcitationEnergy = 0.0;
2127 if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
2128 TResidualExcitationEnergy = SqrtS - SumMasses;
2142 Ptarget = Ptmp; Ptarget.
transform( toLab );
2146 Ptmp.
setE( TNucleonMass );
2147 Pprojectile = Ptmp; Pprojectile.
transform( toLab );
2151 ProjectileResidualMassNumber = TResidualMassNumber;
2152 ProjectileResidualCharge = TResidualCharge;
2153 ProjectileResidualExcitationEnergy = TResidualExcitationEnergy;
2155 Ptmp.
setE( TResidualMass + ProjectileResidualExcitationEnergy );
2157 ProjectileResidual4Momentum = Ptmp;
2163 G4double M2target = Ptarget.mag2();
2165 G4LorentzVector TResidual4Momentum = toCms * ProjectileResidual4Momentum;
2168 TResidualMass += TResidualExcitationEnergy;
2177 G4int NumberOfTries( 0 );
2179 G4bool OuterSuccess(
true );
2183 OuterSuccess =
true;
2189 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
2192 DcorP *= ScaleFactor;
2193 AveragePt2 *= ScaleFactor;
2197 G4cout <<
"ProjectileResidualMassNumber " << ProjectileResidualMassNumber <<
G4endl;
2200 if ( ProjectileResidualMassNumber > 1 ) {
2201 PtNucleon = GaussianPt( AveragePt2, maxPtSquare );
2205 PtResidual = -PtNucleon;
2207 G4double Mprojectile = std::sqrt(
sqr( TNucleonMass ) + PtNucleon.mag2() ) +
2208 std::sqrt(
sqr( TResidualMass ) + PtResidual.mag2() );
2211 G4cout <<
"SqrtS < Mtarget + Mprojectile " << SqrtS <<
" " << Mtarget
2212 <<
" " << Mprojectile <<
" " << Mtarget + Mprojectile <<
G4endl;
2215 M2projectile =
sqr( Mprojectile );
2216 if ( SqrtS < Mtarget + Mprojectile ) {
2217 OuterSuccess =
false;
2221 G4double Xcenter = std::sqrt(
sqr( TNucleonMass ) + PtNucleon.mag2() ) / Mprojectile;
2223 G4bool InerSuccess =
true;
2224 if ( ProjectileResidualMassNumber > 1 ) {
2228 XplusNucleon = Xcenter + tmpX.
x();
2229 if ( XplusNucleon <= 0.0 || XplusNucleon >= 1.0 ) {
2230 InerSuccess =
false;
2233 XplusResidual = 1.0 - XplusNucleon;
2234 }
while ( ! InerSuccess );
2237 XplusResidual = 1.0;
2242 G4cout <<
"TNucleonMass PtNucleon XplusNucleon " << TNucleonMass <<
" " << PtNucleon
2243 <<
" " << XplusNucleon << G4endl
2244 <<
"TResidualMass PtResidual XplusResidual " << TResidualMass <<
" " << PtResidual
2245 <<
" " << XplusResidual <<
G4endl;
2248 M2projectile = (
sqr( TNucleonMass ) + PtNucleon.mag2() ) / XplusNucleon +
2249 (
sqr( TResidualMass ) + PtResidual.mag2() ) / XplusResidual;
2252 G4cout <<
"SqrtS < Mtarget + std::sqrt(M2projectile) " << SqrtS <<
" " << Mtarget
2253 <<
" " << std::sqrt( M2projectile ) <<
" " << Mtarget + std::sqrt( M2projectile )
2257 }
while ( SqrtS < Mtarget + std::sqrt( M2projectile ) );
2260 - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
2262 WplusProjectile = ( S + M2projectile - M2target + std::sqrt( DecayMomentum2 ) )/2.0/SqrtS;
2263 WminusTarget = SqrtS - M2projectile/WplusProjectile;
2265 G4double Pztarget = -WminusTarget/2.0 + M2target/2.0/WminusTarget;
2266 G4double Etarget = WminusTarget/2.0 + M2target/2.0/WminusTarget;
2267 G4double Ytarget = 0.5 * std::log( (Etarget + Pztarget)/(Etarget - Pztarget) );
2270 G4cout <<
"DecayMomentum2 " << DecayMomentum2 << G4endl
2271 <<
"WminusTarget WplusProjectile " << WminusTarget <<
" " << WplusProjectile
2272 << G4endl <<
"YtargetNucleon " << Ytarget <<
G4endl;
2275 G4double Mt2 =
sqr( TNucleonMass ) + PtNucleon.mag2();
2276 G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
2277 G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
2278 G4double YprojectileNucleon = 0.5 * std::log( (E + Pz)/(E - Pz) );
2281 G4cout <<
"YpN Ypr YpN-Ypr " <<
" " << YprojectileNucleon <<
" " << YprojectileNucleus
2282 <<
" " << YprojectileNucleon - YprojectileNucleus << G4endl
2283 <<
"YpN Ytr YpN-Ytr " <<
" " << YprojectileNucleon <<
" " << Ytarget
2284 <<
" " << YprojectileNucleon - Ytarget <<
G4endl;
2287 if ( std::abs( YprojectileNucleon - YprojectileNucleus ) > 2 ||
2288 Ytarget > YprojectileNucleon ) {
2289 OuterSuccess =
false;
2293 }
while ( ! OuterSuccess );
2296 G4double Pztarget = -WminusTarget/2.0 + M2target/2.0/WminusTarget;
2297 G4double Etarget = WminusTarget/2.0 + M2target/2.0/WminusTarget;
2298 Ptarget.setPz( Pztarget ); Ptarget.setE( Etarget );
2299 Ptarget.transform( toLab );
2307 G4double Mt2 =
sqr( TNucleonMass ) + PtNucleon.mag2();
2308 G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
2309 G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
2310 Pprojectile.setPx( PtNucleon.x() ); Pprojectile.setPy( PtNucleon.y() );
2311 Pprojectile.setPz( Pz ); Pprojectile.setE( E );
2312 Pprojectile.transform( toLab );
2316 G4cout <<
"Proj after in Lab " << Pprojectile <<
G4endl;
2320 ProjectileResidualMassNumber = TResidualMassNumber;
2321 ProjectileResidualCharge = TResidualCharge;
2322 ProjectileResidualExcitationEnergy = TResidualExcitationEnergy;
2324 if ( ProjectileResidualMassNumber != 0 ) {
2325 Mt2 =
sqr( TResidualMass ) + PtResidual.mag2();
2326 Pz = WplusProjectile*XplusResidual/2.0 - Mt2/(2.0*WplusProjectile*XplusResidual);
2327 E = WplusProjectile*XplusResidual/2.0 + Mt2/(2.0*WplusProjectile*XplusResidual);
2328 ProjectileResidual4Momentum.setPx( PtResidual.x() );
2329 ProjectileResidual4Momentum.setPy( PtResidual.y() );
2330 ProjectileResidual4Momentum.setPz( Pz );
2331 ProjectileResidual4Momentum.setE( E );
2332 ProjectileResidual4Momentum.transform( toLab );
2350 G4cout <<
"Proj res Init " << ProjectileResidual4Momentum << G4endl
2351 <<
"Targ res Init " << TargetResidual4Momentum << G4endl
2352 <<
"ProjectileResidualMassNumber ProjectileResidualCharge "
2353 << ProjectileResidualMassNumber <<
" " << ProjectileResidualCharge << G4endl
2354 <<
"TargetResidualMassNumber TargetResidualCharge " << TargetResidualMassNumber
2355 <<
" " << TargetResidualCharge <<
G4endl;
2358 G4LorentzVector Psum = ProjectileResidual4Momentum + TargetResidual4Momentum;
2373 G4int PResidualMassNumber = ProjectileResidualMassNumber - 1;
2374 G4int PResidualCharge = ProjectileResidualCharge -
2376 G4double PResidualExcitationEnergy = ProjectileResidualExcitationEnergy +
2377 ExcitationEnergyPerWoundedNucleon;
2378 if ( PResidualMassNumber <= 1 ) {
2379 PResidualExcitationEnergy = 0.0;
2383 if ( PResidualMassNumber != 0 ) {
2385 ->
GetIonMass( PResidualCharge, PResidualMassNumber );
2390 G4int TResidualMassNumber = TargetResidualMassNumber - 1;
2391 G4int TResidualCharge = TargetResidualCharge -
2393 G4double TResidualExcitationEnergy = TargetResidualExcitationEnergy +
2394 ExcitationEnergyPerWoundedNucleon;
2395 if ( TResidualMassNumber <= 1 ) {
2396 TResidualExcitationEnergy = 0.0;
2399 if ( TResidualMassNumber != 0 ) {
2401 ->
GetIonMass( TResidualCharge, TResidualMassNumber );
2406 G4double SumMasses = PNucleonMass + PResidualMass + TNucleonMass + TResidualMass;
2409 G4cout <<
"PNucleonMass PResidualMass TNucleonMass TResidualMass " << PNucleonMass
2410 <<
" " << PResidualMass <<
" " << TNucleonMass <<
" " << TResidualMass << G4endl
2411 <<
"PResidualExcitationEnergy " << PResidualExcitationEnergy << G4endl
2412 <<
"TResidualExcitationEnergy " << TResidualExcitationEnergy <<
G4endl;
2417 if ( ! Annihilation ) {
2420 G4cout <<
"SqrtS < SumMasses " << SqrtS <<
" " << SumMasses <<
G4endl;
2423 if ( SqrtS < SumMasses ) {
2428 G4cout <<
"SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy "
2429 << SqrtS <<
" " << SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy
2433 if ( SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy ) {
2436 if ( PResidualExcitationEnergy <= 0.0 ) {
2437 TResidualExcitationEnergy = SqrtS - SumMasses;
2438 }
else if ( TResidualExcitationEnergy <= 0.0 ) {
2439 PResidualExcitationEnergy = SqrtS - SumMasses;
2441 G4double Fraction = (SqrtS - SumMasses) /
2442 (PResidualExcitationEnergy + TResidualExcitationEnergy);
2443 PResidualExcitationEnergy *= Fraction;
2444 TResidualExcitationEnergy *= Fraction;
2453 if ( Annihilation ) {
2454 if ( SqrtS < SumMasses - TNucleonMass ) {
2457 if ( SqrtS < SumMasses ) {
2459 TNucleonMass = SqrtS - (SumMasses - TNucleonMass);
2461 TResidualExcitationEnergy = 0.0;
2463 if ( SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy ) {
2465 if ( PResidualExcitationEnergy <= 0.0 ) {
2466 TResidualExcitationEnergy = SqrtS - SumMasses;
2467 }
else if ( TResidualExcitationEnergy <= 0.0 ) {
2468 PResidualExcitationEnergy = SqrtS - SumMasses;
2470 G4double Fraction = (SqrtS - SumMasses) /
2471 (PResidualExcitationEnergy + TResidualExcitationEnergy);
2472 PResidualExcitationEnergy *= Fraction;
2473 TResidualExcitationEnergy *= Fraction;
2482 Ptmp.
setE( PNucleonMass );
2483 Pprojectile = Ptmp; Pprojectile.
transform( toLab );
2487 ProjectileResidualMassNumber = PResidualMassNumber;
2488 ProjectileResidualCharge = PResidualCharge;
2489 ProjectileResidualExcitationEnergy = PResidualExcitationEnergy;
2491 Ptmp.
setE( PResidualMass + ProjectileResidualExcitationEnergy );
2493 ProjectileResidual4Momentum = Ptmp;
2496 Ptmp.
setE( TNucleonMass );
2497 Ptarget = Ptmp; Ptarget.
transform( toLab );
2501 TargetResidualMassNumber = TResidualMassNumber;
2502 TargetResidualCharge = TResidualCharge;
2503 TargetResidualExcitationEnergy = TResidualExcitationEnergy;
2505 Ptmp.
setE( TResidualMass + TargetResidualExcitationEnergy );
2507 TargetResidual4Momentum = Ptmp;
2512 G4LorentzVector PResidual4Momentum = toCms * ProjectileResidual4Momentum;
2516 G4cout <<
"YprojectileNucleus XcenterP " << YprojectileNucleus <<
G4endl;
2522 PResidualMass += PResidualExcitationEnergy;
2523 TResidualMass += TResidualExcitationEnergy;
2540 G4int NumberOfTries( 0 );
2542 G4bool OuterSuccess(
true );
2546 OuterSuccess =
true;
2552 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
2555 DcorP *= ScaleFactor;
2556 DcorT *= ScaleFactor;
2557 AveragePt2 *= ScaleFactor;
2564 if ( ProjectileResidualMassNumber > 1 ) {
2565 PtNucleonP = GaussianPt( AveragePt2, maxPtSquare );
2569 PtResidualP = -PtNucleonP;
2571 if ( TargetResidualMassNumber > 1 ) {
2572 PtNucleonT = GaussianPt( AveragePt2, maxPtSquare );
2576 PtResidualT = -PtNucleonT;
2578 G4double Mprojectile = std::sqrt(
sqr( PNucleonMass ) + PtNucleonP.mag2() ) +
2579 std::sqrt(
sqr( PResidualMass ) + PtResidualP.mag2() );
2580 M2projectile =
sqr( Mprojectile );
2582 G4double Mtarget = std::sqrt(
sqr( TNucleonMass ) + PtNucleonT.mag2() ) +
2583 std::sqrt(
sqr( TResidualMass ) + PtResidualT.mag2() );
2584 M2target =
sqr( Mtarget );
2586 if ( SqrtS < Mprojectile + Mtarget ) {
2587 OuterSuccess =
false;
2591 G4bool InerSuccess =
true;
2593 if ( ProjectileResidualMassNumber > 1 ) {
2597 G4double XcenterP = std::sqrt(
sqr( PNucleonMass ) + PtNucleonP.mag2() ) / Mprojectile;
2598 XplusNucleon = XcenterP + tmpX.
x();
2605 if ( XplusNucleon <= 0.0 || XplusNucleon >= 1.0 ) {
2606 InerSuccess =
false;
2609 XplusResidual = 1.0 - XplusNucleon;
2610 }
while ( ! InerSuccess );
2620 XplusResidual = 1.0;
2623 if ( TargetResidualMassNumber > 1 ) {
2628 G4double XcenterT = std::sqrt(
sqr( TNucleonMass ) + PtNucleonT.mag2() ) / Mtarget;
2629 XminusNucleon = XcenterT + tmpX.
x();
2630 if ( XminusNucleon <= 0.0 || XminusNucleon >= 1.0 ) {
2631 InerSuccess =
false;
2634 XminusResidual = 1.0 - XminusNucleon;
2635 }
while ( ! InerSuccess );
2637 XminusNucleon = 1.0;
2638 XminusResidual = 1.0;
2642 G4cout <<
"PtNucleonP " << PtNucleonP <<
" " << PtResidualP << G4endl
2643 <<
"XplusNucleon XplusResidual " << XplusNucleon <<
" " << XplusResidual << G4endl
2644 <<
"PtNucleonT " << PtNucleonT <<
" " << PtResidualT << G4endl
2645 <<
"XminusNucleon XminusResidual " << XminusNucleon <<
" " << XminusResidual
2649 M2projectile = (
sqr( PNucleonMass ) + PtNucleonP.mag2() ) / XplusNucleon +
2650 (
sqr( PResidualMass) + PtResidualP.mag2() ) / XplusResidual;
2651 M2target = (
sqr( TNucleonMass ) + PtNucleonT.mag2() ) / XminusNucleon +
2652 (
sqr( TResidualMass ) + PtResidualT.mag2() ) / XminusResidual;
2654 }
while ( SqrtS < std::sqrt( M2projectile ) + std::sqrt( M2target ) );
2657 - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
2659 WplusProjectile = ( S + M2projectile - M2target + std::sqrt( DecayMomentum2 ) )/2.0/SqrtS;
2660 WminusTarget = SqrtS - M2projectile/WplusProjectile;
2662 G4double Mt2 =
sqr( PNucleonMass ) + PtNucleonP.mag2();
2663 G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
2664 G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
2665 G4double YprojectileNucleon = 0.5 * std::log( (E + Pz)/(E - Pz) );
2667 Mt2 =
sqr( TNucleonMass ) + PtNucleonT.mag2();
2668 Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2669 E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2670 G4double YtargetNucleon = 0.5 * std::log( (E + Pz)/(E - Pz) );
2672 if ( std::abs( YtargetNucleon - YtargetNucleus ) > 2 ||
2673 std::abs( YprojectileNucleon - YprojectileNucleus ) > 2 ||
2674 YprojectileNucleon < YtargetNucleon ) {
2675 OuterSuccess =
false;
2679 }
while ( ! OuterSuccess );
2685 G4double Mt2 =
sqr( PNucleonMass ) + PtNucleonP.mag2();
2686 G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
2687 G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
2689 Pprojectile.setPx( PtNucleonP.x() ); Pprojectile.setPy( PtNucleonP.y() );
2690 Pprojectile.setPz( Pz ); Pprojectile.setE( E );
2691 Pprojectile.transform( toLab );
2695 ProjectileResidualMassNumber = PResidualMassNumber;
2696 ProjectileResidualCharge = PResidualCharge;
2697 ProjectileResidualExcitationEnergy = PResidualExcitationEnergy;
2700 G4cout <<
"PResidualMass PtResidualP " << PResidualMass <<
" " << PtResidualP <<
G4endl;
2703 if ( ProjectileResidualMassNumber != 0 ) {
2704 Mt2 =
sqr( PResidualMass ) + PtResidualP.mag2();
2705 Pz = WplusProjectile*XplusResidual/2.0 - Mt2/(2.0*WplusProjectile*XplusResidual);
2706 E = WplusProjectile*XplusResidual/2.0 + Mt2/(2.0*WplusProjectile*XplusResidual);
2707 ProjectileResidual4Momentum.setPx( PtResidualP.x() );
2708 ProjectileResidual4Momentum.setPy( PtResidualP.y() );
2709 ProjectileResidual4Momentum.setPz( Pz );
2710 ProjectileResidual4Momentum.setE( E );
2711 ProjectileResidual4Momentum.transform( toLab );
2717 G4cout <<
"Pr N R " << Pprojectile << G4endl <<
" "
2718 << ProjectileResidual4Momentum <<
G4endl;
2721 Mt2 =
sqr( TNucleonMass ) + PtNucleonT.mag2();
2722 Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2723 E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2725 Ptarget.
setPx( PtNucleonT.x() ); Ptarget.
setPy( PtNucleonT.y() );
2731 TargetResidualMassNumber = TResidualMassNumber;
2732 TargetResidualCharge = TResidualCharge;
2733 TargetResidualExcitationEnergy = TResidualExcitationEnergy;
2735 if ( TargetResidualMassNumber != 0 ) {
2736 Mt2 =
sqr( TResidualMass ) + PtResidualT.mag2();
2737 Pz = -WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
2738 E = WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
2740 TargetResidual4Momentum.setPx( PtResidualT.x() );
2741 TargetResidual4Momentum.setPy( PtResidualT.y() );
2742 TargetResidual4Momentum.setPz( Pz );
2743 TargetResidual4Momentum.setE( E) ;
2744 TargetResidual4Momentum.transform( toLab );
2750 G4cout <<
"Tr N R " << Ptarget << G4endl <<
" " << TargetResidual4Momentum <<
G4endl;
2772 std::vector< G4VSplitableHadron* > primaries;
2774 while ( theParticipants.
Next() ) {
2778 if ( primaries.end() == std::find( primaries.begin(), primaries.end(),
2785 #ifdef debugBuildString
2786 G4cout <<
"G4FTFModel::BuildStrings()" << G4endl
2787 <<
"Number of projectile strings " << primaries.size() <<
G4endl;
2790 for (
unsigned int ahadron = 0; ahadron < primaries.size(); ahadron++ ) {
2791 G4bool isProjectile(
true );
2794 FirstString = 0; SecondString = 0;
2795 theExcitation->
CreateStrings( primaries[ ahadron ], isProjectile,
2796 FirstString, SecondString, theParameters );
2797 if ( FirstString != 0 ) strings->push_back( FirstString );
2798 if ( SecondString != 0 ) strings->push_back( SecondString );
2800 #ifdef debugBuildString
2801 G4cout <<
"FirstString & SecondString? " << FirstString <<
" " << SecondString << G4endl
2802 <<
"Quarks on the FirstString ends " << FirstString->GetRightParton()->GetPDGcode()
2803 <<
" " << FirstString->GetLeftParton()->GetPDGcode() <<
G4endl;
2808 #ifdef debugBuildString
2809 G4cout <<
"Check 1 string " << strings->operator[](0)->GetRightParton()->GetPDGcode()
2810 <<
" " << strings->operator[](0)->GetLeftParton()->GetPDGcode() << G4endl <<
G4endl;
2818 #ifdef debugBuildString
2819 G4cout <<
"Building of projectile-like strings" <<
G4endl;
2822 G4bool isProjectile =
true;
2823 for (
G4int ahadron = 0; ahadron < NumberOfInvolvedNucleonsOfProjectile; ahadron++ ) {
2825 #ifdef debugBuildString
2826 G4cout <<
"Nucleon #, status, intCount " << ahadron <<
" "
2835 #ifdef debugBuildString
2836 G4cout << G4endl <<
"ahadron aProjectile Status " << ahadron <<
" " << aProjectile
2842 #ifdef debugBuildString
2843 G4cout <<
"Case1 aProjectile->GetStatus() == 0 " <<
G4endl;
2846 FirstString = 0; SecondString = 0;
2848 TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron(),
2849 isProjectile, FirstString, SecondString, theParameters );
2850 if ( FirstString != 0 ) strings->push_back( FirstString );
2851 if ( SecondString != 0 ) strings->push_back( SecondString );
2855 #ifdef debugBuildString
2856 G4cout <<
"Case2 aProjectile->GetStatus() !=0 St==1 SoftCol!=0" <<
G4endl;
2859 FirstString = 0; SecondString = 0;
2861 TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron(),
2862 isProjectile, FirstString, SecondString, theParameters );
2863 if ( FirstString != 0 ) strings->push_back( FirstString );
2864 if ( SecondString != 0 ) strings->push_back( SecondString );
2871 #ifdef debugBuildString
2872 G4cout <<
"Case3 aProjectile->GetStatus() !=0 St==1 SoftCol==0" <<
G4endl;
2875 FirstString = 0; SecondString = 0;
2877 TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron(),
2878 isProjectile, FirstString, SecondString, theParameters );
2879 if ( FirstString != 0 ) strings->push_back( FirstString );
2880 if ( SecondString != 0 ) strings->push_back( SecondString );
2882 #ifdef debugBuildString
2883 G4cout <<
" Strings are built for nucleon marked for an interaction, but"
2884 <<
" the interaction was skipped." <<
G4endl;
2887 }
else if ( aProjectile->
GetStatus() == 2 ) {
2890 #ifdef debugBuildString
2891 G4cout <<
"Case4 aProjectile->GetStatus() !=0 St==2 " <<
G4endl;
2894 FirstString = 0; SecondString = 0;
2896 TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron(),
2897 isProjectile, FirstString, SecondString, theParameters );
2898 if ( FirstString != 0 ) strings->push_back( FirstString );
2899 if ( SecondString != 0 ) strings->push_back( SecondString );
2901 #ifdef debugBuildString
2902 G4cout <<
" Strings are build for involved nucleon." <<
G4endl;
2907 #ifdef debugBuildString
2914 #ifdef debugBuildString
2922 #ifdef debugBuildString
2926 G4bool isProjectile =
false;
2927 for (
G4int ahadron = 0; ahadron < NumberOfInvolvedNucleonsOfTarget; ahadron++ ) {
2930 #ifdef debugBuildString
2931 G4cout <<
"Nucleon #, status, intCount " << aNucleon <<
" " << ahadron <<
" "
2936 FirstString = 0 ; SecondString = 0;
2938 FirstString, SecondString, theParameters );
2939 if ( FirstString != 0 ) strings->push_back( FirstString );
2940 if ( SecondString != 0 ) strings->push_back( SecondString );
2942 #ifdef debugBuildString
2948 FirstString = 0; SecondString = 0;
2950 FirstString, SecondString, theParameters );
2951 if ( FirstString != 0 ) strings->push_back( FirstString );
2952 if ( SecondString != 0 ) strings->push_back( SecondString );
2954 #ifdef debugBuildString
2955 G4cout <<
"2 case A string is build, nucleon was excited." <<
G4endl;
2963 FirstString = 0; SecondString = 0;
2965 FirstString, SecondString, theParameters );
2966 if ( FirstString != 0 ) strings->push_back( FirstString );
2967 if ( SecondString != 0 ) strings->push_back( SecondString );
2969 #ifdef debugBuildString
2974 ! HighEnergyInter ) {
2981 #ifdef debugBuildString
2985 }
else if ( aNucleon->
GetStatus() == 2 ) {
2986 FirstString = 0; SecondString = 0;
2988 FirstString, SecondString, theParameters );
2989 if ( FirstString != 0 ) strings->push_back( FirstString );
2990 if ( SecondString != 0 ) strings->push_back( SecondString );
2992 #ifdef debugBuildString
2998 #ifdef debugBuildString
3005 #ifdef debugBuildString
3006 G4cout << G4endl <<
"theAdditionalString.size() " << theAdditionalString.size()
3010 isProjectile =
true;
3011 if ( theAdditionalString.size() != 0 ) {
3012 for (
unsigned int ahadron = 0; ahadron < theAdditionalString.size(); ahadron++ ) {
3014 FirstString = 0; SecondString = 0;
3015 theExcitation->
CreateStrings( theAdditionalString[ ahadron ], isProjectile,
3016 FirstString, SecondString, theParameters );
3017 if ( FirstString != 0 ) strings->push_back( FirstString );
3018 if ( SecondString != 0 ) strings->push_back( SecondString );
3034 void G4FTFModel::GetResiduals() {
3037 #ifdef debugFTFmodel
3038 G4cout <<
"GetResiduals(): HighEnergyInter? GetProjectileNucleus()?"
3042 if ( HighEnergyInter ) {
3044 #ifdef debugFTFmodel
3045 G4cout <<
"NumberOfInvolvedNucleonsOfTarget "<< NumberOfInvolvedNucleonsOfTarget <<
G4endl;
3048 G4double DeltaExcitationE = TargetResidualExcitationEnergy /
3049 G4double( NumberOfInvolvedNucleonsOfTarget );
3051 G4double( NumberOfInvolvedNucleonsOfTarget );
3053 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
3054 G4Nucleon* aNucleon = TheInvolvedNucleonsOfTarget[i];
3056 #ifdef debugFTFmodel
3069 #ifdef debugFTFmodel
3070 G4cout <<
"NumberOfInvolvedNucleonsOfProjectile " << NumberOfInvolvedNucleonsOfProjectile
3071 << G4endl <<
"ProjectileResidualExcitationEnergy ProjectileResidual4Momentum "
3072 << ProjectileResidualExcitationEnergy <<
" " << ProjectileResidual4Momentum <<
G4endl;
3075 DeltaExcitationE = ProjectileResidualExcitationEnergy /
3076 G4double( NumberOfInvolvedNucleonsOfProjectile );
3077 DeltaPResidualNucleus = ProjectileResidual4Momentum /
3078 G4double( NumberOfInvolvedNucleonsOfProjectile );
3080 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
3081 G4Nucleon* aNucleon = TheInvolvedNucleonsOfProjectile[i];
3083 #ifdef debugFTFmodel
3094 #ifdef debugFTFmodel
3100 #ifdef debugFTFmodel
3101 G4cout <<
"Low energy interaction: Target nucleus --------------" << G4endl
3102 <<
"Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
3103 << TargetResidualMassNumber <<
" " << TargetResidualCharge <<
" "
3104 << TargetResidualExcitationEnergy <<
G4endl;
3107 G4int NumberOfTargetParticipant( 0 );
3108 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
3109 G4Nucleon* aNucleon = TheInvolvedNucleonsOfTarget[i];
3117 if ( NumberOfTargetParticipant != 0 ) {
3118 DeltaExcitationE = TargetResidualExcitationEnergy /
G4double( NumberOfTargetParticipant );
3119 DeltaPResidualNucleus = TargetResidual4Momentum /
G4double( NumberOfTargetParticipant );
3122 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
3123 G4Nucleon* aNucleon = TheInvolvedNucleonsOfTarget[i];
3130 delete targetSplitable;
3131 targetSplitable = 0;
3132 aNucleon->
Hit( targetSplitable );
3137 #ifdef debugFTFmodel
3138 G4cout <<
"NumberOfTargetParticipant " << NumberOfTargetParticipant << G4endl
3139 <<
"TargetResidual4Momentum " << TargetResidual4Momentum <<
G4endl;
3144 #ifdef debugFTFmodel
3145 G4cout <<
"Low energy interaction: Projectile nucleus --------------" << G4endl
3146 <<
"Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
3147 << ProjectileResidualMassNumber <<
" " << ProjectileResidualCharge <<
" "
3148 << ProjectileResidualExcitationEnergy <<
G4endl;
3151 G4int NumberOfProjectileParticipant( 0 );
3152 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
3153 G4Nucleon* aNucleon = TheInvolvedNucleonsOfProjectile[i];
3156 NumberOfProjectileParticipant++;
3159 #ifdef debugFTFmodel
3163 DeltaExcitationE = 0.0;
3166 if ( NumberOfProjectileParticipant != 0 ) {
3167 DeltaExcitationE = ProjectileResidualExcitationEnergy /
3168 G4double( NumberOfProjectileParticipant );
3169 DeltaPResidualNucleus = ProjectileResidual4Momentum /
3170 G4double( NumberOfProjectileParticipant );
3174 for (
G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
3175 G4Nucleon* aNucleon = TheInvolvedNucleonsOfProjectile[i];
3182 delete projectileSplitable;
3183 projectileSplitable = 0;
3184 aNucleon->
Hit( projectileSplitable );
3189 #ifdef debugFTFmodel
3190 G4cout <<
"NumberOfProjectileParticipant " << NumberOfProjectileParticipant << G4endl
3191 <<
"ProjectileResidual4Momentum " << ProjectileResidual4Momentum <<
G4endl;
3196 #ifdef debugFTFmodel
3197 G4cout <<
"End GetResiduals -----------------" <<
G4endl;
3209 if ( AveragePt2 <= 0.0 ) {
3213 ( std::exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
3218 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
3225 desc <<
"please add description here" <<
G4endl;
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
void GetList(const G4ReactionProduct &thePrimary, G4FTFParameters *theParameters)
Hep3Vector boostVector() const
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)
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)
virtual void DoLorentzContraction(const G4LorentzVector &theBoost)=0
virtual G4int GetMassNumber()=0
void SetTimeOfCreation(G4double aTime)
G4int GetPDGEncoding() const
virtual const G4LorentzVector & Get4Momentum() const
G4V3DNucleus * GetTargetNucleus() const
virtual const G4ThreeVector & GetPosition() const
virtual void SetProjectileNucleus(G4V3DNucleus *aNucleus)
HepLorentzVector & rotateZ(double)
void SetStatus(const G4int aStatus)
void SetDefinition(G4ParticleDefinition *aDefinition)
const G4String & GetParticleName() const
G4ParticleDefinition * GetDefinition() const
void SetThisPointer(G4VPartonStringModel *aPointer)
G4VSplitableHadron * GetSplitableHadron() const
G4ParticleDefinition * GetDefinition() const
G4InteractionContent & GetInteraction()
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
G4double GetMaxNumberOfCollisions()
virtual void Init(G4int theZ, G4int theA)
static G4AntiProton * AntiProton()
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double AbsoluteLevel)
G4V3DNucleus * GetProjectileNucleus() const
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
void SetTotalEnergy(const G4double en)
G4double GetTimeOfCreation()
G4V3DNucleus * theProjectileNucleus
static G4Proton * Proton()
HepLorentzRotation & transform(const HepBoost &b)
G4double GetMaxPt2ofNuclearDestruction()
static G4Neutron * Neutron()
void SetBindingEnergy(G4double anEnergy)
void Init(const G4Nucleus &aNucleus, const G4DynamicParticle &aProjectile)
const G4LorentzVector & Get4Momentum() const
virtual void ModelDescription(std::ostream &) const
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
G4double GetTotalEnergy() 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()
G4double GetR2ofNuclearDestruction()
G4double GetExcitationEnergyPerWoundedNucleon()
HepLorentzVector & rotateY(double)
G4VSplitableHadron * GetTarget() const
G4VSplitableHadron * GetProjectile() const
G4ThreeVector GetMomentum() const
virtual G4bool Annihilate(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4VSplitableHadron *&AdditionalString, G4FTFParameters *theParameters) const
G4double GetCofNuclearDestruction()
virtual G4ParticleDefinition * GetDefinition() const
HepLorentzVector & transform(const HepRotation &)
G4double GetProbabilityOfAnnihilation()
virtual G4Nucleon * GetNextNucleon()=0
void Set4Momentum(const G4LorentzVector &a4Momentum)
void Hit(G4VSplitableHadron *aHit)
void setVect(const Hep3Vector &)
G4double GetPDGCharge() const
G4ExcitedStringVector * GetStrings()
static G4AntiNeutron * AntiNeutron()
G4int GetBaryonNumber() const
CLHEP::HepLorentzVector G4LorentzVector
G4double GetDofNuclearDestruction()