Geant4-11
G4CRCoalescence.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27//---------------------------------------------------------------------------
28//
29// ClassName: G4CRCoalescence ("CR" stands for "Cosmic Ray")
30//
31// Author: 2020 Alberto Ribon , based on code written by
32// Diego Mauricio Gomez Coral for the GAPS Collaboration
33//
34// Description: This class can be optionally used in the method:
35//
36// G4TheoFSGenerator::ApplyYourself
37//
38// to coalesce pairs of proton-neutron and antiproton-antineutron
39// into deuterons and antideuterons, respectively, from the list
40// of secondaries produced by a string model.
41// This class can be useful in particular for Cosmic Ray (CR)
42// applications.
43// By default, this class is not used.
44// However, it can be enabled via the UI command:
45//
46// /process/had/enableCRCoalescence true
47//
48// It is assumed that the candidate proton-neutron and
49// antiproton-antideuteron pairs originate from the same
50// spatial position, so the condition for coalescence takes
51// into account only their closeness in momentum space.
52//
53// This class is based entirely on code written by
54// Diego Mauricio Gomez Coral for the GAPS Collaboration.
55// The main application of this work is for cosmic ray physics.
56//
57// Notes:
58// - In its current version, coalescence can occur only for
59// proton projectile (because the coalescence parameters
60// for deuteron and antideuteron are set to non-null values
61// only for the case of proton projectile).
62// - This class is not meant be used for secondaries produces
63// by intranuclear cascade models - such as BERT, BIC and
64// INCL - which should have already a coalescence phase.
65//
66// Modified:
67//
68//----------------------------------------------------------------------------
69//
70#include "G4DynamicParticle.hh"
71#include "G4Proton.hh"
72#include "G4Neutron.hh"
73#include "G4Deuteron.hh"
74#include "G4AntiProton.hh"
75#include "G4AntiNeutron.hh"
76#include "G4CRCoalescence.hh"
77#include "G4ReactionProduct.hh"
78#include "G4IonTable.hh"
80
81
83 fP0_d( 0.0 ), fP0_dbar( 0.0 ), secID( -1 ) {
84 secID = G4PhysicsModelCatalog::GetModelID( "model_G4CRCoalescence" );
85}
86
87
89
90
91void G4CRCoalescence::SetP0Coalescence( const G4HadProjectile &thePrimary, G4String /* model */ ) {
92 // Note by A.R. : in the present version, the coalescence parameters are set only for
93 // proton projectile. If we want to extend this coalescence algorithm
94 // for other applications, besides cosmic rays, we need to set these
95 // coalescence parameters also for all projectiles.
96 // (Note that the method "GenerateDeuterons", instead, can be already used
97 // as it is for all projectiles.)
98 fP0_dbar = 0.0;
99 fP0_d = 0.0;
100 if ( thePrimary.GetDefinition()->GetPDGEncoding() == 2212 ) { // proton
101 G4double mproj = thePrimary.GetDefinition()->GetPDGMass();
102 G4double pz = thePrimary.Get4Momentum().z();
103 G4double ekin = std::sqrt( pz*pz + mproj*mproj ) - mproj;
104 if ( ekin > 10.0 ) {
105 fP0_dbar = 130.0 / ( 1.0 + std::exp( 21.6 - std::log( 0.001*ekin )/0.089 ) ); // set p0 for antideuteron
106 fP0_d = 118.1 * ( 1.0 + std::exp( 5.53 - std::log( 0.001*ekin )/0.43 ) ); // set p0 for deuteron
107 }
108 }
109 //G4cout << "Coalescence parameter p0 deuteron / antideuteron: " << fP0_d << " / " << fP0_dbar << G4endl;
110}
111
112
114 // Deuteron clusters are made with the first nucleon pair that fulfills
115 // the coalescence conditions, starting with the protons.
116 // A deuteron is a pair (i,j) where i is the proton and j the neutron in current event
117 // with the relative momentum less than p0 (i.e. within a sphere of radius p0).
118 // The same applies for antideuteron clusters, with antiprotons and antineutrons,
119 // instead of protons and neutrons, respectively.
120
121 // Vectors of index-position and 3-momentum pairs for, respectively:
122 // protons, neutrons, antiprotons and antineutrons
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 ) { // proton
130 proton.push_back( std::make_pair( i, result->operator[](i)->GetMomentum() ) );
131 result->erase( result->begin() + i );
132 }
133 }
134 for ( unsigned int i = 0; i < result->size(); ++i ) {
135 G4int pdgid = result->operator[](i)->GetDefinition()->GetPDGEncoding();
136 if ( pdgid == 2112 ) { // neutron
137 neutron.push_back( std::make_pair( i, result->operator[](i)->GetMomentum() ) );
138 result->erase( result->begin() + i );
139 }
140 }
141 for ( unsigned int i = 0; i < result->size(); ++i ) {
142 G4int pdgid = result->operator[](i)->GetDefinition()->GetPDGEncoding();
143 if ( pdgid == -2212 ) { // antiproton
144 antiproton.push_back( std::make_pair( i, result->operator[](i)->GetMomentum() ) );
145 result->erase( result->begin() + i );
146 }
147 }
148 for ( unsigned int i = 0; i < result->size(); ++i ) {
149 G4int pdgid = result->operator[](i)->GetDefinition()->GetPDGEncoding();
150 if ( pdgid == -2112 ) { // antineutron
151 antineutron.push_back( std::make_pair( i, result->operator[](i)->GetMomentum() ) );
152 result->erase( result->begin() + i );
153 }
154 }
155
156 for ( unsigned int i = 0; i < proton.size(); ++i ) { // loop over protons
157 if ( proton.at(i).first == -1 ) continue;
158 G4ThreeVector p1 = proton.at(i).second;
159 int partner1 = FindPartner( p1, G4Proton::Proton()->GetPDGMass(), neutron,
160 G4Neutron::Neutron()->GetPDGMass(), 1 );
161 if ( partner1 == -1 ) { // if no partner found, then the proton is a final-state secondary
164 finalp->SetDefinition( prt );
165 G4double massp = prt->GetPDGMass();
166 G4double totalEnergy = std::sqrt( p1.mag()*p1.mag() + massp*massp );
167 finalp->SetMomentum( p1 );
168 finalp->SetTotalEnergy( totalEnergy );
169 finalp->SetMass( massp );
170 result->push_back( finalp );
171 continue;
172 }
173 G4ThreeVector p2 = neutron.at(partner1).second;
174 PushDeuteron( p1, p2, 1, result );
175 neutron.at(partner1).first = -1; // tag the bound neutron
176 }
177
178 for ( unsigned int i = 0; i < neutron.size(); ++i ) { // loop over neutrons
179 if ( neutron.at(i).first == -1 ) continue; // Skip already bound neutron, else it is a final-state secondary
182 finaln->SetDefinition( nrt );
183 G4ThreeVector p2 = neutron.at(i).second;
184 G4double massn = nrt->GetPDGMass();
185 G4double totalEnergy = std::sqrt( p2.mag()*p2.mag() + massn*massn );
186 finaln->SetMomentum( p2 );
187 finaln->SetTotalEnergy( totalEnergy );
188 finaln->SetMass( massn );
189 result->push_back( finaln );
190 }
191
192 for ( unsigned int i = 0; i < antiproton.size(); ++i ) { // loop over antiprotons
193 if ( antiproton.at(i).first == -1 ) continue;
194 G4ThreeVector p1 = antiproton.at(i).second;
195 int partner1 = FindPartner( p1, G4Proton::Proton()->GetPDGMass(), antineutron,
196 G4Neutron::Neutron()->GetPDGMass(), -1 );
197 if ( partner1 == -1 ) { // if no partner found, then the antiproton is a final-state secondary
199 G4ReactionProduct* finalpbar = new G4ReactionProduct;
200 finalpbar->SetDefinition( pbar );
201 G4double massp = pbar->GetPDGMass();
202 G4double totalEnergy = std::sqrt( p1.mag()*p1.mag() + massp*massp );
203 finalpbar->SetMomentum( p1 );
204 finalpbar->SetTotalEnergy( totalEnergy );
205 finalpbar->SetMass( massp );
206 result->push_back( finalpbar );
207 continue;
208 }
209 G4ThreeVector p2 = antineutron.at(partner1).second;
210 PushDeuteron( p1, p2, -1, result );
211 antineutron.at(partner1).first = -1; // tag the bound antineutron
212 }
213
214 for ( unsigned int i = 0; i < antineutron.size(); ++i ) { // loop over antineutrons
215 if ( antineutron.at(i).first == -1 ) continue; // Skip already bound antineutron, else it is a final-state secondary
217 G4ReactionProduct* finalnbar = new G4ReactionProduct;
218 finalnbar->SetDefinition( nbar );
219 G4ThreeVector p2 = antineutron.at(i).second;
220 G4double massn = nbar->GetPDGMass();
221 G4double totalEnergy = std::sqrt( p2.mag()*p2.mag() + massn*massn );
222 finalnbar->SetMomentum( p2 );
223 finalnbar->SetTotalEnergy( totalEnergy );
224 finalnbar->SetMass( massn );
225 result->push_back( finalnbar );
226 }
227}
228
229
230void G4CRCoalescence::PushDeuteron( const G4ThreeVector &p1, const G4ThreeVector &p2, G4int charge, // input
231 G4ReactionProductVector* result ) { // output
232 // Create a deuteron or antideuteron (depending on "charge") object (of type G4ReactionProduct)
233 // from the two input momenta "p1" and "p2", and push it to the vector "result".
234 if ( charge > 0 ) {
236 G4ReactionProduct* finaldeut = new G4ReactionProduct;
237 finaldeut->SetDefinition( deuteron );
238 G4ThreeVector psum = p1 + p2;
239 G4double massd = deuteron->GetPDGMass();
240 G4double totalEnergy = std::sqrt( psum.mag()*psum.mag() + massd*massd );
241 finaldeut->SetMomentum( psum );
242 finaldeut->SetTotalEnergy( totalEnergy );
243 finaldeut->SetMass( massd );
244 finaldeut->SetCreatorModelID( secID );
245 result->push_back( finaldeut );
246 } else {
248 G4ReactionProduct* finalantideut = new G4ReactionProduct;
249 finalantideut->SetDefinition( antideuteron );
250 G4ThreeVector psum = p1 + p2;
251 G4double massd = antideuteron->GetPDGMass();
252 G4double totalEnergy = std::sqrt( psum.mag()*psum.mag() + massd*massd );
253 finalantideut->SetMomentum( psum );
254 finalantideut->SetTotalEnergy( totalEnergy );
255 finalantideut->SetMass( massd );
256 finalantideut->SetCreatorModelID( secID );
257 result->push_back( finalantideut );
258 }
259}
260
261
263 std::vector< std::pair< G4int, G4ThreeVector > > &neutron,
264 G4double m2, G4int charge ) {
265 // Find a nucleon/antinucleon (depending on "charge") partner, from the input list "neutron"
266 // (which is a vector of either neutron or antineutron particles depending on "charge")
267 // within a sphere of radius p0 centered at the input momentum "p1"; exclude already bound
268 // particles (neutrons or antineutrons depending on "charge") of "neutron".
269 for ( unsigned int j = 0; j < neutron.size(); ++j ) {
270 if ( neutron.at(j).first == -1 ) continue; // skip already bound particle
271 G4ThreeVector p2 = neutron.at(j).second;
272 if ( ! Coalescence( p1, m1, p2, m2, charge ) ) continue;
273 return j;
274 }
275 return -1; // no partner found
276}
277
278
280 const G4ThreeVector &p2, G4double m2, G4int charge ) {
281 // Returns true if the momenta of the two nucleons/antinucleons (depending on "charge") are
282 // inside of an sphere of radius p0 (assuming that the two particles are in the same spatial place).
283 return Coalescence( p1.x(), p1.y(), p1.z(), m1, p2.x(), p2.y(), p2.z(), m2, charge );
284}
285
286
288 G4double p2x, G4double p2y, G4double p2z, G4double m2,
289 G4int charge ) {
290 // Returns true if the momenta of the two nucleons/antinucleons (depending on "charge") are
291 // inside of a sphere of radius p0 (assuming that the two particles are in the same spatial place).
292 G4double deltaP = GetPcm( p1x, p1y, p1z, m1, p2x, p2y, p2z, m2 );
293 if ( charge > 0 ) return ( deltaP < fP0_d );
294 else return ( deltaP < fP0_dbar );
295}
296
297
299 const G4ThreeVector& p2, G4double m2 ) {
300 // Momentum in the center-of-mass frame of two particles from LAB values.
301 return GetPcm( p1.x(), p1.y(), p1.z(), m1, p2.x(), p2.y(), p2.z(), m2 );
302}
303
304
306 G4double p2x, G4double p2y, G4double p2z, G4double m2 ) {
307 // Momentum in the center-of-mass frame of two particles from LAB values.
308 G4double scm = GetS( p1x, p1y, p1z, m1, p2x, p2y, p2z, m2 );
309 return std::sqrt( (scm - (m1-m2)*(m1-m2))*(scm - (m1+m2)*(m1+m2)) ) / (2.0*std::sqrt( scm ));
310}
311
312
314 G4double p2x, G4double p2y, G4double p2z, G4double m2 ) {
315 // Square of center-of-mass energy of two particles from LAB values.
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);
319}
std::vector< G4ReactionProduct * > G4ReactionProductVector
static constexpr double m2
Definition: G4SIunits.hh:110
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
double z() const
double x() const
double y() const
double mag() const
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()
Definition: G4Neutron.cc:103
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
G4ParticleDefinition * FindAntiParticle(G4int PDGEncoding)
static G4int GetModelID(const G4int modelIndex)
static G4Proton * Proton()
Definition: G4Proton.cc:92
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)