48 G4bool& incidentHasChanged,
71 G4bool veryForward =
false;
79 G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
80 targetMass*targetMass +
81 2.0*targetMass*etOriginal);
85 if (currentMass == 0.0 && targetMass == 0.0) {
88 currentParticle = *vec[0];
89 targetParticle = *vec[1];
90 for (i = 0; i < (vecLen-2); ++i) *vec[i] = *vec[i+2];
92 for (
G4int j = 0; j < vecLen; j++)
delete vec[j];
95 "G4RPGTwoCluster::ReactionStage : Negative number of particles");
102 incidentHasChanged =
true;
103 targetHasChanged =
true;
117 G4int forwardCount = 1;
123 G4int backwardCount = 1;
129 for (i = 0; i < vecLen; ++i) {
130 if (vec[i]->GetSide() < 0) vec[i]->SetSide(-1);
133 if (vec[i]->GetSide() == -1) {
135 backwardMass += vec[i]->GetMass()/
GeV;
138 forwardMass += vec[i]->GetMass()/
GeV;
143 G4double term1 =
G4Log(centerofmassEnergy*centerofmassEnergy);
144 if(term1 < 0) term1 = 0.0001;
149 xtarg = afc * (a13-1.0) * (2*backwardCount+vecLen+2)/2.0;
151 xtarg = afc * (a13-1.0) * (2*backwardCount);
153 if( xtarg <= 0.0 )xtarg = 0.01;
156 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
160 if( nuclearExcitationCount > 0 )
163 const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
168 for( i=0; i<nuclearExcitationCount; ++i )
186 else if( ran < 0.6819 )
206 G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
207 G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
208 G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
214 ed <<
" While count exceeded " <<
G4endl;
216 while( eAvailable <= 0.0 ) {
223 secondaryDeleted =
false;
224 for( i=(vecLen-1); i>=0; --i )
226 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
228 pMass = vec[i]->GetMass()/
GeV;
229 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
231 forwardEnergy += pMass;
232 forwardMass -= pMass;
233 secondaryDeleted =
true;
236 else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
238 pMass = vec[i]->GetMass()/
GeV;
239 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
241 backwardEnergy += pMass;
242 backwardMass -= pMass;
243 secondaryDeleted =
true;
248 if( secondaryDeleted )
250 delete vec[vecLen-1];
256 if( vecLen == 0 )
return false;
257 if( targetParticle.
GetSide() == -1 )
260 targetParticle = *vec[0];
261 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
263 backwardEnergy += pMass;
264 backwardMass -= pMass;
265 secondaryDeleted =
true;
267 else if( targetParticle.
GetSide() == 1 )
270 targetParticle = *vec[0];
271 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
273 forwardEnergy += pMass;
274 forwardMass -= pMass;
275 secondaryDeleted =
true;
278 if( secondaryDeleted )
280 delete vec[vecLen-1];
285 if( currentParticle.
GetSide() == -1 )
288 currentParticle = *vec[0];
289 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
291 backwardEnergy += pMass;
292 backwardMass -= pMass;
293 secondaryDeleted =
true;
295 else if( currentParticle.
GetSide() == 1 )
298 currentParticle = *vec[0];
299 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
301 forwardEnergy += pMass;
302 forwardMass -= pMass;
303 secondaryDeleted =
true;
305 if( secondaryDeleted )
307 delete vec[vecLen-1];
315 eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
327 const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
328 const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
331 if (forwardCount < 1 || backwardCount < 1)
return false;
334 if (forwardCount > 1) {
339 if( backwardCount > 1 ) {
346 eda <<
" While count exceeded " <<
G4endl;
347 while( rmc+rmd > centerofmassEnergy ) {
354 if( (rmc <= forwardMass) && (rmd <= backwardMass) )
356 G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
362 rmc = 0.1*forwardMass + 0.9*rmc;
363 rmd = 0.1*backwardMass + 0.9*rmd;
368 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
372 pseudoParticle[1].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
380 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
381 pseudoParticle[1].
Lorentz( pseudoParticle[1], pseudoParticle[0] );
382 pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
389 G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
391 pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
392 pf = std::sqrt(
std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
396 pseudoParticle[3].SetMass( rmc*GeV );
397 pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
399 pseudoParticle[4].SetMass( rmd*GeV );
400 pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
408 G4double pin = pseudoParticle[1].GetMomentum().mag()/
GeV;
414 G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
419 pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
420 pf*sintheta*std::sin(phi)*GeV,
422 pseudoParticle[4].SetMomentum( -pseudoParticle[3].GetMomentum());
428 if( nuclearExcitationCount > 0 )
433 if( ekOriginal <= 5.0 )
435 ekit1 *= ekOriginal*ekOriginal/25.0;
436 ekit2 *= ekOriginal*ekOriginal/25.0;
439 for( i=0; i<vecLen; ++i )
441 if( vec[i]->GetSide() == -2 )
444 vec[i]->SetKineticEnergy( kineticE*GeV );
446 G4double totalE = kineticE*GeV + vMass;
447 pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
449 G4double sint = std::sqrt(1.0-cost*cost);
451 vec[i]->SetMomentum( pp*sint*std::cos(phi)*MeV,
452 pp*sint*std::sin(phi)*MeV,
454 vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
463 currentParticle.
SetMomentum( pseudoParticle[3].GetMomentum() );
464 currentParticle.
SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
466 targetParticle.
SetMomentum( pseudoParticle[4].GetMomentum() );
467 targetParticle.
SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
469 pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
470 pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
471 pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
473 pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
474 pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
475 pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
479 if( forwardCount > 1 )
483 G4bool constantCrossSection =
true;
485 if( currentParticle.
GetSide() == 1 )
486 tempV.
SetElement( tempLen++, ¤tParticle );
487 if( targetParticle.
GetSide() == 1 )
488 tempV.
SetElement( tempLen++, &targetParticle );
489 for( i=0; i<vecLen; ++i )
491 if( vec[i]->GetSide() == 1 )
497 vec[i]->SetSide( -1 );
505 constantCrossSection, tempV, tempLen );
506 if( currentParticle.
GetSide() == 1 )
507 currentParticle.
Lorentz( currentParticle, pseudoParticle[5] );
508 if( targetParticle.
GetSide() == 1 )
509 targetParticle.
Lorentz( targetParticle, pseudoParticle[5] );
510 for( i=0; i<vecLen; ++i )
512 if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
517 if( backwardCount > 1 )
521 G4bool constantCrossSection =
true;
523 if( currentParticle.
GetSide() == -1 )
524 tempV.
SetElement( tempLen++, ¤tParticle );
525 if( targetParticle.
GetSide() == -1 )
526 tempV.
SetElement( tempLen++, &targetParticle );
527 for( i=0; i<vecLen; ++i )
529 if( vec[i]->GetSide() == -1 )
535 vec[i]->SetSide( -2 );
536 vec[i]->SetKineticEnergy( 0.0 );
537 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
545 constantCrossSection, tempV, tempLen );
546 if( currentParticle.
GetSide() == -1 )
547 currentParticle.
Lorentz( currentParticle, pseudoParticle[6] );
548 if( targetParticle.
GetSide() == -1 )
549 targetParticle.
Lorentz( targetParticle, pseudoParticle[6] );
550 for( i=0; i<vecLen; ++i )
552 if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
561 currentParticle.
Lorentz( currentParticle, pseudoParticle[2] );
562 targetParticle.
Lorentz( targetParticle, pseudoParticle[2] );
563 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
584 for( i=0; i<vecLen; ++i )
586 if( vec[i]->GetDefinition() == leadingStrangeParticle.
GetDefinition() )
597 if( ((leadMass < protonMass) && (targetParticle.
GetMass()/MeV < protonMass)) ||
598 ((leadMass >= protonMass) && (targetParticle.
GetMass()/MeV >= protonMass)) )
611 targetHasChanged =
true;
626 incidentHasChanged =
true;
634 std::pair<G4int, G4int> finalStateNucleons =
637 G4int protonsInFinalState = finalStateNucleons.first;
638 G4int neutronsInFinalState = finalStateNucleons.second;
640 G4int numberofFinalStateNucleons =
641 protonsInFinalState + neutronsInFinalState;
647 numberofFinalStateNucleons++;
649 numberofFinalStateNucleons =
std::max(1, numberofFinalStateNucleons);
659 pseudoParticle[4].SetMass( mOriginal*GeV );
660 pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
661 pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
666 if(numberofFinalStateNucleons == 1) diff = 0;
667 pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
668 pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
669 pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
672 pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/
GeV;
674 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
675 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
676 pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
681 tempR[0] = currentParticle;
682 tempR[1] = targetParticle;
683 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
687 G4bool constantCrossSection =
true;
689 for( i=0; i<vecLen+2; ++i )tempV.
SetElement( tempLen++, &tempR[i] );
695 pseudoParticle[5].GetTotalEnergy()/MeV,
696 constantCrossSection, tempV, tempLen );
699 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
701 constantCrossSection, tempV, tempLen );
703 theoreticalKinetic = 0.0;
704 for( i=0; i<vecLen+2; ++i )
706 pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
707 pseudoParticle[7].SetMass( tempV[i]->GetMass() );
708 pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
709 pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
710 theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/
GeV;
717 theoreticalKinetic -=
719 for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/
GeV;
723 for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/
GeV;
728 if( simulatedKinetic != 0.0 )
730 wgt = (theoreticalKinetic)/simulatedKinetic;
734 if( pp1 < 0.001*MeV ) {
744 if( pp1 < 0.001*MeV ) {
751 for( i=0; i<vecLen; ++i )
753 vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
754 pp = vec[i]->GetTotalMomentum()/
MeV;
755 pp1 = vec[i]->GetMomentum().mag()/
MeV;
758 vec[i]->SetMomentum( iso.
x(), iso.
y(), iso.
z() );
760 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
766 Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
767 modifiedOriginal, originalIncident, targetNucleus,
768 currentParticle, targetParticle, vec, vecLen );
774 if( atomicWeight >= 1.5 )
805 if( epnb >= pnCutOff )
807 npnb =
G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
808 if( numberofFinalStateNucleons + npnb > atomicWeight )
809 npnb =
G4int(atomicWeight - numberofFinalStateNucleons);
810 npnb =
std::min( npnb, 127-vecLen );
812 if( edta >= dtaCutOff )
814 ndta =
G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
815 ndta =
std::min( ndta, 127-vecLen );
817 if (npnb == 0 && ndta == 0) npnb = 1;
822 PinNucleus, NinNucleus, targetNucleus,
832 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
835 currentParticle.
SetTOF( 1.0 );
static G4PionMinus * PionMinus()
void Initialize(G4int items)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static G4Lambda * Lambda()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
std::ostringstream G4ExceptionDescription
void AddBlackTrackParticles(const G4double, const G4int, const G4double, const G4int, const G4ReactionProduct &, G4int, G4int, const G4Nucleus &, G4FastVector< G4ReactionProduct, 256 > &, G4int &)
G4double A13(G4double A) const
static constexpr double MeV
G4double GetDTABlackTrackEnergy() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetMomentum(const G4double x, const G4double y, const G4double z)
static G4PionPlus * PionPlus()
const G4ParticleDefinition * GetDefinition() const
static G4Proton * Proton()
G4double GetPDGMass() const
G4double G4Log(G4double x)
void SetTotalEnergy(const G4double en)
void SetMass(const G4double mas)
void SetTOF(const G4double t)
G4double GetTotalEnergy() const
void SetKineticEnergy(const G4double en)
G4double GetKineticEnergy() const
static G4Pow * GetInstance()
G4int GetBaryonNumber() const
static constexpr double twopi
G4ThreeVector Isotropic(const G4double &)
G4double GetTotalMomentum() const
G4ThreeVector GetMomentum() const
void SetSide(const G4int sid)
void Rotate(const G4double numberofFinalStateNucleons, const G4ThreeVector &temp, const G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, const G4Nucleus &targetNucleus, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
G4double GetAnnihilationPNBlackTrackEnergy() const
void SetElement(G4int anIndex, Type *anElement)
G4double GetAnnihilationDTABlackTrackEnergy() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
const G4ParticleDefinition * GetDefinition() const
static G4Neutron * Neutron()
std::pair< G4int, G4int > GetFinalStateNucleons(const G4DynamicParticle *originalTarget, const G4FastVector< G4ReactionProduct, 256 > &vec, const G4int &vecLen)
static G4PionZero * PionZero()
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
G4long G4Poisson(G4double mean)
void SetNewlyAdded(const G4bool f)
G4double GetPNBlackTrackEnergy() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
static constexpr double GeV
T min(const T t1, const T t2)
brief Return the smallest of the two arguments