| Geant4 Cross Reference |
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"
79 #include "G4PhysicsModelCatalog.hh"
80
81
82 G4CRCoalescence::G4CRCoalescence() : G4HadronicInteraction("G4CRCoalescence" ),
83 fP0_d( 0.0 ), fP0_dbar( 0.0 ), secID( -1 ) {
84 secID = G4PhysicsModelCatalog::GetModelID( "model_G4CRCoalescence" );
85 }
86
87
88 G4CRCoalescence::~G4CRCoalescence() {}
89
90
91 void 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
113 void G4CRCoalescence::GenerateDeuterons( G4ReactionProductVector* result ) {
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
162 G4ParticleDefinition* prt = G4ParticleTable::GetParticleTable()->FindParticle( "proton" );
163 G4ReactionProduct* finalp = new G4ReactionProduct;
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
180 G4ParticleDefinition* nrt = G4ParticleTable::GetParticleTable()->FindParticle( "neutron" );
181 G4ReactionProduct* finaln = new G4ReactionProduct;
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
198 G4ParticleDefinition* pbar = G4ParticleTable::GetParticleTable()->FindAntiParticle( "proton" );
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
216 G4ParticleDefinition* nbar = G4ParticleTable::GetParticleTable()->FindAntiParticle( "neutron" );
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
230 void 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 ) {
235 G4ParticleDefinition* deuteron = G4ParticleTable::GetParticleTable()->FindParticle( "deuteron" );
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 {
247 G4ParticleDefinition* antideuteron = G4ParticleTable::GetParticleTable()->FindAntiParticle( "deuteron" );
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
262 G4int G4CRCoalescence::FindPartner( const G4ThreeVector &p1, G4double m1,
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
279 G4bool G4CRCoalescence::Coalescence( const G4ThreeVector &p1, G4double m1,
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
287 G4bool G4CRCoalescence::Coalescence( G4double p1x, G4double p1y, G4double p1z, G4double m1,
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
298 G4double G4CRCoalescence::GetPcm( const G4ThreeVector& p1, G4double m1,
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
305 G4double G4CRCoalescence::GetPcm( G4double p1x, G4double p1y, G4double p1z, G4double m1,
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
313 G4double G4CRCoalescence::GetS( G4double p1x, G4double p1y, G4double p1z, G4double m1,
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 }
320