51 nCutMax(7),ThresholdParameter(0.000*
GeV),
52 QGSMThreshold(3*
GeV),theNucleonRadius(1.5*
fermi),
alpha(-0.5),beta(2.5)
67 ThresholdParameter(right.ThresholdParameter), QGSMThreshold(right.QGSMThreshold),
68 theNucleonRadius(right.theNucleonRadius)
109 #ifdef debugQGSParticipants
119 const G4int maxNumberOfLoops = 1000;
120 G4int loopCounter = 0;
123 const G4int maxNumberOfInternalLoops = 1000;
124 G4int internalLoopCounter = 0;
146 }
while( (!Success) && ++internalLoopCounter < maxNumberOfInternalLoops );
148 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
153 #ifdef debugQGSParticipants
154 G4cout<<G4endl<<
"PerformDiffractiveCollisions(); if they happend." <<
G4endl;
159 #ifdef debugQGSParticipants
166 #ifdef debugQGSParticipants
167 G4cout<<
"Perform non-Diffractive Collisions if they happend. Determine Parton Momenta and so on." <<
G4endl;
174 }
while( (!Success) && ++loopCounter < maxNumberOfLoops );
176 if ( loopCounter >= maxNumberOfLoops ) {
178 #ifdef debugQGSParticipants
188 #ifdef debugQGSParticipants
194 #ifdef debugQGSParticipants
205 #ifdef debugQGSParticipants
212 if ( (aNucleon != 0 ) && (aNucleon->
GetStatus() >= 1) )
delete aNucleon;
220 if ( aNucleon )
delete aNucleon;
224 #ifdef debugQGSParticipants
225 G4cout<<
"Delition of target nucleons from soft interactions "<<
theTargets.size()
244 if( pProjectile )
delete pProjectile;
256 if ( (splaNucleon != 0) && (splaNucleon->
GetStatus() >=1) )
delete splaNucleon;
257 aNucleon->
Hit(
nullptr);
296 #ifdef debugQGSParticipants
307 G4double SS=(aPrimaryMomentum + aNucleonMomentum).mag2();
316 #ifdef debugQGSParticipants
317 G4cout <<
"QGSM - BAD situation: pNucleon is NULL ! Leaving immediately!" <<
G4endl;
328 theNucleusOuterR = 0.;
339 const G4int maxNumberOfLoops = 1000;
341 G4int NumberOfTries = 0;
344 G4int loopCounter = -1;
345 while( (
theInteractions.size() == 0) && ++loopCounter < maxNumberOfLoops )
350 std::pair<G4double, G4double> theImpactParameter;
353 if( NumberOfTries == 100*(NumberOfTries/100) ) Scale /=2.0;
356 G4double impactX = theImpactParameter.first;
357 G4double impactY = theImpactParameter.second;
359 #ifdef debugQGSParticipants
366 G4int nucleonCount = -1;
373 if(Power <= 0.)
break;
382 G4double Pprd(0.), Ptrd(0.), Pdd(0.);
384 G4int NcutPomerons(0);
387 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr);
388 #ifdef debugQGSParticipants
389 G4cout<<
"Nucleon & its impact parameter: "<<nucleonCount<<
" "<<std::sqrt(Distance2)/
fermi<<
" (fm)"<<
G4endl;
391 G4cout<<
"Probability of PrD, TrD, DD: "<<Pprd<<
" "<<Ptrd<<
" "<<Pdd<<
G4endl;
392 G4cout<<
"Probability of NonDiff, QuarkExc.: "<<Pnd<<
" "<<Pnvr<<
" in inel. inter."<<
G4endl;
399 G4int InteractionType(0);
413 if( rndNumber < Ptrd ) {InteractionType =
TrD; }
414 else if( rndNumber < Ptrd + Pnd) {InteractionType =
NonD; NcutPomerons =
Regge->
ncPomerons();}
417 if( (InteractionType ==
NonD) && (NcutPomerons == 0))
continue;
421 tNucleon->
Hit(aTargetSPB);
423 #ifdef debugQGSParticipants
425 G4cout<<
"Target nucleon - "<<nucleonCount<<
" "
427 G4cout<<
"Interaction type:"<<InteractionType
428 <<
" (0 -PrD, 1 - TrD, 2 - DD, 3 - NonD, 4 - Qexc)"<<
G4endl;
430 <<
" (0 -ALL, 1 - WITHOUT_R, 2 - NON_DIFF)"<<
G4endl;
431 if( InteractionType ==
NonD )
432 G4cout<<
"Number of cutted pomerons: "<<NcutPomerons<<
G4endl;
435 if((InteractionType ==
PrD) || (InteractionType ==
TrD) || (InteractionType ==
DD) ||
436 (InteractionType ==
Qexc))
438 #ifdef debugQGSParticipants
452 aInteraction->
SetStatus(InteractionType);
457 #ifdef debugQGSParticipants
458 G4cout<<
"Non-diffractive interaction occurs, max NcutPomerons "<<NcutPomerons<<
G4endl;
464 for(nCuts = 0; nCuts < NcutPomerons; nCuts++)
466 if(
G4UniformRand() < Power/MaxPower ){Vncut++; Power--;
if(Power <= 0.)
break;}
470 if( nCuts == 0 ) {
delete aTargetSPB; tNucleon->
Hit(
nullptr);
continue;}
473 #ifdef debugQGSParticipants
474 G4cout<<
"Number of cuts in the interaction "<<nCuts<<
G4endl;
488 aInteraction->
SetStatus(InteractionType);
494 #ifdef debugQGSParticipants
500 if ( loopCounter >= maxNumberOfLoops ) {
501 #ifdef debugQGSParticipants
510 std::vector<G4InteractionContent*>::iterator i;
538 aTargetNucleon->
Hit(
nullptr);
551 #ifdef debugQGSParticipants
573 #ifdef debugQGSParticipants
575 <<
"Stored # of wounded nucleons of target "
585 #ifdef debugQGSParticipants
594 for (
G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
617 #ifdef debugQGSParticipants
618 G4cout<<
"Target nucleon involved in reggeon cascading No "<<TrgNuc<<
" "
626 Neighbour->
Hit( targetSplitable );
632 anInteraction->
SetTarget(targetSplitable);
644 #ifdef debugQGSParticipants
654 G4bool isProjectileNucleus =
false;
656 isProjectileNucleus =
true;
659 #ifdef debugPutOnMassShell
661 if ( isProjectileNucleus ) {
G4cout <<
"PutOnMassShell for Nucleus_Nucleus " <<
G4endl;}
665 if ( Pprojectile.z() < 0.0 ) {
677 #ifdef debugPutOnMassShell
685 if ( ! isOk )
return false;
694 if ( ! isProjectileNucleus ) {
695 Mprojectile = Pprojectile.mag();
696 M2projectile = Pprojectile.mag2();
697 SumMasses += Mprojectile + 20.0*
MeV;
700 #ifdef debugPutOnMassShell
701 G4cout <<
"Projectile : ";
707 if ( ! isOk )
return false;
714 #ifdef debugPutOnMassShell
717 G4cout <<
"Psum " << Psum/
GeV <<
" GeV" << G4endl <<
"SqrtS " << SqrtS/
GeV <<
" GeV" << G4endl
718 <<
"SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/
GeV <<
" "
719 << PrResidualMass/
GeV <<
" " << TargetResidualMass/
GeV <<
" GeV" <<
G4endl;
723 if ( SqrtS < SumMasses ) {
730 G4double savedSumMasses = SumMasses;
731 if ( isProjectileNucleus ) {
732 SumMasses -= std::sqrt(
sqr( PrResidualMass ) + PprojResidual.
perp2() );
734 + PprojResidual.
perp2() );
736 SumMasses -= std::sqrt(
sqr( TargetResidualMass ) + PtargetResidual.
perp2() );
738 + PtargetResidual.
perp2() );
740 if ( SqrtS < SumMasses ) {
741 SumMasses = savedSumMasses;
742 if ( isProjectileNucleus ) {
749 if ( isProjectileNucleus ) {
753 #ifdef debugPutOnMassShell
754 if ( isProjectileNucleus ) {
755 G4cout <<
"PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/
GeV <<
" "
758 G4cout <<
"TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/
GeV <<
" GeV "
760 <<
"Sum masses " << SumMasses/
GeV <<
G4endl;
764 if ( isProjectileNucleus && thePrNucleus->
GetMassNumber() != 1 ) {
773 if ( ! isOk )
return false;
786 if ( Ptmp.pz() <= 0.0 ) {
793 if ( isProjectileNucleus ) {
795 YprojectileNucleus = Ptmp.
rapidity();
797 Ptmp = toCms*Ptarget;
798 G4double YtargetNucleus = Ptmp.rapidity();
802 if ( isProjectileNucleus ) {
809 #ifdef debugPutOnMassShell
810 if ( isProjectileNucleus ) {
811 G4cout <<
"Y projectileNucleus " << YprojectileNucleus <<
G4endl;
813 G4cout <<
"Y targetNucleus " << YtargetNucleus << G4endl
815 <<
" DcorP DcorT " << DcorP <<
" " << DcorT <<
" AveragePt2 " << AveragePt2 <<
G4endl;
822 G4int NumberOfTries = 0;
824 G4bool OuterSuccess =
true;
826 const G4int maxNumberOfLoops = 1000;
827 G4int loopCounter = 0;
829 G4double sqrtM2proj = 0.0, sqrtM2target = 0.0;
831 const G4int maxNumberOfTries = 1000;
834 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
840 DcorP *= ScaleFactor;
841 DcorT *= ScaleFactor;
842 AveragePt2 *= ScaleFactor;
844 if ( isProjectileNucleus ) {
847 thePrNucleus, PprojResidual,
855 theTargetNucleus, PtargetResidual,
860 if ( M2proj < 0.0 ) {
864 <<
") M2proj=" << M2proj <<
" -> sets it to 0.0 !" <<
G4endl;
865 G4Exception(
"G4QGSParticipants::PutOnMassShell(): negative projectile squared mass!",
869 sqrtM2proj = std::sqrt( M2proj );
870 if ( M2target < 0.0 ) {
874 <<
") M2target=" << M2target <<
" -> sets it to 0.0 !" <<
G4endl;
875 G4Exception(
"G4QGSParticipants::PutOnMassShell(): negative target squared mass!",
879 sqrtM2target = std::sqrt( M2target );
881 #ifdef debugPutOnMassShell
882 G4cout <<
"SqrtS, Mp+Mt, Mp, Mt " << SqrtS/
GeV <<
" "
883 << ( sqrtM2proj + sqrtM2target )/
GeV <<
" "
884 << sqrtM2proj/
GeV <<
" " << sqrtM2target/
GeV << G4endl;
887 if ( ! isOk )
return false;
888 }
while ( ( SqrtS < ( sqrtM2proj + sqrtM2target ) ) &&
889 ++NumberOfTries < maxNumberOfTries );
890 if ( NumberOfTries >= maxNumberOfTries ) {
893 if ( isProjectileNucleus ) {
894 isOk =
CheckKinematics( S, SqrtS, M2proj, M2target, YprojectileNucleus,
true,
897 WminusTarget, WplusProjectile, OuterSuccess );
902 WminusTarget, WplusProjectile, OuterSuccess );
903 if ( ! isOk )
return false;
904 }
while ( ( ! OuterSuccess ) &&
905 ++loopCounter < maxNumberOfLoops );
906 if ( loopCounter >= maxNumberOfLoops ) {
916 if ( ! isProjectileNucleus ) {
918 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
919 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
920 Pprojectile.setPz( Pzprojectile );
921 Pprojectile.setE( Eprojectile );
923 #ifdef debugPutOnMassShell
927 Pprojectile.transform( toLab );
933 #ifdef debugPutOnMassShell
944 #ifdef debugPutOnMassShell
948 if ( ! isOk )
return false;
952 #ifdef debugPutOnMassShell
962 #ifdef debugPutOnMassShell
966 if ( ! isOk )
return false;
970 #ifdef debugPutOnMassShell
984 if ( AveragePt2 <= 0.0 ) {
992 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
1001 G4double& residualExcitationEnergy,
1003 G4int& residualMassNumber,
1004 G4int& residualCharge ) {
1016 if ( ! nucleus )
return false;
1041 sumMasses += 20.0*
MeV;
1045 residualMassNumber--;
1052 #ifdef debugPutOnMassShell
1053 G4cout <<
"ExcitationEnergyPerWoundedNucleon " << ExcitationEPerWoundedNucleon <<
" MeV"<<
G4endl
1054 <<
"\t Residual Charge, MassNumber " << residualCharge <<
" " << residualMassNumber
1055 <<
G4endl <<
"\t Initial Momentum " << nucleusMomentum/
GeV<<
" GeV"
1056 <<
G4endl <<
"\t Residual Momentum " << residualMomentum/
GeV<<
" GeV"<<
G4endl;
1059 residualMomentum.
setPz( 0.0 );
1060 residualMomentum.
setE( 0.0 );
1061 if ( residualMassNumber == 0 ) {
1063 residualExcitationEnergy = 0.0;
1066 GetIonMass( residualCharge, residualMassNumber );
1067 if ( residualMassNumber == 1 ) {
1068 residualExcitationEnergy = 0.0;
1071 sumMasses += std::sqrt(
sqr( residualMass ) + residualMomentum.
perp2() );
1080 const G4int numberOfInvolvedNucleons,
1098 if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 )
return false;
1102 const G4double probDeltaIsobar = 0.10;
1104 G4int maxNumberOfDeltas =
G4int( (sqrtS - sumMasses)/(400.0*
MeV) );
1105 G4int numberOfDeltas = 0;
1107 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1110 if (
G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
1112 if ( ! involvedNucleons[i] )
continue;
1119 G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4;
1128 if ( sqrtS < sumMasses + massDelta - massNuc ) {
1132 sumMasses += ( massDelta - massNuc );
1151 const G4int residualMassNumber,
1152 const G4int numberOfInvolvedNucleons,
1166 if ( ! nucleus )
return false;
1168 if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
1176 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1177 G4Nucleon* aNucleon = involvedNucleons[i];
1178 if ( ! aNucleon )
continue;
1182 const G4int maxNumberOfLoops = 1000;
1183 G4int loopCounter = 0;
1190 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1191 G4Nucleon* aNucleon = involvedNucleons[i];
1192 if ( ! aNucleon )
continue;
1198 if ( x < 0.0 || x > 1.0 ) {
1208 if ( xSum < 0.0 || xSum > 1.0 ) success =
false;
1210 if ( ! success )
continue;
1212 G4double deltaPx = ( ptSum.
x() - pResidual.
x() ) / numberOfInvolvedNucleons;
1213 G4double deltaPy = ( ptSum.
y() - pResidual.
y() ) / numberOfInvolvedNucleons;
1215 if ( residualMassNumber == 0 ) {
1216 delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
1223 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1224 G4Nucleon* aNucleon = involvedNucleons[i];
1225 if ( ! aNucleon )
continue;
1228 if ( residualMassNumber == 0 ) {
1229 if ( x <= 0.0 || x > 1.0 ) {
1234 if ( x <= 0.0 || x > 1.0 || xSum <= 0.0 || xSum > 1.0 ) {
1242 +
sqr( px ) +
sqr( py ) ) / x;
1247 if ( success && residualMassNumber != 0 ) {
1248 mass2 += (
sqr( residualMass ) + pResidual.
perp2() ) / xSum;
1251 #ifdef debugPutOnMassShell
1255 }
while ( ( ! success ) &&
1256 ++loopCounter < maxNumberOfLoops );
1257 if ( loopCounter >= maxNumberOfLoops ) {
1273 const G4bool isProjectileNucleus,
1274 const G4int numberOfInvolvedNucleons,
1289 G4double decayMomentum2 =
sqr( sValue ) +
sqr( projectileMass2 ) +
sqr( targetMass2 )
1290 - 2.0*sValue*projectileMass2 - 2.0*sValue*targetMass2
1291 - 2.0*projectileMass2*targetMass2;
1292 targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
1294 projectileWplus = sqrtS - targetMass2/targetWminus;
1295 G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
1296 G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
1298 if (projectileE - projectilePz > 0.) {
1299 projectileY = 0.5 *
G4Log( (projectileE + projectilePz)/
1300 (projectileE - projectilePz) );
1302 G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
1303 G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
1304 G4double targetY = 0.5 *
G4Log( (targetE + targetPz)/(targetE - targetPz) );
1306 #ifdef debugPutOnMassShell
1307 G4cout <<
"decayMomentum2 " << decayMomentum2 <<
G4endl
1308 <<
"\t targetWminus projectileWplus " << targetWminus <<
" " << projectileWplus <<
G4endl
1309 <<
"\t projectileY targetY " << projectileY <<
" " << targetY <<
G4endl;
1312 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1313 G4Nucleon* aNucleon = involvedNucleons[i];
1314 if ( ! aNucleon )
continue;
1319 G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
1320 G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*
x);
1321 if ( isProjectileNucleus ) {
1322 pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*
x);
1323 e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*
x);
1327 #ifdef debugPutOnMassShell
1328 G4cout <<
"i nY pY nY-AY AY " << i <<
" " << nucleonY <<
" " << projectileY <<
G4endl;
1331 if (
std::abs( nucleonY - nucleusY ) > 2 ||
1332 ( isProjectileNucleus && targetY > nucleonY ) ||
1333 ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
1346 const G4bool isProjectileNucleus,
1349 const G4int residualMassNumber,
1350 const G4int numberOfInvolvedNucleons,
1368 for (
G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1369 G4Nucleon* aNucleon = involvedNucleons[i];
1370 if ( ! aNucleon )
continue;
1372 residual3Momentum -= tmp.
vect();
1376 G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w *
x );
1377 G4double e = w * x / 2.0 + mt2 / ( 2.0 * w *
x );
1379 if ( isProjectileNucleus ) pz *= -1.0;
1386 #ifdef debugPutOnMassShell
1387 G4cout <<
"Target involved nucleon No, name, 4Mom "
1392 G4double residualMt2 =
sqr( residualMass ) +
sqr( residual3Momentum.
x() )
1393 +
sqr( residual3Momentum.
y() );
1395 #ifdef debugPutOnMassShell
1396 G4cout <<
G4endl<<
"w residual3Momentum.z() " << w <<
" " << residual3Momentum.
z() <<
G4endl;
1401 if ( residualMassNumber != 0 ) {
1402 residualPz = -w * residual3Momentum.
z() / 2.0 +
1403 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
1404 residualE = w * residual3Momentum.
z() / 2.0 +
1405 residualMt2 / ( 2.0 * w * residual3Momentum.
z() );
1407 if ( isProjectileNucleus ) residualPz *= -1.0;
1410 residual4Momentum.
setPx( residual3Momentum.
x() );
1411 residual4Momentum.
setPy( residual3Momentum.
y() );
1412 residual4Momentum.
setPz( residualPz );
1413 residual4Momentum.
setE( residualE );
1421 #ifdef debugQGSParticipants
1430 #ifdef debugQGSParticipants
1431 G4cout<<
"Interaction # and its status "
1436 if ( (InterStatus ==
PrD) || (InterStatus ==
TrD) || (InterStatus ==
DD))
1438 #ifdef debugQGSParticipants
1439 G4cout<<
"Simulation of diffractive interaction. "<<InterStatus <<
" PrD/TrD/DD/ND/Qech - 0,1,2,3,4"<<
G4endl;
1444 #ifdef debugQGSParticipants
1445 G4cout<<
"The proj. before inter "
1452 if ( InterStatus ==
PrD )
1455 if ( InterStatus ==
TrD )
1458 if ( InterStatus ==
DD )
1461 #ifdef debugQGSParticipants
1469 if ( InterStatus ==
Qexc )
1471 #ifdef debugQGSParticipants
1472 G4cout<<
"Simulation of interaction with quark exchange."<<
G4endl;
1476 #ifdef debugQGSParticipants
1485 #ifdef debugQGSParticipants
1500 const G4double aHugeValue = 1.0e+10;
1502 #ifdef debugQGSParticipants
1512 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(350.), VqqM_tr(700);
1515 #ifdef debugQGSParticipants
1516 G4cout<<
"Projectile 4 momentum "<<Psum<<G4endl
1517 <<
"Target nucleon momenta at start"<<
G4endl;
1520 std::vector<G4VSplitableHadron*>::iterator i;
1525 Psum += (*i)->Get4Momentum();
1526 #ifdef debugQGSParticipants
1527 G4cout<<
"Nusleus nucleon # and its 4Mom. "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1536 toCms.
rotateZ( -1*Ptmp.phi() );
1537 toCms.
rotateY( -1*Ptmp.theta() );
1542 #ifdef debugQGSParticipants
1544 G4cout<<
"Projectile 4 Mom "<<Projectile4Momentum<<
G4endl;
1552 (*i)->Set4Momentum( tmp );
1553 #ifdef debugQGSParticipants
1554 G4cout<<
"Target nucleon # and 4Mom "<<
" "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1556 Target4Momentum +=
tmp;
1563 #ifdef debugQGSParticipants
1564 G4cout<<
"Sum of target nucleons 4 momentum "<<Target4Momentum<<G4endl<<
G4endl;
1565 G4cout<<
"Target nucleons mom: px, py, z_1, m_i"<<
G4endl;
1585 if ( Mass2 < 0.0 ) {
1588 <<
"Â 4-momentum " << Psum <<
G4endl;
1589 ed <<
"LorentzVector tmp " << tmp <<
" with Mass2 " << Mass2 <<
G4endl;
1590 G4Exception(
"G4QGSParticipants::DeterminePartonMomenta(): 4-momentum with negative mass!",
1593 Mass = std::sqrt( Mass2 );
1597 (*i)->Set4Momentum(tmp);
1598 #ifdef debugQGSParticipants
1599 G4cout<<
"Target nucleons # and mom: "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1612 G4double ProjSumMt(0.), ProjSumMt2perX(0.);
1613 G4double TargSumMt(0.), TargSumMt2perX(0.);
1629 const G4int maxNumberOfAttempts = 1000;
1632 attempt++;
if( attempt == 100*(attempt/100) ) {SigPt/=2.;}
1634 ProjSumMt=0.; ProjSumMt2perX=0.;
1635 TargSumMt=0.; TargSumMt2perX=0.;
1639 #ifdef debugQGSParticipants
1640 G4cout<<
"attempt ------------------------ "<<attempt<<
G4endl;
1647 G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1650 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1653 #ifdef debugQGSParticipants
1658 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1659 Mt = std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1663 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1666 NumberOfUnsampledSeaQuarks--;
1667 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1672 #ifdef debugQGSParticipants
1678 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1679 Mt = std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1683 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1686 NumberOfUnsampledSeaQuarks--;
1687 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1690 #ifdef debugQGSParticipants
1697 #ifdef debugQGSParticipants
1702 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1703 Mt = std::sqrt(aPtVector.
mag2()+
sqr(VqM_pr));
1707 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1710 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1716 #ifdef debugQGSParticipants
1722 Mt = std::sqrt(
sqr(SumPx) +
sqr(SumPy) +
sqr(VaqM_pr) );
1726 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1729 #ifdef debugQGSParticipants
1730 G4cout<<
" "<<tmp<<
" "<<SumZ+(1.-SumZ)<<
" (z-fraction)"<<G4endl;
1740 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1741 #ifdef debugQGSParticipants
1742 G4cout<<
"nSeaPair of target N "<<nSeaPair<<G4endl
1743 <<
"Target nucleon 4Mom "<<(*i)->Get4Momentum()<<
G4endl;
1746 SumPx = (*i)->Get4Momentum().px() * (-1.);
1747 SumPy = (*i)->Get4Momentum().py() * (-1.);
1751 NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1754 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1756 aParton = (*i)->GetNextParton();
1757 #ifdef debugQGSParticipants
1762 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1763 Mt=std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1767 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1769 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1771 NumberOfUnsampledSeaQuarks--;
1772 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1776 aParton = (*i)->GetNextAntiParton();
1777 #ifdef debugQGSParticipants
1783 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1784 Mt=std::sqrt(aPtVector.
mag2()+
sqr(Qmass));
1788 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1790 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1792 NumberOfUnsampledSeaQuarks--;
1793 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1796 #ifdef debugQGSParticipants
1802 aParton = (*i)->GetNextParton();
1803 #ifdef debugQGSParticipants
1808 SumPx += aPtVector.
x(); SumPy += aPtVector.
y();
1809 Mt=std::sqrt(aPtVector.
mag2()+
sqr(VqM_tr));
1813 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1815 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1817 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1822 aParton = (*i)->GetNextAntiParton();
1823 #ifdef debugQGSParticipants
1825 G4cout<<
" "<<tmp<<
" "<<SumZw<<
" (sum z-fracs) "<<SumZ<<
" (total z-sum) "<<
G4endl;
1829 Mt=std::sqrt(
sqr(SumPx) +
sqr(SumPy) +
sqr(VqqM_tr) );
1832 tmp.
setPz((*i)->Get4Momentum().pz()*(1.0 - SumZ));
1834 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1837 #ifdef debugQGSParticipants
1838 G4cout<<
" "<<tmp<<
" "<<SumZw<<
" "<<1.0<<
" "<<(*i)->Get4Momentum().
pz()<<
G4endl;
1843 if( ProjSumMt + TargSumMt > SqrtS ) {
1844 Success =
false;
continue;}
1845 if( std::sqrt(ProjSumMt2perX) + std::sqrt(TargSumMt2perX) > SqrtS ) {
1846 Success =
false;
continue;}
1848 }
while( (!Success) &&
1849 attempt < maxNumberOfAttempts );
1851 if ( attempt >= maxNumberOfAttempts ) {
1858 - 2.0*S*ProjSumMt2perX - 2.0*S*TargSumMt2perX - 2.0*ProjSumMt2perX*TargSumMt2perX;
1860 G4double targetWminus=( S - ProjSumMt2perX + TargSumMt2perX + std::sqrt( DecayMomentum2 ))/2.0/SqrtS;
1861 G4double projectileWplus = SqrtS - TargSumMt2perX/targetWminus;
1867 #ifdef debugQGSParticipants
1868 G4cout<<
"Backward transformation ===================="<<
G4endl;
1872 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1875 #ifdef debugQGSParticipants
1880 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1881 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
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));
1897 #ifdef debugQGSParticipants
1904 #ifdef debugQGSParticipants
1908 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1909 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1916 #ifdef debugQGSParticipants
1921 Tmp.
setPz(projectileWplus*z/2.0 - Tmp.
e()/(2.0*z*projectileWplus));
1922 Tmp.
setE( projectileWplus*z/2.0 + Tmp.
e()/(2.0*z*projectileWplus));
1927 #ifdef debugQGSParticipants
1937 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1938 #ifdef debugQGSParticipants
1939 G4cout<<
"nSeaPair of target and N# "<<nSeaPair<<
" "<<NuclNo<<
G4endl;
1942 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1944 aParton = (*i)->GetNextParton();
1945 #ifdef debugQGSParticipants
1949 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1950 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1955 aParton = (*i)->GetNextAntiParton();
1956 #ifdef debugQGSParticipants
1961 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1962 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1966 #ifdef debugQGSParticipants
1973 aParton = (*i)->GetNextParton();
1974 #ifdef debugQGSParticipants
1978 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1979 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1985 aParton = (*i)->GetNextAntiParton();
1986 #ifdef debugQGSParticipants
1991 Tmp.
setPz(-targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1992 Tmp.
setE( targetWminus*z/2.0 + Tmp.
e()/(2.0*z*targetWminus));
1996 #ifdef debugQGSParticipants
2014 const G4int maxNumberOfLoops = 1000;
2015 G4int loopCounter = 0;
2021 }
while( ( r12 > 1.) &&
2022 ++loopCounter < maxNumberOfLoops );
2023 if ( loopCounter >= maxNumberOfLoops ) {
2035 #ifdef debugQGSParticipants
2040 #ifdef debugQGSParticipants
2041 G4cout<<
"BAD situation: theProjectileSplitable is NULL ! Returning immediately"<<
G4endl;
2046 #ifdef debugQGSParticipants
2058 #ifdef debugQGSParticipants
2075 for (
G4int i = 0; i < N_EnvTarg; i++ ) {
2080 #ifdef debugQGSParticipants
2090 #ifdef debugQGSParticipants
2108 #ifdef debugQGSParticipants
2109 G4cout<<
"Strings created in soft interactions"<<
G4endl;
2111 std::vector<G4InteractionContent*>::iterator i;
2119 #ifdef debugQGSParticipants
2120 G4cout<<
"An interaction # and soft coll. # "<<IntNo<<
" "
2133 #ifdef debugQGSParticipants
2148 #ifdef debugQGSParticipants
2156 #ifdef debugQGSParticipants
2172 #ifdef debugQGSParticipants
2173 G4cout <<
"Sum of strings 4 momenta " << str4Mom << G4endl<<
G4endl;
2182 #ifdef debugQGSParticipants
2186 #ifdef debugQGSParticipants
2197 #ifdef debugQGSParticipants
2221 residualMomentum +=
tmp;
2243 const G4int maxNumberOfLoops = 1000;
2244 G4int loopCounter = 0;
2261 if(SumMasses > Mass) {Chigh=
C;}
2264 }
while( (Chigh-Clow > 0.01) &&
2265 ++loopCounter < maxNumberOfLoops );
2266 if ( loopCounter >= maxNumberOfLoops ) {
2267 #ifdef debugQGSParticipants
2288 #ifdef debugQGSParticipants
2289 G4cout <<
"End GetResiduals -----------------" <<
G4endl;
2298 std::vector<G4InteractionContent*>::iterator i;
2313 #ifdef debugQGSParticipants
2321 #ifdef debugQGSParticipants
2328 #ifdef debugQGSParticipants
2336 #ifdef debugQGSParticipants
2348 #ifdef debugQGSParticipants
2365 if ((!(aPrimaryMomentum.e()>-1)) && (!(aPrimaryMomentum.e()<1)) )
2368 "G4GammaParticipants::SelectInteractions: primary nan energy.");
2370 G4double S = (aPrimaryMomentum + aTargetNMomentum).mag2();
2374 #ifdef debugQGSParticipants
2375 G4cout <<
"Gamma projectile - SelectInteractions " <<
G4endl;
2377 G4cout <<
"SqrtS ThresholdMass ModelMode " <<std::sqrt(S)<<
" "<<ThresholdMass<<
" "<<
ModelMode<<
G4endl;
2392 G4int totalCuts = 0;
2401 {
if(NucleonNo == theCurrent)
break; NucleonNo++; }
2406 pNucleon->
Hit(aTarget);