52 nCutMax(7),ThresholdParameter(0.000*
GeV),
69 ThresholdParameter(right.ThresholdParameter), QGSMThreshold(right.QGSMThreshold),
70 theNucleonRadius(right.theNucleonRadius)
113 #ifdef debugQGSParticipants
123 const G4int maxNumberOfLoops = 1000;
124 G4int loopCounter = 0;
127 const G4int maxNumberOfInternalLoops = 1000;
128 G4int internalLoopCounter = 0;
156 }
while( (!Success) && ++internalLoopCounter < maxNumberOfInternalLoops );
158 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
163 #ifdef debugQGSParticipants
164 G4cout<<G4endl<<
"PerformDiffractiveCollisions(); if they happend." <<
G4endl;
169 #ifdef debugQGSParticipants
176 #ifdef debugQGSParticipants
177 G4cout<<
"Perform non-Diffractive Collisions if they happend. Determine Parton Momenta and so on." <<
G4endl;
184 }
while( (!Success) && ++loopCounter < maxNumberOfLoops );
186 if ( loopCounter >= maxNumberOfLoops ) {
188 #ifdef debugQGSParticipants
198 #ifdef debugQGSParticipants
204 #ifdef debugQGSParticipants
215 #ifdef debugQGSParticipants
222 if ( (aNucleon != 0 ) && (aNucleon->
GetStatus() >= 1) )
delete aNucleon;
230 if ( aNucleon )
delete aNucleon;
234 #ifdef debugQGSParticipants
235 G4cout<<
"Delition of target nucleons from soft interactions "<<
theTargets.size()
254 if( pProjectile )
delete pProjectile;
266 if ( (splaNucleon != 0) && (splaNucleon->
GetStatus() >=1) )
delete splaNucleon;
267 aNucleon->
Hit(
nullptr);
298 #ifdef debugQGSParticipants
309 G4double SS=(aPrimaryMomentum + aNucleonMomentum).mag2();
318 #ifdef debugQGSParticipants
319 G4cout <<
"QGSM - BAD situation: pNucleon is NULL ! Leaving immediately!" <<
G4endl;
332 const G4int maxNumberOfLoops = 1000;
333 G4int loopCounter = -1;
334 while( (
theInteractions.size() == 0) && ++loopCounter < maxNumberOfLoops )
339 std::pair<G4double, G4double> theImpactParameter;
342 G4double impactX = theImpactParameter.first;
343 G4double impactY = theImpactParameter.second;
345 #ifdef debugQGSParticipants
352 G4int nucleonCount = -1;
361 if(Power <= 0.)
break;
370 G4double Pprd(0.), Ptrd(0.), Pdd(0.);
372 G4int NcutPomerons(0);
375 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr);
376 #ifdef debugQGSParticipants
377 G4cout<<
"Nucleon & its impact parameter: "<<nucleonCount<<
" "<<std::sqrt(Distance2)/
fermi<<
" (fm)"<<
G4endl;
379 G4cout<<
"Probability of PrD, TrD, DD: "<<Pprd<<
" "<<Ptrd<<
" "<<Pdd<<
G4endl;
380 G4cout<<
"Probability of NonDiff, QuarkExc.: "<<Pnd<<
" "<<Pnvr<<
" in inel. inter."<<
G4endl;
392 G4int InteractionType(0);
408 if( rndNumber < Ptrd ) {InteractionType =
TrD; }
409 else if( rndNumber < Ptrd + Pnd) {InteractionType =
NonD; NcutPomerons =
Regge->
ncPomerons();}
412 if( (InteractionType ==
NonD) && (NcutPomerons == 0))
continue;
416 tNucleon->
Hit(aTargetSPB);
418 #ifdef debugQGSParticipants
420 G4cout<<
"Target nucleon - "<<nucleonCount<<
" "
422 G4cout<<
"Interaction type:"<<InteractionType
423 <<
" (0 -PrD, 1 - TrD, 2 - DD, 3 - NonD, 4 - Qexc)"<<
G4endl;
425 <<
" (0 -ALL, 1 - WITHOUT_R, 2 - NON_DIFF)"<<
G4endl;
426 if( InteractionType ==
NonD )
427 G4cout<<
"Number of cutted pomerons: "<<NcutPomerons<<
G4endl;
430 if((InteractionType ==
PrD) || (InteractionType ==
TrD) || (InteractionType ==
DD) ||
431 (InteractionType ==
Qexc))
433 #ifdef debugQGSParticipants
447 aInteraction->
SetStatus(InteractionType);
452 #ifdef debugQGSParticipants
453 G4cout<<
"Non-diffractive interaction occurs, max NcutPomerons "<<NcutPomerons<<
G4endl;
459 for(nCuts = 0; nCuts < NcutPomerons; nCuts++)
461 if(
G4UniformRand() < Power/MaxPower ){Vncut++; Power--;
if(Power <= 0.)
break;}
466 if( nCuts == 0 ) {
delete aTargetSPB; tNucleon->
Hit(
nullptr);
continue;}
469 #ifdef debugQGSParticipants
470 G4cout<<
"Number of cuts in the interaction "<<nCuts<<
G4endl;
485 aInteraction->
SetStatus(InteractionType);
491 #ifdef debugQGSParticipants
497 if ( loopCounter >= maxNumberOfLoops ) {
498 #ifdef debugQGSParticipants
507 std::vector<G4InteractionContent*>::iterator i;
534 aTargetNucleon->
Hit(
nullptr);
548 #ifdef debugQGSParticipants
570 #ifdef debugQGSParticipants
572 <<
"Stored # of wounded nucleons of target "
582 #ifdef debugQGSParticipants
591 for (
G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
614 #ifdef debugQGSParticipants
615 G4cout<<
"Target nucleon involved in reggeon cascading No "<<TrgNuc<<
" "
623 Neighbour->
Hit( targetSplitable );
629 anInteraction->
SetTarget(targetSplitable);
641 #ifdef debugQGSParticipants
651 G4bool isProjectileNucleus =
false;
653 isProjectileNucleus =
true;
656 #ifdef debugPutOnMassShell
658 if ( isProjectileNucleus ) {
G4cout <<
"PutOnMassShell for Nucleus_Nucleus " <<
G4endl;}
662 if ( Pprojectile.z() < 0.0 ) {
674 #ifdef debugPutOnMassShell
682 if ( ! isOk )
return false;
691 if ( ! isProjectileNucleus ) {
692 Mprojectile = Pprojectile.mag();
693 M2projectile = Pprojectile.mag2();
694 SumMasses += Mprojectile + 20.0*
MeV;
697 #ifdef debugPutOnMassShell
698 G4cout <<
"Projectile : ";
704 if ( ! isOk )
return false;
711 #ifdef debugPutOnMassShell
714 G4cout <<
"Psum " << Psum/
GeV <<
" GeV" << G4endl <<
"SqrtS " << SqrtS/
GeV <<
" GeV" << G4endl
715 <<
"SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/
GeV <<
" "
716 << PrResidualMass/
GeV <<
" " << TargetResidualMass/
GeV <<
" GeV" <<
G4endl;
720 if ( SqrtS < SumMasses ) {
727 G4double savedSumMasses = SumMasses;
728 if ( isProjectileNucleus ) {
729 SumMasses -= std::sqrt(
sqr( PrResidualMass ) + PprojResidual.
perp2() );
731 + PprojResidual.
perp2() );
733 SumMasses -= std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.
perp2() );
735 + PtargetResidual.
perp2() );
737 if ( SqrtS < SumMasses ) {
738 SumMasses = savedSumMasses;
739 if ( isProjectileNucleus ) {
746 if ( isProjectileNucleus ) {
750 #ifdef debugPutOnMassShell
751 if ( isProjectileNucleus ) {
752 G4cout <<
"PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/
GeV <<
" "
755 G4cout <<
"TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/
GeV <<
" GeV "
757 <<
"Sum masses " << SumMasses/
GeV <<
G4endl;
761 if ( isProjectileNucleus && thePrNucleus->
GetMassNumber() != 1 ) {
770 if ( ! isOk )
return false;
783 if ( Ptmp.pz() <= 0.0 ) {
790 if ( isProjectileNucleus ) {
792 YprojectileNucleus = Ptmp.
rapidity();
794 Ptmp = toCms*Ptarget;
795 G4double YtargetNucleus = Ptmp.rapidity();
799 if ( isProjectileNucleus ) {
806 #ifdef debugPutOnMassShell
807 if ( isProjectileNucleus ) {
808 G4cout <<
"Y projectileNucleus " << YprojectileNucleus <<
G4endl;
810 G4cout <<
"Y targetNucleus " << YtargetNucleus << G4endl
812 <<
" DcorP DcorT " << DcorP <<
" " << DcorT <<
" AveragePt2 " << AveragePt2 <<
G4endl;
819 G4int NumberOfTries = 0;
821 G4bool OuterSuccess =
true;
823 const G4int maxNumberOfLoops = 1000;
824 G4int loopCounter = 0;
826 G4double sqrtM2proj = 0.0, sqrtM2target = 0.0;
828 const G4int maxNumberOfTries = 1000;
831 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
837 DcorP *= ScaleFactor;
838 DcorT *= ScaleFactor;
839 AveragePt2 *= ScaleFactor;
841 if ( isProjectileNucleus ) {
844 thePrNucleus, PprojResidual,
852 theTargetNucleus, PtargetResidual,
857 if ( M2proj < 0.0 ) {
861 <<
") M2proj=" << M2proj <<
" -> sets it to 0.0 !" <<
G4endl;
862 G4Exception(
"G4QGSParticipants::PutOnMassShell(): negative projectile squared mass!",
866 sqrtM2proj = std::sqrt( M2proj );
867 if ( M2target < 0.0 ) {
871 <<
") M2target=" << M2target <<
" -> sets it to 0.0 !" <<
G4endl;
872 G4Exception(
"G4QGSParticipants::PutOnMassShell(): negative target squared mass!",
876 sqrtM2target = std::sqrt( M2target );
878 #ifdef debugPutOnMassShell
879 G4cout <<
"SqrtS, Mp+Mt, Mp, Mt " << SqrtS/
GeV <<
" "
880 << ( sqrtM2proj + sqrtM2target )/
GeV <<
" "
881 << sqrtM2proj/
GeV <<
" " << sqrtM2target/
GeV << G4endl;
884 if ( ! isOk )
return false;
885 }
while ( ( SqrtS < ( sqrtM2proj + sqrtM2target ) ) &&
886 ++NumberOfTries < maxNumberOfTries );
887 if ( NumberOfTries >= maxNumberOfTries ) {
890 if ( isProjectileNucleus ) {
891 isOk =
CheckKinematics( S, SqrtS, M2proj, M2target, YprojectileNucleus,
true,
894 WminusTarget, WplusProjectile, OuterSuccess );
899 WminusTarget, WplusProjectile, OuterSuccess );
900 if ( ! isOk )
return false;
901 }
while ( ( ! OuterSuccess ) &&
902 ++loopCounter < maxNumberOfLoops );
903 if ( loopCounter >= maxNumberOfLoops ) {
913 if ( ! isProjectileNucleus ) {
915 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
916 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
917 Pprojectile.setPz( Pzprojectile );
918 Pprojectile.setE( Eprojectile );
920 #ifdef debugPutOnMassShell
924 Pprojectile.transform( toLab );
930 #ifdef debugPutOnMassShell
941 #ifdef debugPutOnMassShell
945 if ( ! isOk )
return false;
949 #ifdef debugPutOnMassShell
959 #ifdef debugPutOnMassShell
963 if ( ! isOk )
return false;
967 #ifdef debugPutOnMassShell
981 if ( AveragePt2 <= 0.0 ) {
985 (
G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
990 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
1001 G4int& residualMassNumber,
1002 G4int& residualCharge ) {
1014 if ( ! nucleus )
return false;
1039 sumMasses += 20.0*
MeV;
1042 residualExcitationEnergy += -ExcitationEPerWoundedNucleon*
1044 residualMassNumber--;
1051 #ifdef debugPutOnMassShell
1052 G4cout <<
"ExcitationEnergyPerWoundedNucleon " << ExcitationEPerWoundedNucleon <<
" MeV"<<
G4endl
1053 <<
"\t Residual Charge, MassNumber " << residualCharge <<
" " << residualMassNumber
1054 <<
G4endl <<
"\t Initial Momentum " << nucleusMomentum/
GeV<<
" GeV"
1055 <<
G4endl <<
"\t Residual Momentum " << residualMomentum/
GeV<<
" GeV"<<
G4endl;
1058 residualMomentum.
setPz( 0.0 );
1059 residualMomentum.
setE( 0.0 );
1060 if ( residualMassNumber == 0 ) {
1062 residualExcitationEnergy = 0.0;
1065 GetIonMass( residualCharge, residualMassNumber );
1066 if ( residualMassNumber == 1 ) {
1067 residualExcitationEnergy = 0.0;
1070 sumMasses += std::sqrt(
sqr( residualMass ) + residualMomentum.
perp2() );
1079 const G4int numberOfInvolvedNucleons,
1097 if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 )
return false;
1101 const G4double probDeltaIsobar = 0.10;
1103 G4int maxNumberOfDeltas =
G4int( (sqrtS - sumMasses)/(400.0*
MeV) );
1104 G4int numberOfDeltas = 0;
1106 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1109 if (
G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
1111 if ( ! involvedNucleons[i] )
continue;
1118 G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4;
1127 if ( sqrtS < sumMasses + massDelta - massNuc ) {
1131 sumMasses += ( massDelta - massNuc );
1150 const G4int residualMassNumber,
1151 const G4int numberOfInvolvedNucleons,
1165 if ( ! nucleus )
return false;
1167 if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
1175 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1176 G4Nucleon* aNucleon = involvedNucleons[i];
1177 if ( ! aNucleon )
continue;
1181 const G4int maxNumberOfLoops = 1000;
1182 G4int loopCounter = 0;
1189 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1190 G4Nucleon* aNucleon = involvedNucleons[i];
1191 if ( ! aNucleon )
continue;
1197 if ( x < 0.0 || x > 1.0 ) {
1207 if ( xSum < 0.0 || xSum > 1.0 ) success =
false;
1209 if ( ! success )
continue;
1211 G4double deltaPx = ( ptSum.
x() - pResidual.
x() ) / numberOfInvolvedNucleons;
1212 G4double deltaPy = ( ptSum.
y() - pResidual.
y() ) / numberOfInvolvedNucleons;
1214 if ( residualMassNumber == 0 ) {
1215 delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
1222 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1223 G4Nucleon* aNucleon = involvedNucleons[i];
1224 if ( ! aNucleon )
continue;
1227 if ( residualMassNumber == 0 ) {
1228 if ( x <= 0.0 || x > 1.0 ) {
1233 if ( x <= 0.0 || x > 1.0 || xSum <= 0.0 || xSum > 1.0 ) {
1241 +
sqr( px ) +
sqr( py ) ) / x;
1246 if ( success && residualMassNumber != 0 ) {
1247 mass2 += (
sqr( residualMass ) + pResidual.
perp2() ) / xSum;
1250 #ifdef debugPutOnMassShell
1254 }
while ( ( ! success ) &&
1255 ++loopCounter < maxNumberOfLoops );
1256 if ( loopCounter >= maxNumberOfLoops ) {
1272 const G4bool isProjectileNucleus,
1273 const G4int numberOfInvolvedNucleons,
1288 G4double decayMomentum2 =
sqr( sValue ) +
sqr( projectileMass2 ) +
sqr( targetMass2 )
1289 - 2.0*sValue*projectileMass2 - 2.0*sValue*targetMass2
1290 - 2.0*projectileMass2*targetMass2;
1291 targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
1293 projectileWplus = sqrtS - targetMass2/targetWminus;
1294 G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
1295 G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
1297 if(projectileE - projectilePz > 0.) {
1298 projectileY = 0.5 *
G4Log( (projectileE + projectilePz)/
1299 (projectileE - projectilePz) );
1301 G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
1302 G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
1303 G4double targetY = 0.5 *
G4Log( (targetE + targetPz)/(targetE - targetPz) );
1305 #ifdef debugPutOnMassShell
1306 G4cout <<
"decayMomentum2 " << decayMomentum2 <<
G4endl
1307 <<
"\t targetWminus projectileWplus " << targetWminus <<
" " << projectileWplus <<
G4endl
1308 <<
"\t projectileY targetY " << projectileY <<
" " << targetY <<
G4endl;
1311 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1312 G4Nucleon* aNucleon = involvedNucleons[i];
1313 if ( ! aNucleon )
continue;
1318 G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
1319 G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
1320 if ( isProjectileNucleus ) {
1321 pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*
x);
1322 e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*
x);
1326 #ifdef debugPutOnMassShell
1327 G4cout <<
"i nY pY nY-AY AY " << i <<
" " << nucleonY <<
" " << projectileY <<
G4endl;
1330 if ( std::abs( nucleonY - nucleusY ) > 2 ||
1331 ( isProjectileNucleus && targetY > nucleonY ) ||
1332 ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
1345 const G4bool isProjectileNucleus,
1348 const G4int residualMassNumber,
1349 const G4int numberOfInvolvedNucleons,
1367 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1368 G4Nucleon* aNucleon = involvedNucleons[i];
1369 if ( ! aNucleon )
continue;
1371 residual3Momentum -= tmp.
vect();
1375 G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w *
x );
1376 G4double e = w * x / 2.0 + mt2 / ( 2.0 * w *
x );
1378 if ( isProjectileNucleus ) pz *= -1.0;
1385 #ifdef debugPutOnMassShell // Uzhi 14.05.2015
1386 G4cout <<
"Target involved nucleon No, name, 4Mom "
1391 G4double residualMt2 =
sqr( residualMass ) +
sqr( residual3Momentum.
x() )
1392 +
sqr( residual3Momentum.
y() );
1394 #ifdef debugPutOnMassShell
1395 G4cout <<
G4endl<<
"w residual3Momentum.z() " << w <<
" " << residual3Momentum.
z() <<
G4endl;
1400 if ( residualMassNumber != 0 ) {
1401 residualPz = -w * residual3Momentum.
z() / 2.0 +
1402 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
1403 residualE = w * residual3Momentum.
z() / 2.0 +
1404 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
1406 if ( isProjectileNucleus ) residualPz *= -1.0;
1409 residual4Momentum.
setPx( residual3Momentum.
x() );
1410 residual4Momentum.
setPy( residual3Momentum.
y() );
1411 residual4Momentum.
setPz( residualPz );
1412 residual4Momentum.
setE( residualE );
1420 #ifdef debugQGSParticipants
1429 #ifdef debugQGSParticipants
1430 G4cout<<
"Interaction # and its status "
1435 if( (InterStatus ==
PrD) || (InterStatus ==
TrD) || (InterStatus ==
DD))
1437 #ifdef debugQGSParticipants
1438 G4cout<<
"Simulation of diffractive interaction. "<<InterStatus
1439 <<
" PrD/TrD/DD/ND/Qech - 0,1,2,3,4"<<
G4endl;
1444 #ifdef debugQGSParticipants
1445 G4cout<<
"The proj. before inter "
1448 G4cout<<
"The targ. before inter "
1453 if( InterStatus ==
PrD )
1456 if( InterStatus ==
TrD )
1459 if( InterStatus ==
DD )
1462 #ifdef debugQGSParticipants
1463 G4cout<<
"The proj. after inter "
1466 G4cout<<
"The targ. after inter "
1472 if( InterStatus ==
Qexc )
1474 #ifdef debugQGSParticipants
1475 G4cout<<
"Simulation of interaction with quark exchange."<<
G4endl;
1479 #ifdef debugQGSParticipants
1480 G4cout<<
"The proj. before inter "
1483 G4cout<<
"The targ. before inter "
1490 #ifdef debugQGSParticipants
1491 G4cout<<
"The proj. after inter "
1494 G4cout<<
"The targ. after inter "
1507 const G4double aHugeValue = 1.0e+10;
1509 #ifdef debugQGSParticipants
1519 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(350.), VqqM_tr(700);
1523 #ifdef debugQGSParticipants
1524 G4cout<<
"Projectile 4 momentum "<<Psum<<G4endl
1525 <<
"Target nucleon momenta at start"<<
G4endl;
1528 std::vector<G4VSplitableHadron*>::iterator i;
1533 Psum += (*i)->Get4Momentum();
1534 #ifdef debugQGSParticipants
1535 G4cout<<
"Nusleus nucleon # and its 4Mom. "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1544 toCms.
rotateZ( -1*Ptmp.phi() );
1545 toCms.
rotateY( -1*Ptmp.theta() );
1550 #ifdef debugQGSParticipants
1552 G4cout<<
"Projectile 4 Mom "<<Projectile4Momentum<<
G4endl;
1560 (*i)->Set4Momentum( tmp );
1561 #ifdef debugQGSParticipants
1562 G4cout<<
"Target nucleon # and 4Mom "<<
" "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1564 Target4Momentum +=
tmp;
1571 #ifdef debugQGSParticipants
1572 G4cout<<
"Sum of target nucleons 4 momentum "<<Target4Momentum<<G4endl<<
G4endl;
1573 G4cout<<
"Target nucleons mom: px, py, z_1, m_i"<<
G4endl;
1593 if ( Mass2 < 0.0 ) {
1596 <<
" 4-momentum " << Psum <<
G4endl;
1597 ed <<
"LorentzVector tmp " << tmp <<
" with Mass2 " << Mass2 <<
G4endl;
1598 G4Exception(
"G4QGSParticipants::DeterminePartonMomenta(): 4-momentum with negative mass!",
1601 Mass = std::sqrt( Mass2 );
1605 (*i)->Set4Momentum(tmp);
1606 #ifdef debugQGSParticipants
1607 G4cout<<
"Target nucleons # and mom: "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1621 G4double ProjSumMt(0.), ProjSumMt2perX(0.);
1622 G4double TargSumMt(0.), TargSumMt2perX(0.);
1638 const G4int maxNumberOfAttempts = 1000;
1641 attempt++;
if( attempt == 100*(attempt/100) ) {SigPt/=2.;}
1644 ProjSumMt=0.; ProjSumMt2perX=0.;
1645 TargSumMt=0.; TargSumMt2perX=0.;
1649 #ifdef debugQGSParticipants
1650 G4cout<<
"attempt ------------------------ "<<attempt<<
G4endl;
1657 G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1660 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1663 #ifdef debugQGSParticipants
1668 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1669 Mt = std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1674 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1677 NumberOfUnsampledSeaQuarks--;
1678 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1683 #ifdef debugQGSParticipants
1689 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1690 Mt = std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1695 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1698 NumberOfUnsampledSeaQuarks--;
1699 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1702 #ifdef debugQGSParticipants
1709 #ifdef debugQGSParticipants
1714 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1715 Mt = std::sqrt(aPtVector.
mag2()+
sqr(VqM_pr));
1720 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1723 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1729 #ifdef debugQGSParticipants
1734 Mt = std::sqrt(aPtVector.
mag2()+
sqr(VaqM_pr));
1738 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1741 #ifdef debugQGSParticipants
1742 G4cout<<
" "<<tmp<<
" "<<SumZ+(1.-SumZ)<<
" (z-fraction)"<<G4endl;
1751 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1752 #ifdef debugQGSParticipants
1753 G4cout<<
"nSeaPair of target N "<<nSeaPair<<G4endl
1754 <<
"Target nucleon 4Mom "<<(*i)->Get4Momentum()<<
G4endl;
1757 SumPx = (*i)->Get4Momentum().px() * (-1.);
1758 SumPy = (*i)->Get4Momentum().py() * (-1.);
1762 NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1765 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1767 aParton = (*i)->GetNextParton();
1768 #ifdef debugQGSParticipants
1773 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1774 Mt=std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1779 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1781 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1783 NumberOfUnsampledSeaQuarks--;
1784 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1788 aParton = (*i)->GetNextAntiParton();
1789 #ifdef debugQGSParticipants
1795 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1796 Mt=std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1801 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1803 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1805 NumberOfUnsampledSeaQuarks--;
1806 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1809 #ifdef debugQGSParticipants
1815 aParton = (*i)->GetNextParton();
1816 #ifdef debugQGSParticipants
1821 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1822 Mt=std::sqrt(aPtVector.
mag2()+
sqr(VqM_tr));
1827 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1829 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1831 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1836 aParton = (*i)->GetNextAntiParton();
1837 #ifdef debugQGSParticipants
1839 G4cout<<
" "<<tmp<<
" "<<SumZw<<
" (sum z-fracs) "<<SumZ<<
" (total z-sum) "<<
G4endl;
1842 Mt=std::sqrt(aPtVector.
mag2()+
sqr(VqqM_tr));
1845 tmp.
setPz((*i)->Get4Momentum().pz()*(1.0 - SumZ));
1847 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1850 #ifdef debugQGSParticipants
1851 G4cout<<
" "<<tmp<<
" "<<SumZw<<
" "<<1.0<<
" "<<(*i)->Get4Momentum().
pz()<<
G4endl;
1856 if( ProjSumMt + TargSumMt > SqrtS ) {
1857 Success =
false;
continue;}
1858 if( std::sqrt(ProjSumMt2perX) + std::sqrt(TargSumMt2perX) > SqrtS ) {
1859 Success =
false;
continue;}
1861 }
while( (!Success) &&
1862 attempt < maxNumberOfAttempts );
1864 if ( attempt >= maxNumberOfAttempts ) {
1870 - 2.0*S*ProjSumMt2perX - 2.0*S*TargSumMt2perX - 2.0*ProjSumMt2perX*TargSumMt2perX;
1872 G4double targetWminus=( S - ProjSumMt2perX + TargSumMt2perX + std::sqrt( DecayMomentum2 ))/2.0/SqrtS;
1873 G4double projectileWplus = SqrtS - TargSumMt2perX/targetWminus;
1879 #ifdef debugQGSParticipants
1880 G4cout<<
"Backward transformation ===================="<<
G4endl;
1884 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1887 #ifdef debugQGSParticipants
1892 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1893 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1899 #ifdef debugQGSParticipants
1904 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1905 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1909 #ifdef debugQGSParticipants
1916 #ifdef debugQGSParticipants
1920 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1921 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1928 #ifdef debugQGSParticipants
1933 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1934 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1939 #ifdef debugQGSParticipants
1948 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1949 #ifdef debugQGSParticipants
1950 G4cout<<
"nSeaPair of target and N# "<<nSeaPair<<
" "<<NuclNo<<
G4endl;
1953 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1955 aParton = (*i)->GetNextParton();
1956 #ifdef debugQGSParticipants
1960 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1961 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1966 aParton = (*i)->GetNextAntiParton();
1967 #ifdef debugQGSParticipants
1972 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1973 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1977 #ifdef debugQGSParticipants
1984 aParton = (*i)->GetNextParton();
1985 #ifdef debugQGSParticipants
1989 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1990 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1996 aParton = (*i)->GetNextAntiParton();
1997 #ifdef debugQGSParticipants
2002 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
2003 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
2007 #ifdef debugQGSParticipants
2026 const G4int maxNumberOfLoops = 1000;
2027 G4int loopCounter = 0;
2033 }
while( ( r12 > 1.) &&
2034 ++loopCounter < maxNumberOfLoops );
2035 if ( loopCounter >= maxNumberOfLoops ) {
2047 #ifdef debugQGSParticipants
2052 #ifdef debugQGSParticipants
2053 G4cout<<
"BAD situation: theProjectileSplitable is NULL ! Returning immediately"<<
G4endl;
2058 #ifdef debugQGSParticipants
2071 #ifdef debugQGSParticipants
2089 for (
G4int i = 0; i < N_EnvTarg; i++ ) {
2094 #ifdef debugQGSParticipants
2105 #ifdef debugQGSParticipants
2123 #ifdef debugQGSParticipants
2124 G4cout<<
"Strings created in soft interactions"<<
G4endl;
2126 std::vector<G4InteractionContent*>::iterator i;
2134 #ifdef debugQGSParticipants
2135 G4cout<<
"An interaction # and soft coll. # "<<IntNo<<
" "
2148 #ifdef debugQGSParticipants
2166 #ifdef debugQGSParticipants
2176 #ifdef debugQGSParticipants
2191 #ifdef debugQGSParticipants
2192 G4cout <<
"Sum of strings 4 momenta " << str4Mom << G4endl<<
G4endl;
2201 #ifdef debugQGSParticipants
2202 G4cout <<
"GetResiduals(): GetProjectileNucleus()? "
2206 #ifdef debugQGSParticipants
2218 #ifdef debugQGSParticipants
2242 residualMomentum +=
tmp;
2264 const G4int maxNumberOfLoops = 1000;
2265 G4int loopCounter = 0;
2282 if(SumMasses > Mass) {Chigh=
C;}
2285 }
while( (Chigh-Clow > 0.01) &&
2286 ++loopCounter < maxNumberOfLoops );
2287 if ( loopCounter >= maxNumberOfLoops ) {
2288 #ifdef debugQGSParticipants
2310 #ifdef debugQGSParticipants
2311 G4cout <<
"End GetResiduals -----------------" <<
G4endl;
2320 std::vector<G4InteractionContent*>::iterator i;
2335 #ifdef debugQGSParticipants
2343 #ifdef debugQGSParticipants
2350 #ifdef debugQGSParticipants
2358 #ifdef debugQGSParticipants
2370 #ifdef debugQGSParticipants
2387 if((!(aPrimaryMomentum.e()>-1)) && (!(aPrimaryMomentum.e()<1)) )
2390 "G4GammaParticipants::SelectInteractions: primary nan energy.");
2392 G4double S = (aPrimaryMomentum + aTargetNMomentum).mag2();
2396 #ifdef debugQGSParticipants
2397 G4cout <<
"Gamma projectile - SelectInteractions " <<
G4endl;
2399 G4cout <<
"SqrtS ThresholdMass ModelMode " <<std::sqrt(S)<<
" "<<ThresholdMass<<
" "<<
ModelMode<<
G4endl;
2414 G4int totalCuts = 0;
2428 {
if(NucleonNo == theCurrent)
break; NucleonNo++;}
2433 pNucleon->
Hit(aTarget);
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
void SetStatus(const G4int aStatus)
G4bool CheckKinematics(const G4double sValue, const G4double sqrtS, const G4double projectileMass2, const G4double targetMass2, const G4double nucleusY, const G4bool isProjectileNucleus, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4double &targetWminus, G4double &projectileWplus, G4bool &success)
G4double GetMaxPt2ofNuclearDestruction()
G4double GetExcitationEnergyPerWoundedNucleon()
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4QGSMSplitableHadron * theProjectileSplitable
void SetDofNuclearDestruction(const G4double aValue)
const G4ParticleDefinition * GetDefinition() const
HepLorentzVector & transform(const HepRotation &)
void SetStatus(G4int aValue)
CLHEP::Hep3Vector G4ThreeVector
G4Parton * GetParton2(void)
void SetExcitationEnergyPerWoundedNucleon(const G4double aValue)
std::ostringstream G4ExceptionDescription
std::pair< G4double, G4double > ChooseImpactXandY(G4double maxImpact)
G4bool FinalizeKinematics(const G4double w, const G4bool isProjectileNucleus, const G4LorentzRotation &boostFromCmsToLab, const G4double residualMass, const G4int residualMassNumber, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4LorentzVector &residual4Momentum)
static G4ParticleTable * GetParticleTable()
G4ParticleDefinition * GetDefinition()
static constexpr double MeV
void SetMomentum(G4LorentzVector &aMomentum)
void SetBindingEnergy(G4double anEnergy)
G4int NumberOfInvolvedNucleonsOfProjectile
const G4double QGSMThreshold
std::vector< G4InteractionContent * > theInteractions
G4Nucleon * GetTargetNucleon() const
G4double GetDofNuclearDestruction()
G4VSplitableHadron * GetTarget() const
void SetMomentum(const G4double x, const G4double y, const G4double z)
Hep3Vector findBoostToCM() const
const G4String & GetParticleName() const
G4IonTable * GetIonTable() const
G4double GetPDGCharge() const
void setVect(const Hep3Vector &)
void SetPt2ofNuclearDestruction(const G4double aValue)
const G4ParticleDefinition * GetDefinition() const
virtual G4Parton * GetNextParton()=0
static G4PionPlus * PionPlusDefinition()
G4double GetPDGMass() const
G4double G4Log(G4double x)
G4VSplitableHadron * GetSplitableHadron() const
virtual void SetProjectileNucleus(G4V3DNucleus *aNucleus)
void SetTotalEnergy(const G4double en)
G4VSplitableHadron * GetProjectile() const
G4LorentzVector ProjectileResidual4Momentum
virtual ~G4QGSParticipants()
void BuildInteractions(const G4ReactionProduct &thePrimary)
const G4ThreeVector & GetPosition() const
G4double GetTotalEnergy() const
virtual G4VSplitableHadron * SelectInteractions(const G4ReactionProduct &thePrimary)
G4LorentzVector TargetResidual4Momentum
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
static constexpr double fermi
G4int GetNumberOfSoftCollisions()
void SetTimeOfCreation(G4double aTime)
virtual G4Parton * GetNextAntiParton()
G4V3DNucleus * GetProjectileNucleus() const
void SetNumberOfDiffractiveCollisions(int)
static G4Pow * GetInstance()
G4Nucleon * TheInvolvedNucleonsOfTarget[250]
G4double GetBindingEnergy() const
virtual void Init(G4int theA, G4int theZ)=0
void SetNumberOfSoftCollisions(int)
const G4ThreeVector & GetPosition() const
G4double powA(G4double A, G4double y) const
G4int GetBaryonNumber() const
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction) const
void PerformDiffractiveCollisions()
G4Parton * GetParton1(void)
G4V3DNucleus * GetTargetNucleus() const
G4int ProjectileResidualCharge
static const G4double alpha
G4bool GenerateDeltaIsobar(const G4double sqrtS, const G4int numberOfInvolvedNucleons, G4Nucleon *involvedNucleons[], G4double &sumMasses)
static constexpr double twopi
static G4PionZero * PionZeroDefinition()
void GetList(const G4ReactionProduct &thePrimary)
G4bool ComputeNucleusProperties(G4V3DNucleus *nucleus, G4LorentzVector &nucleusMomentum, G4LorentzVector &residualMomentum, G4double &sumMasses, G4double &residualExcitationEnergy, G4double &residualMass, G4int &residualMassNumber, G4int &residualCharge)
HepLorentzVector & rotateY(double)
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
virtual G4Parton * GetNextAntiParton()=0
G4ThreadLocal G4int G4QGSParticipants_NPart
G4ThreeVector GetMomentum() const
virtual G4int GetCharge()=0
virtual G4Parton * GetNextParton()
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction=TRUE) const
virtual void DoLorentzBoost(G4ThreeVector aBoost)
static G4PionMinus * PionMinusDefinition()
virtual G4int GetMassNumber()=0
std::vector< G4VSplitableHadron * > theTargets
G4int ProjectileResidualMassNumber
virtual G4Nucleon * GetNextNucleon()=0
G4double G4ParticleHPJENDLHEData::G4double result
void StoreInvolvedNucleon()
virtual const G4LorentzVector & Get4Momentum() const
static G4KaonMinus * KaonMinusDefinition()
G4int GetSoftCollisionCount()
static G4Gamma * GammaDefinition()
virtual const G4ParticleDefinition * GetDefinition() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double SampleX(G4double anXmin, G4int nSea, G4int theTotalSea, G4double aBeta)
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)
HepLorentzVector & rotateZ(double)
G4double TargetResidualExcitationEnergy
G4ThreeVector theCurrentVelocity
void Hit(G4VSplitableHadron *aHit)
void SetTarget(G4VSplitableHadron *aTarget)
G4double GetCofNuclearDestruction()
G4SingleDiffractiveExcitation theSingleDiffExcitation
G4int NumberOfInvolvedNucleonsOfTarget
G4Nucleon * TheInvolvedNucleonsOfProjectile[250]
void SetDefinition(const G4ParticleDefinition *aDefinition)
G4double ProjectileResidualExcitationEnergy
void SetCofNuclearDestruction(const G4double aValue)
const G4double ThresholdParameter
G4double GetR2ofNuclearDestruction()
G4int GetPDGEncoding() const
G4GLOB_DLL std::ostream G4cout
void PerformSoftCollisions()
G4ReactionProduct theProjectile
Hep3Vector boostVector() const
const G4LorentzVector & Get4Momentum() const
void PrepareInitialState(const G4ReactionProduct &thePrimary)
G4QGSDiffractiveExcitation theDiffExcitaton
static G4KaonPlus * KaonPlusDefinition()
std::vector< G4PartonPair * > thePartonPairs
void SetCollisionCount(G4int aCount)
G4double GetPt2ofNuclearDestruction()
void SetMaxPt2ofNuclearDestruction(const G4double aValue)
void SetTargetNucleon(G4Nucleon *aNucleon)
const G4LorentzVector & Get4Momentum() const
const G4double theNucleonRadius
G4QuarkExchange theQuarkExchange
virtual G4double GetOuterRadius()=0
virtual G4bool StartLoop()=0
virtual void SortNucleonsIncZ()=0
static constexpr double GeV
void GetProbabilities(G4double B, G4int Mode, G4double &Pint, G4double &Pprd, G4double &Ptrd, G4double &Pdd, G4double &Pnd, G4double &Pnvr)
void SetR2ofNuclearDestruction(const G4double aValue)
HepLorentzRotation & transform(const HepBoost &b)
G4double GetTimeOfCreation()
G4int TargetResidualMassNumber
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4int TargetResidualCharge
void Set4Momentum(const G4LorentzVector &aMomentum)
void IncrementCollisionCount(G4int aCount)
HepLorentzVector & boost(double, double, double)
G4bool DeterminePartonMomenta()
G4V3DNucleus * theNucleus