46 G4bool& incidentHasChanged,
69 G4bool veryForward =
false;
77 G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
78 targetMass*targetMass +
79 2.0*targetMass*etOriginal);
83 if (currentMass == 0.0 && targetMass == 0.0) {
86 currentParticle = *vec[0];
87 targetParticle = *vec[1];
88 for (i = 0; i < (vecLen-2); ++i) *vec[i] = *vec[i+2];
90 for (
G4int j = 0; j < vecLen; j++)
delete vec[j];
93 "G4RPGTwoCluster::ReactionStage : Negative number of particles");
100 incidentHasChanged =
true;
101 targetHasChanged =
true;
115 G4int forwardCount = 1;
121 G4int backwardCount = 1;
127 for (i = 0; i < vecLen; ++i) {
128 if (vec[i]->GetSide() < 0) vec[i]->SetSide(-1);
131 if (vec[i]->GetSide() == -1) {
133 backwardMass += vec[i]->GetMass()/
GeV;
136 forwardMass += vec[i]->GetMass()/
GeV;
141 G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
142 if(term1 < 0) term1 = 0.0001;
143 const G4double afc = 0.312 + 0.2 * std::log(term1);
146 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
148 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
149 if( xtarg <= 0.0 )xtarg = 0.01;
152 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
156 if( nuclearExcitationCount > 0 )
159 const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
164 for( i=0; i<nuclearExcitationCount; ++i )
182 else if( ran < 0.6819 )
202 G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
203 G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
204 G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
208 while( eAvailable <= 0.0 )
210 secondaryDeleted =
false;
211 for( i=(vecLen-1); i>=0; --i )
213 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
215 pMass = vec[i]->GetMass()/
GeV;
216 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
218 forwardEnergy += pMass;
219 forwardMass -= pMass;
220 secondaryDeleted =
true;
223 else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
225 pMass = vec[i]->GetMass()/
GeV;
226 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
228 backwardEnergy += pMass;
229 backwardMass -= pMass;
230 secondaryDeleted =
true;
235 if( secondaryDeleted )
237 delete vec[vecLen-1];
243 if( vecLen == 0 )
return false;
244 if( targetParticle.
GetSide() == -1 )
247 targetParticle = *vec[0];
248 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
250 backwardEnergy += pMass;
251 backwardMass -= pMass;
252 secondaryDeleted =
true;
254 else if( targetParticle.
GetSide() == 1 )
257 targetParticle = *vec[0];
258 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
260 forwardEnergy += pMass;
261 forwardMass -= pMass;
262 secondaryDeleted =
true;
265 if( secondaryDeleted )
267 delete vec[vecLen-1];
272 if( currentParticle.
GetSide() == -1 )
275 currentParticle = *vec[0];
276 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
278 backwardEnergy += pMass;
279 backwardMass -= pMass;
280 secondaryDeleted =
true;
282 else if( currentParticle.
GetSide() == 1 )
285 currentParticle = *vec[0];
286 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
288 forwardEnergy += pMass;
289 forwardMass -= pMass;
290 secondaryDeleted =
true;
292 if( secondaryDeleted )
294 delete vec[vecLen-1];
302 eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
314 const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
315 const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
318 if (forwardCount < 1 || backwardCount < 1)
return false;
321 if (forwardCount > 1) {
323 rmc += std::pow(-std::log(1.0-
G4UniformRand()),1./cpar[ntc])/gpar[ntc];
326 if( backwardCount > 1 ) {
328 rmd += std::pow(-std::log(1.0-
G4UniformRand()),1./cpar[ntc])/gpar[ntc];
331 while( rmc+rmd > centerofmassEnergy )
333 if( (rmc <= forwardMass) && (rmd <= backwardMass) )
335 G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
341 rmc = 0.1*forwardMass + 0.9*rmc;
342 rmd = 0.1*backwardMass + 0.9*rmd;
347 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
351 pseudoParticle[1].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
359 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
360 pseudoParticle[1].
Lorentz( pseudoParticle[1], pseudoParticle[0] );
361 pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
368 G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
370 pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
371 pf = std::sqrt(
std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
375 pseudoParticle[3].SetMass( rmc*GeV );
376 pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
378 pseudoParticle[4].SetMass( rmd*GeV );
379 pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
387 G4double pin = pseudoParticle[1].GetMomentum().mag()/
GeV;
389 G4double factor = 1.0 - std::exp(-dtb);
393 G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
398 pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
399 pf*sintheta*std::sin(phi)*GeV,
401 pseudoParticle[4].SetMomentum( -pseudoParticle[3].GetMomentum());
407 if( nuclearExcitationCount > 0 )
412 if( ekOriginal <= 5.0 )
414 ekit1 *= ekOriginal*ekOriginal/25.0;
415 ekit2 *= ekOriginal*ekOriginal/25.0;
418 for( i=0; i<vecLen; ++i )
420 if( vec[i]->GetSide() == -2 )
423 vec[i]->SetKineticEnergy( kineticE*GeV );
425 G4double totalE = kineticE*GeV + vMass;
426 pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
428 G4double sint = std::sqrt(1.0-cost*cost);
430 vec[i]->SetMomentum( pp*sint*std::cos(phi)*MeV,
431 pp*sint*std::sin(phi)*MeV,
433 vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
442 currentParticle.
SetMomentum( pseudoParticle[3].GetMomentum() );
443 currentParticle.
SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
445 targetParticle.
SetMomentum( pseudoParticle[4].GetMomentum() );
446 targetParticle.
SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
448 pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
449 pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
450 pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
452 pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
453 pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
454 pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
458 if( forwardCount > 1 )
462 G4bool constantCrossSection =
true;
464 if( currentParticle.
GetSide() == 1 )
465 tempV.
SetElement( tempLen++, ¤tParticle );
466 if( targetParticle.
GetSide() == 1 )
467 tempV.
SetElement( tempLen++, &targetParticle );
468 for( i=0; i<vecLen; ++i )
470 if( vec[i]->GetSide() == 1 )
476 vec[i]->SetSide( -1 );
484 constantCrossSection, tempV, tempLen );
485 if( currentParticle.
GetSide() == 1 )
486 currentParticle.
Lorentz( currentParticle, pseudoParticle[5] );
487 if( targetParticle.
GetSide() == 1 )
488 targetParticle.
Lorentz( targetParticle, pseudoParticle[5] );
489 for( i=0; i<vecLen; ++i )
491 if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
496 if( backwardCount > 1 )
500 G4bool constantCrossSection =
true;
502 if( currentParticle.
GetSide() == -1 )
503 tempV.
SetElement( tempLen++, ¤tParticle );
504 if( targetParticle.
GetSide() == -1 )
505 tempV.
SetElement( tempLen++, &targetParticle );
506 for( i=0; i<vecLen; ++i )
508 if( vec[i]->GetSide() == -1 )
514 vec[i]->SetSide( -2 );
515 vec[i]->SetKineticEnergy( 0.0 );
516 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
524 constantCrossSection, tempV, tempLen );
525 if( currentParticle.
GetSide() == -1 )
526 currentParticle.
Lorentz( currentParticle, pseudoParticle[6] );
527 if( targetParticle.
GetSide() == -1 )
528 targetParticle.
Lorentz( targetParticle, pseudoParticle[6] );
529 for( i=0; i<vecLen; ++i )
531 if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
540 currentParticle.
Lorentz( currentParticle, pseudoParticle[2] );
541 targetParticle.
Lorentz( targetParticle, pseudoParticle[2] );
542 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
563 for( i=0; i<vecLen; ++i )
565 if( vec[i]->GetDefinition() == leadingStrangeParticle.
GetDefinition() )
576 if( ((leadMass < protonMass) && (targetParticle.
GetMass()/MeV < protonMass)) ||
577 ((leadMass >= protonMass) && (targetParticle.
GetMass()/MeV >= protonMass)) )
590 targetHasChanged =
true;
605 incidentHasChanged =
true;
613 std::pair<G4int, G4int> finalStateNucleons =
616 G4int protonsInFinalState = finalStateNucleons.first;
617 G4int neutronsInFinalState = finalStateNucleons.second;
619 G4int numberofFinalStateNucleons =
620 protonsInFinalState + neutronsInFinalState;
626 numberofFinalStateNucleons++;
628 numberofFinalStateNucleons =
std::max(1, numberofFinalStateNucleons);
638 pseudoParticle[4].SetMass( mOriginal*GeV );
639 pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
640 pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
645 if(numberofFinalStateNucleons == 1) diff = 0;
646 pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
647 pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
648 pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
651 pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/
GeV;
653 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
654 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
655 pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
660 tempR[0] = currentParticle;
661 tempR[1] = targetParticle;
662 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
666 G4bool constantCrossSection =
true;
668 for( i=0; i<vecLen+2; ++i )tempV.
SetElement( tempLen++, &tempR[i] );
674 pseudoParticle[5].GetTotalEnergy()/MeV,
675 constantCrossSection, tempV, tempLen );
678 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
680 constantCrossSection, tempV, tempLen );
682 theoreticalKinetic = 0.0;
683 for( i=0; i<vecLen+2; ++i )
685 pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
686 pseudoParticle[7].SetMass( tempV[i]->GetMass() );
687 pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
688 pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
689 theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/
GeV;
696 theoreticalKinetic -=
698 for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/
GeV;
702 for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/
GeV;
707 if( simulatedKinetic != 0.0 )
709 wgt = (theoreticalKinetic)/simulatedKinetic;
713 if( pp1 < 0.001*MeV ) {
723 if( pp1 < 0.001*MeV ) {
730 for( i=0; i<vecLen; ++i )
732 vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
733 pp = vec[i]->GetTotalMomentum()/
MeV;
734 pp1 = vec[i]->GetMomentum().mag()/
MeV;
737 vec[i]->SetMomentum( iso.
x(), iso.
y(), iso.
z() );
739 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
745 Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
746 modifiedOriginal, originalIncident, targetNucleus,
747 currentParticle, targetParticle, vec, vecLen );
753 if( atomicWeight >= 1.5 )
784 if( epnb >= pnCutOff )
786 npnb =
G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
787 if( numberofFinalStateNucleons + npnb > atomicWeight )
788 npnb =
G4int(atomicWeight - numberofFinalStateNucleons);
789 npnb =
std::min( npnb, 127-vecLen );
791 if( edta >= dtaCutOff )
793 ndta =
G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
794 ndta =
std::min( ndta, 127-vecLen );
796 if (npnb == 0 && ndta == 0) npnb = 1;
801 PinNucleus, NinNucleus, targetNucleus,
811 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
814 currentParticle.
SetTOF( 1.0 );
void SetElement(G4int anIndex, Type *anElement)
G4long G4Poisson(G4double mean)
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 GetTotalMomentum() const
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4double GetAnnihilationPNBlackTrackEnergy() const
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetSide(const G4int sid)
G4double GetDTABlackTrackEnergy() const
void Initialize(G4int items)
G4ParticleDefinition * GetDefinition() const
void SetNewlyAdded(const G4bool f)
std::pair< G4int, G4int > GetFinalStateNucleons(const G4DynamicParticle *originalTarget, const G4FastVector< G4ReactionProduct, 256 > &vec, const G4int &vecLen)
void SetMass(const G4double mas)
const G4ParticleDefinition * GetDefinition() const
G4double GetAnnihilationDTABlackTrackEnergy() const
void SetTotalEnergy(const G4double en)
static G4Proton * Proton()
static G4PionPlus * PionPlus()
static G4Neutron * Neutron()
static G4PionZero * PionZero()
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4PionMinus * PionMinus()
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
G4ThreeVector GetMomentum() const
void AddBlackTrackParticles(const G4double, const G4int, const G4double, const G4int, const G4ReactionProduct &, G4int, G4int, const G4Nucleus &, G4FastVector< G4ReactionProduct, 256 > &, G4int &)
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
void SetTOF(const G4double t)
static G4Lambda * Lambda()
G4ThreeVector Isotropic(const G4double &)
G4double GetPNBlackTrackEnergy() const
G4int GetBaryonNumber() const