83 fP0_d( 0.0 ), fP0_dbar( 0.0 ), secID( -1 ) {
103 G4double ekin = std::sqrt( pz*pz + mproj*mproj ) - mproj;
105 fP0_dbar = 130.0 / ( 1.0 + std::exp( 21.6 - std::log( 0.001*ekin )/0.089 ) );
106 fP0_d = 118.1 * ( 1.0 + std::exp( 5.53 - std::log( 0.001*ekin )/0.43 ) );
123 std::vector< std::pair< G4int, G4ThreeVector > >
proton;
124 std::vector< std::pair< G4int, G4ThreeVector > >
neutron;
125 std::vector< std::pair< G4int, G4ThreeVector > > antiproton;
126 std::vector< std::pair< G4int, G4ThreeVector > > antineutron;
127 for (
unsigned int i = 0; i < result->size(); ++i ) {
128 G4int pdgid = result->operator[](i)->GetDefinition()->GetPDGEncoding();
129 if ( pdgid == 2212 ) {
130 proton.push_back( std::make_pair( i, result->operator[](i)->GetMomentum() ) );
131 result->erase( result->begin() + i );
134 for (
unsigned int i = 0; i < result->size(); ++i ) {
135 G4int pdgid = result->operator[](i)->GetDefinition()->GetPDGEncoding();
136 if ( pdgid == 2112 ) {
137 neutron.push_back( std::make_pair( i, result->operator[](i)->GetMomentum() ) );
138 result->erase( result->begin() + i );
141 for (
unsigned int i = 0; i < result->size(); ++i ) {
142 G4int pdgid = result->operator[](i)->GetDefinition()->GetPDGEncoding();
143 if ( pdgid == -2212 ) {
144 antiproton.push_back( std::make_pair( i, result->operator[](i)->GetMomentum() ) );
145 result->erase( result->begin() + i );
148 for (
unsigned int i = 0; i < result->size(); ++i ) {
149 G4int pdgid = result->operator[](i)->GetDefinition()->GetPDGEncoding();
150 if ( pdgid == -2112 ) {
151 antineutron.push_back( std::make_pair( i, result->operator[](i)->GetMomentum() ) );
152 result->erase( result->begin() + i );
156 for (
unsigned int i = 0; i <
proton.size(); ++i ) {
157 if (
proton.at(i).first == -1 )
continue;
161 if ( partner1 == -1 ) {
166 G4double totalEnergy = std::sqrt( p1.
mag()*p1.
mag() + massp*massp );
170 result->push_back( finalp );
175 neutron.at(partner1).first = -1;
178 for (
unsigned int i = 0; i <
neutron.size(); ++i ) {
179 if (
neutron.at(i).first == -1 )
continue;
185 G4double totalEnergy = std::sqrt( p2.
mag()*p2.
mag() + massn*massn );
189 result->push_back( finaln );
192 for (
unsigned int i = 0; i < antiproton.size(); ++i ) {
193 if ( antiproton.at(i).first == -1 )
continue;
197 if ( partner1 == -1 ) {
202 G4double totalEnergy = std::sqrt( p1.
mag()*p1.
mag() + massp*massp );
206 result->push_back( finalpbar );
211 antineutron.at(partner1).first = -1;
214 for (
unsigned int i = 0; i < antineutron.size(); ++i ) {
215 if ( antineutron.at(i).first == -1 )
continue;
221 G4double totalEnergy = std::sqrt( p2.
mag()*p2.
mag() + massn*massn );
225 result->push_back( finalnbar );
240 G4double totalEnergy = std::sqrt( psum.
mag()*psum.
mag() + massd*massd );
245 result->push_back( finaldeut );
252 G4double totalEnergy = std::sqrt( psum.
mag()*psum.
mag() + massd*massd );
255 finalantideut->
SetMass( massd );
257 result->push_back( finalantideut );
263 std::vector< std::pair< G4int, G4ThreeVector > > &
neutron,
269 for (
unsigned int j = 0; j <
neutron.size(); ++j ) {
270 if (
neutron.at(j).first == -1 )
continue;
293 if ( charge > 0 )
return ( deltaP <
fP0_d );
301 return GetPcm( p1.
x(), p1.
y(), p1.
z(), m1, p2.
x(), p2.
y(), p2.
z(),
m2 );
309 return std::sqrt( (scm - (m1-
m2)*(m1-
m2))*(scm - (m1+
m2)*(m1+
m2)) ) / (2.0*std::sqrt( scm ));
316 G4double E1 = std::sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1 );
317 G4double E2 = std::sqrt( p2x*p2x + p2y*p2y + p2z*p2z +
m2*
m2 );
318 return (E1+E2)*(E1+E2) - (p1x+p2x)*(p1x+p2x) - (p1y+p2y)*(p1y+p2y) - (p1z+p2z)*(p1z+p2z);
std::vector< G4ReactionProduct * > G4ReactionProductVector
static constexpr double m2
void GenerateDeuterons(G4ReactionProductVector *result)
~G4CRCoalescence() override
G4double GetS(G4double p1x, G4double p1y, G4double p1z, G4double m1, G4double p2x, G4double p2y, G4double p2z, G4double m2)
G4bool Coalescence(const G4ThreeVector &p1, G4double m1, const G4ThreeVector &p2, G4double m2, G4int charge)
G4double GetPcm(const G4ThreeVector &p1, G4double m1, const G4ThreeVector &p2, G4double m2)
G4int FindPartner(const G4ThreeVector &p1, G4double m1, std::vector< std::pair< G4int, G4ThreeVector > > &neutron, G4double m2, G4int charge)
void PushDeuteron(const G4ThreeVector &p1, const G4ThreeVector &p2, G4int charge, G4ReactionProductVector *result)
void SetP0Coalescence(const G4HadProjectile &thePrimary, G4String)
const G4ParticleDefinition * GetDefinition() const
const G4LorentzVector & Get4Momentum() const
static G4Neutron * Neutron()
G4double GetPDGMass() const
G4int GetPDGEncoding() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
G4ParticleDefinition * FindAntiParticle(G4int PDGEncoding)
static G4int GetModelID(const G4int modelIndex)
static G4Proton * Proton()
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
void SetCreatorModelID(const G4int mod)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetMass(const G4double mas)