50 : theDiffExcitaton(), ModelMode(SOFT), nCutMax(7),
51 ThresholdParameter(0.0*
GeV), QGSMThreshold(3.0*
GeV),
53 theProjectileSplitable(nullptr), Regge(nullptr),
54 InteractionMode(
ALL),
alpha(-0.5),
beta(2.5), sigmaPt(0.0),
55 NumberOfInvolvedNucleonsOfTarget(0), NumberOfInvolvedNucleonsOfProjectile(0),
56 ProjectileResidual4Momentum(
G4LorentzVector()), ProjectileResidualMassNumber(0),
57 ProjectileResidualCharge(0), ProjectileResidualExcitationEnergy(0.0),
58 TargetResidual4Momentum(
G4LorentzVector()), TargetResidualMassNumber(0),
59 TargetResidualCharge(0), TargetResidualExcitationEnergy(0.0),
60 CofNuclearDestruction(0.0), R2ofNuclearDestruction(0.0),
61 ExcitationEnergyPerWoundedNucleon(0.0), DofNuclearDestruction(0.0),
62 Pt2ofNuclearDestruction(0.0), MaxPt2ofNuclearDestruction(0.0)
64 for (
size_t i=0; i < 250; ++i) {
80 :
G4VParticipants(), ModelMode(right.ModelMode), nCutMax(right.nCutMax),
81 ThresholdParameter(right.ThresholdParameter),
82 QGSMThreshold(right.QGSMThreshold),
83 theNucleonRadius(right.theNucleonRadius),
84 theCurrentVelocity(right.theCurrentVelocity),
85 theProjectileSplitable(right.theProjectileSplitable),
86 Regge(right.Regge), InteractionMode(right.InteractionMode),
88 NumberOfInvolvedNucleonsOfTarget(right.NumberOfInvolvedNucleonsOfTarget),
89 NumberOfInvolvedNucleonsOfProjectile(right.NumberOfInvolvedNucleonsOfProjectile),
90 ProjectileResidual4Momentum(right.ProjectileResidual4Momentum),
91 ProjectileResidualMassNumber(right.ProjectileResidualMassNumber),
92 ProjectileResidualCharge(right.ProjectileResidualCharge),
93 ProjectileResidualExcitationEnergy(right.ProjectileResidualExcitationEnergy),
94 TargetResidual4Momentum(right.TargetResidual4Momentum),
95 TargetResidualMassNumber(right.TargetResidualMassNumber),
96 TargetResidualCharge(right.TargetResidualCharge),
97 TargetResidualExcitationEnergy(right.TargetResidualExcitationEnergy),
98 CofNuclearDestruction(right.CofNuclearDestruction),
99 R2ofNuclearDestruction(right.R2ofNuclearDestruction),
100 ExcitationEnergyPerWoundedNucleon(right.ExcitationEnergyPerWoundedNucleon),
101 DofNuclearDestruction(right.DofNuclearDestruction),
102 Pt2ofNuclearDestruction(right.Pt2ofNuclearDestruction),
103 MaxPt2ofNuclearDestruction(right.MaxPt2ofNuclearDestruction)
105 for (
size_t i=0; i < 250; ++i) {
149 #ifdef debugQGSParticipants
159 const G4int maxNumberOfLoops = 1000;
160 G4int loopCounter = 0;
163 const G4int maxNumberOfInternalLoops = 1000;
164 G4int internalLoopCounter = 0;
186 }
while( (!Success) && ++internalLoopCounter < maxNumberOfInternalLoops );
188 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
193 #ifdef debugQGSParticipants
199 #ifdef debugQGSParticipants
206 #ifdef debugQGSParticipants
207 G4cout<<
"Perform non-Diffractive Collisions if they happend. Determine Parton Momenta and so on." <<
G4endl;
214 }
while( (!Success) && ++loopCounter < maxNumberOfLoops );
216 if ( loopCounter >= maxNumberOfLoops ) {
218 #ifdef debugQGSParticipants
228 #ifdef debugQGSParticipants
234 #ifdef debugQGSParticipants
245 #ifdef debugQGSParticipants
252 if ( (aNucleon != 0 ) && (aNucleon->
GetStatus() >= 1) )
delete aNucleon;
260 if ( aNucleon )
delete aNucleon;
264 #ifdef debugQGSParticipants
265 G4cout<<
"Delition of target nucleons from soft interactions "<<
theTargets.size()
284 if( pProjectile )
delete pProjectile;
296 if ( (splaNucleon != 0) && (splaNucleon->
GetStatus() >=1) )
delete splaNucleon;
297 aNucleon->
Hit(
nullptr);
336 #ifdef debugQGSParticipants
347 G4double SS=(aPrimaryMomentum + aNucleonMomentum).mag2();
356 #ifdef debugQGSParticipants
357 G4cout <<
"QGSM - BAD situation: pNucleon is NULL ! Leaving immediately!" <<
G4endl;
368 theNucleusOuterR = 0.;
378 const G4int maxNumberOfLoops = 1000;
380 G4int NumberOfTries = 0;
383 G4int loopCounter = -1;
384 while( (
theInteractions.size() == 0) && ++loopCounter < maxNumberOfLoops )
389 std::pair<G4double, G4double> theImpactParameter;
392 if( NumberOfTries == 100*(NumberOfTries/100) ) Scale /=2.0;
395 G4double impactX = theImpactParameter.first;
396 G4double impactY = theImpactParameter.second;
398 #ifdef debugQGSParticipants
405 G4int nucleonCount = -1;
412 if(Power <= 0.)
break;
421 G4double Pprd(0.), Ptrd(0.), Pdd(0.);
423 G4int NcutPomerons(0);
426 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr);
427 #ifdef debugQGSParticipants
428 G4cout<<
"Nucleon & its impact parameter: "<<nucleonCount<<
" "<<std::sqrt(Distance2)/
fermi<<
" (fm)"<<
G4endl;
430 G4cout<<
"Probability of PrD, TrD, DD: "<<Pprd<<
" "<<Ptrd<<
" "<<Pdd<<
G4endl;
431 G4cout<<
"Probability of NonDiff, QuarkExc.: "<<Pnd<<
" "<<Pnvr<<
" in inel. inter."<<
G4endl;
438 G4int InteractionType(0);
452 if( rndNumber < Ptrd ) {InteractionType =
TrD; }
453 else if( rndNumber < Ptrd + Pnd) {InteractionType =
NonD; NcutPomerons =
Regge->
ncPomerons();}
456 if( (InteractionType ==
NonD) && (NcutPomerons == 0))
continue;
460 tNucleon->
Hit(aTargetSPB);
462 #ifdef debugQGSParticipants
464 G4cout<<
"Target nucleon - "<<nucleonCount<<
" "
466 G4cout<<
"Interaction type:"<<InteractionType
467 <<
" (0 -PrD, 1 - TrD, 2 - DD, 3 - NonD, 4 - Qexc)"<<
G4endl;
469 <<
" (0 -ALL, 1 - WITHOUT_R, 2 - NON_DIFF)"<<
G4endl;
470 if( InteractionType ==
NonD )
471 G4cout<<
"Number of cutted pomerons: "<<NcutPomerons<<
G4endl;
474 if((InteractionType ==
PrD) || (InteractionType ==
TrD) || (InteractionType ==
DD) ||
475 (InteractionType ==
Qexc))
477 #ifdef debugQGSParticipants
491 aInteraction->
SetStatus(InteractionType);
496 #ifdef debugQGSParticipants
497 G4cout<<
"Non-diffractive interaction occurs, max NcutPomerons "<<NcutPomerons<<
G4endl;
503 for(nCuts = 0; nCuts < NcutPomerons; nCuts++)
505 if(
G4UniformRand() < Power/MaxPower ){Vncut++; Power--;
if(Power <= 0.)
break;}
509 if( nCuts == 0 ) {
delete aTargetSPB; tNucleon->
Hit(
nullptr);
continue;}
511 #ifdef debugQGSParticipants
512 G4cout<<
"Number of cuts in the interaction "<<nCuts<<
G4endl;
526 aInteraction->
SetStatus(InteractionType);
532 #ifdef debugQGSParticipants
538 if ( loopCounter >= maxNumberOfLoops ) {
539 #ifdef debugQGSParticipants
548 std::vector<G4InteractionContent*>::iterator i;
576 aTargetNucleon->
Hit(
nullptr);
607 #ifdef debugQGSParticipants
609 <<
"Stored # of wounded nucleons of target "
619 #ifdef debugQGSParticipants
628 for (
G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
651 #ifdef debugQGSParticipants
652 G4cout<<
"Target nucleon involved in reggeon cascading No "<<TrgNuc<<
" "
660 Neighbour->
Hit( targetSplitable );
666 anInteraction->
SetTarget(targetSplitable);
678 #ifdef debugQGSParticipants
688 G4bool isProjectileNucleus =
false;
690 isProjectileNucleus =
true;
693 #ifdef debugPutOnMassShell
695 if ( isProjectileNucleus ) {
G4cout <<
"PutOnMassShell for Nucleus_Nucleus " <<
G4endl;}
699 if ( Pprojectile.
z() < 0.0 ) {
711 #ifdef debugPutOnMassShell
719 if ( ! isOk )
return false;
728 if ( ! isProjectileNucleus ) {
729 Mprojectile = Pprojectile.
mag();
730 M2projectile = Pprojectile.
mag2();
731 SumMasses += Mprojectile + 20.0*
MeV;
734 #ifdef debugPutOnMassShell
735 G4cout <<
"Projectile : ";
741 if ( ! isOk )
return false;
748 #ifdef debugPutOnMassShell
752 <<
"SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/
GeV <<
" "
753 << PrResidualMass/
GeV <<
" " << TargetResidualMass/
GeV <<
" GeV" <<
G4endl;
757 if ( SqrtS < SumMasses ) {
764 G4double savedSumMasses = SumMasses;
765 if ( isProjectileNucleus ) {
766 SumMasses -= std::sqrt(
sqr( PrResidualMass ) + PprojResidual.
perp2() );
768 + PprojResidual.
perp2() );
770 SumMasses -= std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.
perp2() );
772 + PtargetResidual.
perp2() );
774 if ( SqrtS < SumMasses ) {
775 SumMasses = savedSumMasses;
776 if ( isProjectileNucleus ) {
783 if ( isProjectileNucleus ) {
787 #ifdef debugPutOnMassShell
788 if ( isProjectileNucleus ) {
789 G4cout <<
"PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/
GeV <<
" "
792 G4cout <<
"TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/
GeV <<
" GeV "
794 <<
"Sum masses " << SumMasses/
GeV <<
G4endl;
798 if ( isProjectileNucleus && thePrNucleus->
GetMassNumber() != 1 ) {
807 if ( ! isOk )
return false;
820 if ( Ptmp.
pz() <= 0.0 ) {
827 if ( isProjectileNucleus ) {
829 YprojectileNucleus = Ptmp.
rapidity();
831 Ptmp = toCms*Ptarget;
836 if ( isProjectileNucleus ) {
843 #ifdef debugPutOnMassShell
844 if ( isProjectileNucleus ) {
845 G4cout <<
"Y projectileNucleus " << YprojectileNucleus <<
G4endl;
847 G4cout <<
"Y targetNucleus " << YtargetNucleus <<
G4endl
849 <<
" DcorP DcorT " << DcorP <<
" " << DcorT <<
" AveragePt2 " << AveragePt2 <<
G4endl;
856 G4int NumberOfTries = 0;
858 G4bool OuterSuccess =
true;
860 const G4int maxNumberOfLoops = 1000;
861 G4int loopCounter = 0;
863 G4double sqrtM2proj = 0.0, sqrtM2target = 0.0;
865 const G4int maxNumberOfTries = 1000;
868 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
874 DcorP *= ScaleFactor;
875 DcorT *= ScaleFactor;
876 AveragePt2 *= ScaleFactor;
878 if ( isProjectileNucleus ) {
881 thePrNucleus, PprojResidual,
889 theTargetNucleus, PtargetResidual,
894 if ( M2proj < 0.0 ) {
895 if( M2proj < -0.000001 ) {
899 <<
") M2proj=" << M2proj <<
" -> sets it to 0.0 !" <<
G4endl;
900 G4Exception(
"G4QGSParticipants::PutOnMassShell(): negative projectile squared mass!",
905 sqrtM2proj = std::sqrt( M2proj );
906 if ( M2target < 0.0 ) {
910 <<
") M2target=" << M2target <<
" -> sets it to 0.0 !" <<
G4endl;
911 G4Exception(
"G4QGSParticipants::PutOnMassShell(): negative target squared mass!",
915 sqrtM2target = std::sqrt( M2target );
917 #ifdef debugPutOnMassShell
918 G4cout <<
"SqrtS, Mp+Mt, Mp, Mt " << SqrtS/
GeV <<
" "
919 << ( sqrtM2proj + sqrtM2target )/
GeV <<
" "
923 if ( ! isOk )
return false;
924 }
while ( ( SqrtS < ( sqrtM2proj + sqrtM2target ) ) &&
925 ++NumberOfTries < maxNumberOfTries );
926 if ( NumberOfTries >= maxNumberOfTries ) {
929 if ( isProjectileNucleus ) {
930 isOk =
CheckKinematics(
S, SqrtS, M2proj, M2target, YprojectileNucleus,
true,
933 WminusTarget, WplusProjectile, OuterSuccess );
938 WminusTarget, WplusProjectile, OuterSuccess );
939 if ( ! isOk )
return false;
940 }
while ( ( ! OuterSuccess ) &&
941 ++loopCounter < maxNumberOfLoops );
942 if ( loopCounter >= maxNumberOfLoops ) {
952 if ( ! isProjectileNucleus ) {
954 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
955 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
956 Pprojectile.
setPz( Pzprojectile );
957 Pprojectile.
setE( Eprojectile );
959 #ifdef debugPutOnMassShell
969 #ifdef debugPutOnMassShell
980 #ifdef debugPutOnMassShell
984 if ( ! isOk )
return false;
988 #ifdef debugPutOnMassShell
998 #ifdef debugPutOnMassShell
1002 if ( ! isOk )
return false;
1006 #ifdef debugPutOnMassShell
1020 if ( AveragePt2 > 0.0 ) {
1021 G4double x = maxPtSquare/AveragePt2;
1025 Pt = std::sqrt( Pt2 );
1029 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
1038 G4double& residualExcitationEnergy,
1040 G4int& residualMassNumber,
1041 G4int& residualCharge ) {
1053 if ( ! nucleus )
return false;
1078 sumMasses += 20.0*
MeV;
1082 residualMassNumber--;
1089 #ifdef debugPutOnMassShell
1090 G4cout <<
"ExcitationEnergyPerWoundedNucleon " << ExcitationEPerWoundedNucleon <<
" MeV"<<
G4endl
1091 <<
"\t Residual Charge, MassNumber " << residualCharge <<
" " << residualMassNumber
1092 <<
G4endl <<
"\t Initial Momentum " << nucleusMomentum/
GeV<<
" GeV"
1093 <<
G4endl <<
"\t Residual Momentum " << residualMomentum/
GeV<<
" GeV"<<
G4endl;
1096 residualMomentum.
setPz( 0.0 );
1097 residualMomentum.
setE( 0.0 );
1098 if ( residualMassNumber == 0 ) {
1100 residualExcitationEnergy = 0.0;
1103 GetIonMass( residualCharge, residualMassNumber );
1104 if ( residualMassNumber == 1 ) {
1105 residualExcitationEnergy = 0.0;
1108 sumMasses += std::sqrt(
sqr( residualMass ) + residualMomentum.
perp2() );
1117 const G4int numberOfInvolvedNucleons,
1135 if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 )
return false;
1139 const G4double probDeltaIsobar = 0.10;
1141 G4int maxNumberOfDeltas =
G4int( (sqrtS - sumMasses)/(400.0*
MeV) );
1142 G4int numberOfDeltas = 0;
1144 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1147 if (
G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
1149 if ( ! involvedNucleons[i] )
continue;
1156 G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4;
1165 if ( sqrtS < sumMasses + massDelta - massNuc ) {
1169 sumMasses += ( massDelta - massNuc );
1188 const G4int residualMassNumber,
1189 const G4int numberOfInvolvedNucleons,
1203 if ( ! nucleus )
return false;
1205 if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
1213 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1214 G4Nucleon* aNucleon = involvedNucleons[i];
1215 if ( ! aNucleon )
continue;
1219 const G4int maxNumberOfLoops = 1000;
1220 G4int loopCounter = 0;
1227 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1228 G4Nucleon* aNucleon = involvedNucleons[i];
1229 if ( ! aNucleon )
continue;
1235 if ( x < 0.0 || x > 1.0 ) {
1245 if ( xSum < 0.0 || xSum > 1.0 ) success =
false;
1247 if ( ! success )
continue;
1249 G4double deltaPx = ( ptSum.
x() - pResidual.
x() ) / numberOfInvolvedNucleons;
1250 G4double deltaPy = ( ptSum.
y() - pResidual.
y() ) / numberOfInvolvedNucleons;
1252 if ( residualMassNumber == 0 ) {
1253 delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
1260 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1261 G4Nucleon* aNucleon = involvedNucleons[i];
1262 if ( ! aNucleon )
continue;
1265 if ( residualMassNumber == 0 ) {
1266 if ( x <= 0.0 || x > 1.0 ) {
1271 if ( x <= 0.0 || x > 1.0 || xSum <= 0.0 || xSum > 1.0 ) {
1279 +
sqr( px ) +
sqr( py ) ) / x;
1284 if ( success && residualMassNumber != 0 ) {
1285 mass2 += (
sqr( residualMass ) + pResidual.
perp2() ) / xSum;
1288 #ifdef debugPutOnMassShell
1292 }
while ( ( ! success ) &&
1293 ++loopCounter < maxNumberOfLoops );
1294 if ( loopCounter >= maxNumberOfLoops ) {
1310 const G4bool isProjectileNucleus,
1311 const G4int numberOfInvolvedNucleons,
1326 G4double decayMomentum2 =
sqr( sValue ) +
sqr( projectileMass2 ) +
sqr( targetMass2 )
1327 - 2.0*sValue*projectileMass2 - 2.0*sValue*targetMass2
1328 - 2.0*projectileMass2*targetMass2;
1329 targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
1331 projectileWplus = sqrtS - targetMass2/targetWminus;
1332 G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
1333 G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
1335 if (projectileE - projectilePz > 0.) {
1336 projectileY = 0.5 *
G4Log( (projectileE + projectilePz)/
1337 (projectileE - projectilePz) );
1339 G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
1340 G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
1341 G4double targetY = 0.5 *
G4Log( (targetE + targetPz)/(targetE - targetPz) );
1343 #ifdef debugPutOnMassShell
1344 G4cout <<
"decayMomentum2 " << decayMomentum2 <<
G4endl
1345 <<
"\t targetWminus projectileWplus " << targetWminus <<
" " << projectileWplus <<
G4endl
1346 <<
"\t projectileY targetY " << projectileY <<
" " << targetY <<
G4endl;
1349 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1350 G4Nucleon* aNucleon = involvedNucleons[i];
1351 if ( ! aNucleon )
continue;
1356 G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*x);
1357 G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*x);
1358 if ( isProjectileNucleus ) {
1359 pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*x);
1360 e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*x);
1364 #ifdef debugPutOnMassShell
1365 G4cout <<
"i nY pY nY-AY AY " << i <<
" " << nucleonY <<
" " << projectileY <<
G4endl;
1368 if ( std::abs( nucleonY - nucleusY ) > 2 ||
1369 ( isProjectileNucleus && targetY > nucleonY ) ||
1370 ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
1383 const G4bool isProjectileNucleus,
1386 const G4int residualMassNumber,
1387 const G4int numberOfInvolvedNucleons,
1405 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1406 G4Nucleon* aNucleon = involvedNucleons[i];
1407 if ( ! aNucleon )
continue;
1409 residual3Momentum -= tmp.
vect();
1413 G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w * x );
1414 G4double e = w * x / 2.0 + mt2 / ( 2.0 * w * x );
1416 if ( isProjectileNucleus ) pz *= -1.0;
1423 #ifdef debugPutOnMassShell
1424 G4cout <<
"Target involved nucleon No, name, 4Mom "
1429 G4double residualMt2 =
sqr( residualMass ) +
sqr( residual3Momentum.
x() )
1430 +
sqr( residual3Momentum.
y() );
1432 #ifdef debugPutOnMassShell
1433 G4cout <<
G4endl<<
"w residual3Momentum.z() " << w <<
" " << residual3Momentum.
z() <<
G4endl;
1438 if ( residualMassNumber != 0 ) {
1439 residualPz = -w * residual3Momentum.
z() / 2.0 +
1440 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
1441 residualE = w * residual3Momentum.
z() / 2.0 +
1442 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
1444 if ( isProjectileNucleus ) residualPz *= -1.0;
1447 residual4Momentum.
setPx( residual3Momentum.
x() );
1448 residual4Momentum.
setPy( residual3Momentum.
y() );
1449 residual4Momentum.
setPz( residualPz );
1450 residual4Momentum.
setE( residualE );
1458 #ifdef debugQGSParticipants
1467 #ifdef debugQGSParticipants
1468 G4cout<<
"Interaction # and its status "
1473 if ( (InterStatus ==
PrD) || (InterStatus ==
TrD) || (InterStatus ==
DD))
1475 #ifdef debugQGSParticipants
1476 G4cout<<
"Simulation of diffractive interaction. "<<InterStatus <<
" PrD/TrD/DD/ND/Qech - 0,1,2,3,4"<<
G4endl;
1481 #ifdef debugQGSParticipants
1482 G4cout<<
"The proj. before inter "
1489 if ( InterStatus ==
PrD )
1492 if ( InterStatus ==
TrD )
1495 if ( InterStatus ==
DD )
1498 #ifdef debugQGSParticipants
1506 if ( InterStatus ==
Qexc )
1508 #ifdef debugQGSParticipants
1509 G4cout<<
"Simulation of interaction with quark exchange."<<
G4endl;
1513 #ifdef debugQGSParticipants
1522 #ifdef debugQGSParticipants
1537 const G4double aHugeValue = 1.0e+10;
1539 #ifdef debugQGSParticipants
1549 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(350.), VqqM_tr(700);
1552 #ifdef debugQGSParticipants
1554 <<
"Target nucleon momenta at start"<<
G4endl;
1557 std::vector<G4VSplitableHadron*>::iterator i;
1562 Psum += (*i)->Get4Momentum();
1563 #ifdef debugQGSParticipants
1564 G4cout<<
"Nusleus nucleon # and its 4Mom. "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1579 #ifdef debugQGSParticipants
1581 G4cout<<
"Projectile 4 Mom "<<Projectile4Momentum<<
G4endl;
1589 (*i)->Set4Momentum( tmp );
1590 #ifdef debugQGSParticipants
1591 G4cout<<
"Target nucleon # and 4Mom "<<
" "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1593 Target4Momentum += tmp;
1600 #ifdef debugQGSParticipants
1602 G4cout<<
"Target nucleons mom: px, py, z_1, m_i"<<
G4endl;
1622 if ( Mass2 < 0.0 ) {
1625 <<
"Â 4-momentum " << Psum <<
G4endl;
1626 ed <<
"LorentzVector tmp " << tmp <<
" with Mass2 " << Mass2 <<
G4endl;
1627 G4Exception(
"G4QGSParticipants::DeterminePartonMomenta(): 4-momentum with negative mass!",
1630 Mass = std::sqrt( Mass2 );
1634 (*i)->Set4Momentum(tmp);
1635 #ifdef debugQGSParticipants
1636 G4cout<<
"Target nucleons # and mom: "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1649 G4double ProjSumMt(0.), ProjSumMt2perX(0.);
1650 G4double TargSumMt(0.), TargSumMt2perX(0.);
1666 const G4int maxNumberOfAttempts = 1000;
1669 attempt++;
if( attempt == 100*(attempt/100) ) {SigPt/=2.;}
1671 ProjSumMt=0.; ProjSumMt2perX=0.;
1672 TargSumMt=0.; TargSumMt2perX=0.;
1676 #ifdef debugQGSParticipants
1677 G4cout<<
"attempt ------------------------ "<<attempt<<
G4endl;
1684 G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1687 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1690 #ifdef debugQGSParticipants
1695 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1696 Mt = std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1700 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1703 NumberOfUnsampledSeaQuarks--;
1704 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1709 #ifdef debugQGSParticipants
1715 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1716 Mt = std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1720 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1723 NumberOfUnsampledSeaQuarks--;
1724 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1727 #ifdef debugQGSParticipants
1734 #ifdef debugQGSParticipants
1739 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1740 Mt = std::sqrt(aPtVector.
mag2()+
sqr(VqM_pr));
1744 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1747 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1753 #ifdef debugQGSParticipants
1759 Mt = std::sqrt(
sqr(SumPx) +
sqr(SumPy) +
sqr(VaqM_pr) );
1763 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1766 #ifdef debugQGSParticipants
1767 G4cout<<
" "<<tmp<<
" "<<SumZ+(1.-SumZ)<<
" (z-fraction)"<<
G4endl;
1777 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1778 #ifdef debugQGSParticipants
1780 <<
"Target nucleon 4Mom "<<(*i)->Get4Momentum()<<
G4endl;
1783 SumPx = (*i)->Get4Momentum().px() * (-1.);
1784 SumPy = (*i)->Get4Momentum().py() * (-1.);
1787 NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1790 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1792 aParton = (*i)->GetNextParton();
1793 #ifdef debugQGSParticipants
1798 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1799 Mt=std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1803 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1805 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1806 NumberOfUnsampledSeaQuarks--;
1807 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1811 aParton = (*i)->GetNextAntiParton();
1812 #ifdef debugQGSParticipants
1818 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1819 Mt=std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1823 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1825 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1826 NumberOfUnsampledSeaQuarks--;
1827 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1830 #ifdef debugQGSParticipants
1836 aParton = (*i)->GetNextParton();
1837 #ifdef debugQGSParticipants
1842 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1843 Mt=std::sqrt(aPtVector.
mag2()+
sqr(VqM_tr));
1847 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1849 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1850 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1855 aParton = (*i)->GetNextAntiParton();
1856 #ifdef debugQGSParticipants
1858 G4cout<<
" "<<tmp<<
" "<<SumZ<<
" (total z-sum) "<<
G4endl;
1862 Mt=std::sqrt(
sqr(SumPx) +
sqr(SumPy) +
sqr(VqqM_tr) );
1865 tmp.
setPz((*i)->Get4Momentum().pz()*(1.0 - SumZ));
1866 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1869 #ifdef debugQGSParticipants
1870 G4cout<<
" "<<tmp<<
" "<<1.0<<
" "<<(*i)->Get4Momentum().
pz()<<
G4endl;
1875 if( ProjSumMt + TargSumMt > SqrtS ) {
1876 Success =
false;
continue;}
1877 if( std::sqrt(ProjSumMt2perX) + std::sqrt(TargSumMt2perX) > SqrtS ) {
1878 Success =
false;
continue;}
1880 }
while( (!Success) &&
1881 attempt < maxNumberOfAttempts );
1883 if ( attempt >= maxNumberOfAttempts ) {
1890 - 2.0*
S*ProjSumMt2perX - 2.0*
S*TargSumMt2perX - 2.0*ProjSumMt2perX*TargSumMt2perX;
1892 G4double targetWminus=(
S - ProjSumMt2perX + TargSumMt2perX + std::sqrt( DecayMomentum2 ))/2.0/SqrtS;
1893 G4double projectileWplus = SqrtS - TargSumMt2perX/targetWminus;
1899 #ifdef debugQGSParticipants
1900 G4cout<<
"Backward transformation ===================="<<
G4endl;
1904 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1907 #ifdef debugQGSParticipants
1912 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1913 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1919 #ifdef debugQGSParticipants
1924 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1925 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1929 #ifdef debugQGSParticipants
1936 #ifdef debugQGSParticipants
1940 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1941 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1948 #ifdef debugQGSParticipants
1953 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1954 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1959 #ifdef debugQGSParticipants
1969 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1970 #ifdef debugQGSParticipants
1971 G4cout<<
"nSeaPair of target and N# "<<nSeaPair<<
" "<<NuclNo<<
G4endl;
1974 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1976 aParton = (*i)->GetNextParton();
1977 #ifdef debugQGSParticipants
1981 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1982 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1987 aParton = (*i)->GetNextAntiParton();
1988 #ifdef debugQGSParticipants
1993 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1994 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1998 #ifdef debugQGSParticipants
2005 aParton = (*i)->GetNextParton();
2006 #ifdef debugQGSParticipants
2010 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
2011 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
2017 aParton = (*i)->GetNextAntiParton();
2018 #ifdef debugQGSParticipants
2023 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
2024 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
2028 #ifdef debugQGSParticipants
2045 const G4int maxNumberOfLoops = 1000;
2046 G4int loopCounter = 0;
2052 }
while( ( r12 > 1.) &&
2053 ++loopCounter < maxNumberOfLoops );
2054 if ( loopCounter >= maxNumberOfLoops ) {
2066 #ifdef debugQGSParticipants
2071 #ifdef debugQGSParticipants
2072 G4cout<<
"BAD situation: theProjectileSplitable is NULL ! Returning immediately"<<
G4endl;
2077 #ifdef debugQGSParticipants
2089 #ifdef debugQGSParticipants
2106 for (
G4int i = 0; i < N_EnvTarg; i++ ) {
2111 #ifdef debugQGSParticipants
2121 #ifdef debugQGSParticipants
2139 #ifdef debugQGSParticipants
2140 G4cout<<
"Strings created in soft interactions"<<
G4endl;
2142 std::vector<G4InteractionContent*>::iterator i;
2150 #ifdef debugQGSParticipants
2151 G4cout<<
"An interaction # and soft coll. # "<<IntNo<<
" "
2164 #ifdef debugQGSParticipants
2179 #ifdef debugQGSParticipants
2187 #ifdef debugQGSParticipants
2203 #ifdef debugQGSParticipants
2213 #ifdef debugQGSParticipants
2217 #ifdef debugQGSParticipants
2228 #ifdef debugQGSParticipants
2252 residualMomentum +=tmp;
2274 const G4int maxNumberOfLoops = 1000;
2275 G4int loopCounter = 0;
2292 if(SumMasses > Mass) {Chigh=
C;}
2295 }
while( (Chigh-Clow > 0.01) &&
2296 ++loopCounter < maxNumberOfLoops );
2297 if ( loopCounter >= maxNumberOfLoops ) {
2298 #ifdef debugQGSParticipants
2319 #ifdef debugQGSParticipants
2320 G4cout <<
"End GetResiduals -----------------" <<
G4endl;
2329 std::vector<G4InteractionContent*>::iterator i;
2344 #ifdef debugQGSParticipants
2352 #ifdef debugQGSParticipants
2359 #ifdef debugQGSParticipants
2367 #ifdef debugQGSParticipants
2379 #ifdef debugQGSParticipants
2396 if ((!(aPrimaryMomentum.
e()>-1)) && (!(aPrimaryMomentum.
e()<1)) )
2399 "G4GammaParticipants::SelectInteractions: primary nan energy.");
2401 G4double S = (aPrimaryMomentum + aTargetNMomentum).mag2();
2405 #ifdef debugQGSParticipants
2406 G4cout <<
"Gamma projectile - SelectInteractions " <<
G4endl;
2431 {
if(NucleonNo == theCurrent)
break; NucleonNo++; }
2436 pNucleon->
Hit(aTarget);
G4double C(G4double temp)
G4double S(G4double temp)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static const G4double alpha
G4double G4Log(G4double x)
G4ThreadLocal G4int G4QGSParticipants_NPart
static constexpr double twopi
static constexpr double fermi
static constexpr double GeV
static constexpr double MeV
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
Hep3Vector boostVector() const
HepLorentzVector & boost(double, double, double)
void setVect(const Hep3Vector &)
Hep3Vector findBoostToCM() const
HepLorentzVector & transform(const HepRotation &)
static G4Gamma * GammaDefinition()
void SetNumberOfDiffractiveCollisions(int)
void SetTargetNucleon(G4Nucleon *aNucleon)
G4Nucleon * GetTargetNucleon() const
void SetNumberOfSoftCollisions(int)
G4VSplitableHadron * GetProjectile() const
G4int GetNumberOfSoftCollisions()
void SetTarget(G4VSplitableHadron *aTarget)
void SetStatus(G4int aValue)
G4VSplitableHadron * GetTarget() const
static G4KaonMinus * KaonMinusDefinition()
static G4KaonPlus * KaonPlusDefinition()
const G4ThreeVector & GetPosition() const
void SetPosition(const G4ThreeVector aPosition)
void SetMomentum(G4LorentzVector &aMomentum)
G4VSplitableHadron * GetSplitableHadron() const
virtual const G4LorentzVector & Get4Momentum() const
void Hit(G4VSplitableHadron *aHit)
G4double GetBindingEnergy() const
void SetBindingEnergy(G4double anEnergy)
virtual const G4ParticleDefinition * GetDefinition() const
G4double GetPDGMass() const
G4int GetPDGEncoding() const
G4double GetPDGCharge() const
G4int GetBaryonNumber() const
const G4String & GetParticleName() const
G4IonTable * GetIonTable() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
G4Parton * GetParton2(void)
G4Parton * GetParton1(void)
G4ParticleDefinition * GetDefinition()
void Set4Momentum(const G4LorentzVector &aMomentum)
const G4LorentzVector & Get4Momentum() const
static G4PionMinus * PionMinusDefinition()
static G4PionPlus * PionPlusDefinition()
static G4PionZero * PionZeroDefinition()
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction=TRUE) const
virtual G4Parton * GetNextParton()
virtual G4Parton * GetNextAntiParton()
const G4double ThresholdParameter
G4int TargetResidualMassNumber
G4V3DNucleus * GetTargetNucleus() const
G4double GetExcitationEnergyPerWoundedNucleon()
G4double GetMaxPt2ofNuclearDestruction()
G4LorentzVector ProjectileResidual4Momentum
G4double ProjectileResidualExcitationEnergy
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)
void SetR2ofNuclearDestruction(const G4double aValue)
G4double GetDofNuclearDestruction()
G4double SampleX(G4double anXmin, G4int nSea, G4int theTotalSea, G4double aBeta)
G4ThreeVector theCurrentVelocity
void SetDofNuclearDestruction(const G4double aValue)
G4QuarkExchange theQuarkExchange
G4SingleDiffractiveExcitation theSingleDiffExcitation
void SetExcitationEnergyPerWoundedNucleon(const G4double aValue)
virtual void DoLorentzBoost(G4ThreeVector aBoost)
std::vector< G4InteractionContent * > theInteractions
G4double GetR2ofNuclearDestruction()
G4bool DeterminePartonMomenta()
G4Nucleon * TheInvolvedNucleonsOfTarget[250]
void SetMaxPt2ofNuclearDestruction(const G4double aValue)
virtual G4VSplitableHadron * SelectInteractions(const G4ReactionProduct &thePrimary)
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)
G4bool FinalizeKinematics(const G4double w, const G4bool isProjectileNucleus, const G4LorentzRotation &boostFromCmsToLab, const G4double residualMass, const G4int residualMassNumber, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4LorentzVector &residual4Momentum)
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
void SetPt2ofNuclearDestruction(const G4double aValue)
G4double TargetResidualExcitationEnergy
void PerformDiffractiveCollisions()
G4bool GenerateDeltaIsobar(const G4double sqrtS, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4double &sumMasses)
G4QGSDiffractiveExcitation theDiffExcitaton
G4bool ComputeNucleusProperties(G4V3DNucleus *nucleus, G4LorentzVector &nucleusMomentum, G4LorentzVector &residualMomentum, G4double &sumMasses, G4double &residualExcitationEnergy, G4double &residualMass, G4int &residualMassNumber, G4int &residualCharge)
G4QGSMSplitableHadron * theProjectileSplitable
void StoreInvolvedNucleon()
G4Nucleon * TheInvolvedNucleonsOfProjectile[250]
void GetList(const G4ReactionProduct &thePrimary)
void PrepareInitialState(const G4ReactionProduct &thePrimary)
G4ReactionProduct theProjectile
virtual ~G4QGSParticipants()
void SetCofNuclearDestruction(const G4double aValue)
const G4double QGSMThreshold
G4LorentzVector TargetResidual4Momentum
G4V3DNucleus * GetProjectileNucleus() const
G4int NumberOfInvolvedNucleonsOfProjectile
const G4double theNucleonRadius
G4double GetPt2ofNuclearDestruction()
void BuildInteractions(const G4ReactionProduct &thePrimary)
G4int ProjectileResidualCharge
G4int TargetResidualCharge
std::vector< G4PartonPair * > thePartonPairs
std::vector< G4VSplitableHadron * > theTargets
G4int ProjectileResidualMassNumber
void PerformSoftCollisions()
G4int NumberOfInvolvedNucleonsOfTarget
G4double GetCofNuclearDestruction()
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
const G4ParticleDefinition * GetDefinition() const
G4double GetTotalEnergy() const
G4ThreeVector GetMomentum() const
void GetProbabilities(G4double B, G4int Mode, G4double &Pint, G4double &Pprd, G4double &Ptrd, G4double &Pdd, G4double &Pnd, G4double &Pnvr)
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction) const
virtual void SortNucleonsIncZ()=0
virtual G4double GetOuterRadius()=0
virtual G4Nucleon * GetNextNucleon()=0
virtual G4int GetCharge()=0
virtual G4bool StartLoop()=0
virtual void Init(G4int theA, G4int theZ, G4int numberOfLambdas=0)=0
std::pair< G4double, G4double > ChooseImpactXandY(G4double maxImpact)
virtual G4int GetMassNumber()=0
virtual void SetProjectileNucleus(G4V3DNucleus *aNucleus)
G4V3DNucleus * theNucleus
void SetTimeOfCreation(G4double aTime)
void SetStatus(const G4int aStatus)
void SetCollisionCount(G4int aCount)
const G4ParticleDefinition * GetDefinition() const
G4double GetTimeOfCreation()
void Set4Momentum(const G4LorentzVector &a4Momentum)
virtual G4Parton * GetNextParton()=0
virtual G4Parton * GetNextAntiParton()=0
void SetDefinition(const G4ParticleDefinition *aDefinition)
const G4LorentzVector & Get4Momentum() const
const G4ThreeVector & GetPosition() const
void IncrementCollisionCount(G4int aCount)
G4int GetSoftCollisionCount()