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