| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //--------------------------------------------
28 //
29 // ClassName: G4CRCoalescence ("CR" stand
30 //
31 // Author: 2020 Alberto Ribon , based on
32 // Diego Mauricio Gomez Coral fo
33 //
34 // Description: This class can be optionally
35 //
36 // G4TheoFSGenerator::ApplyYou
37 //
38 // to coalesce pairs of proton-n
39 // into deuterons and antideuter
40 // of secondaries produced by a
41 // This class can be useful in p
42 // applications.
43 // By default, this class is not
44 // However, it can be enabled vi
45 //
46 // /process/had/enableCRCoales
47 //
48 // It is assumed that the candid
49 // antiproton-antideuteron pairs
50 // spatial position, so the cond
51 // into account only their close
52 //
53 // This class is based entirely
54 // Diego Mauricio Gomez Coral fo
55 // The main application of this
56 //
57 // Notes:
58 // - In its current version, co
59 // proton projectile (because
60 // for deuteron and antideute
61 // only for the case of proto
62 // - This class is not meant be
63 // by intranuclear cascade mo
64 // INCL - which should have a
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() : G4Hadroni
83 fP0_d( 0.0 ), fP0_dbar( 0.0 ), secID( -1 ) {
84 secID = G4PhysicsModelCatalog::GetModelID( "
85 }
86
87
88 G4CRCoalescence::~G4CRCoalescence() {}
89
90
91 void G4CRCoalescence::SetP0Coalescence( const
92 // Note by A.R. : in the present version, th
93 // proton projectile. If we w
94 // for other applications, be
95 // coalescence parameters als
96 // (Note that the method "Gen
97 // as it is for all projecti
98 fP0_dbar = 0.0;
99 fP0_d = 0.0;
100 if ( thePrimary.GetDefinition()->GetPDGEncod
101 G4double mproj = thePrimary.GetDefinition(
102 G4double pz = thePrimary.Get4Momentum().z(
103 G4double ekin = std::sqrt( pz*pz + mproj*m
104 if ( ekin > 10.0 ) {
105 fP0_dbar = 130.0 / ( 1.0 + std::exp( 21.
106 fP0_d = 118.1 * ( 1.0 + std::exp( 5.5
107 }
108 }
109 //G4cout << "Coalescence parameter p0 deuter
110 }
111
112
113 void G4CRCoalescence::GenerateDeuterons( G4Rea
114 // Deuteron clusters are made with the first
115 // the coalescence conditions, starting with
116 // A deuteron is a pair (i,j) where i is the
117 // with the relative momentum less than p0 (
118 // The same applies for antideuteron cluster
119 // instead of protons and neutrons, respecti
120
121 // Vectors of index-position and 3-momentum
122 // protons, neutrons, antiprotons and antine
123 std::vector< std::pair< G4int, G4ThreeVector
124 std::vector< std::pair< G4int, G4ThreeVector
125 std::vector< std::pair< G4int, G4ThreeVector
126 std::vector< std::pair< G4int, G4ThreeVector
127 for ( unsigned int i = 0; i < result->size()
128 G4int pdgid = result->operator[](i)->GetDe
129 if ( pdgid == 2212 ) { // proton
130 proton.push_back( std::make_pair( i, res
131 result->erase( result->begin() + i );
132 }
133 }
134 for ( unsigned int i = 0; i < result->size()
135 G4int pdgid = result->operator[](i)->GetDe
136 if ( pdgid == 2112 ) { // neutron
137 neutron.push_back( std::make_pair( i, re
138 result->erase( result->begin() + i );
139 }
140 }
141 for ( unsigned int i = 0; i < result->size()
142 G4int pdgid = result->operator[](i)->GetDe
143 if ( pdgid == -2212 ) { // antiproton
144 antiproton.push_back( std::make_pair( i,
145 result->erase( result->begin() + i );
146 }
147 }
148 for ( unsigned int i = 0; i < result->size()
149 G4int pdgid = result->operator[](i)->GetDe
150 if ( pdgid == -2112 ) { // antineutron
151 antineutron.push_back( std::make_pair( i
152 result->erase( result->begin() + i );
153 }
154 }
155
156 for ( unsigned int i = 0; i < proton.size();
157 if ( proton.at(i).first == -1 ) continue;
158 G4ThreeVector p1 = proton.at(i).second;
159 int partner1 = FindPartner( p1, G4Proton::
160 G4Neutron::Neutron()->GetPDGMass(), 1
161 if ( partner1 == -1 ) { // if no partner
162 G4ParticleDefinition* prt = G4ParticleTa
163 G4ReactionProduct* finalp = new G4Reacti
164 finalp->SetDefinition( prt );
165 G4double massp = prt->GetPDGMass();
166 G4double totalEnergy = std::sqrt( p1.mag
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).se
174 PushDeuteron( p1, p2, 1, result );
175 neutron.at(partner1).first = -1; // tag t
176 }
177
178 for ( unsigned int i = 0; i < neutron.size()
179 if ( neutron.at(i).first == -1 ) continue;
180 G4ParticleDefinition* nrt = G4ParticleTabl
181 G4ReactionProduct* finaln = new G4Reaction
182 finaln->SetDefinition( nrt );
183 G4ThreeVector p2 = neutron.at(i).second;
184 G4double massn = nrt->GetPDGMass();
185 G4double totalEnergy = std::sqrt( p2.mag()
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.siz
193 if ( antiproton.at(i).first == -1 ) contin
194 G4ThreeVector p1 = antiproton.at(i).second
195 int partner1 = FindPartner( p1, G4Proton::
196 G4Neutron::Neutron()->GetPDGMass(), -1
197 if ( partner1 == -1 ) { // if no partner
198 G4ParticleDefinition* pbar = G4ParticleT
199 G4ReactionProduct* finalpbar = new G4Rea
200 finalpbar->SetDefinition( pbar );
201 G4double massp = pbar->GetPDGMass();
202 G4double totalEnergy = std::sqrt( p1.mag
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
210 PushDeuteron( p1, p2, -1, result );
211 antineutron.at(partner1).first = -1; // t
212 }
213
214 for ( unsigned int i = 0; i < antineutron.si
215 if ( antineutron.at(i).first == -1 ) conti
216 G4ParticleDefinition* nbar = G4ParticleTab
217 G4ReactionProduct* finalnbar = new G4React
218 finalnbar->SetDefinition( nbar );
219 G4ThreeVector p2 = antineutron.at(i).secon
220 G4double massn = nbar->GetPDGMass();
221 G4double totalEnergy = std::sqrt( p2.mag()
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 G4Th
231 G4ReactionProductVector* result )
232 // Create a deuteron or antideuteron (depend
233 // from the two input momenta "p1" and "p2",
234 if ( charge > 0 ) {
235 G4ParticleDefinition* deuteron = G4Particl
236 G4ReactionProduct* finaldeut = new G4React
237 finaldeut->SetDefinition( deuteron );
238 G4ThreeVector psum = p1 + p2;
239 G4double massd = deuteron->GetPDGMass();
240 G4double totalEnergy = std::sqrt( psum.mag
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 = G4Par
248 G4ReactionProduct* finalantideut = new G4R
249 finalantideut->SetDefinition( antideuteron
250 G4ThreeVector psum = p1 + p2;
251 G4double massd = antideuteron->GetPDGMass(
252 G4double totalEnergy = std::sqrt( psum.mag
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 G4Th
263 std::vector< std::pair< G4int, G4T
264 G4double m2, G4int charge ) {
265 // Find a nucleon/antinucleon (depending on
266 // (which is a vector of either neutron or a
267 // within a sphere of radius p0 centered at
268 // particles (neutrons or antineutrons depen
269 for ( unsigned int j = 0; j < neutron.size()
270 if ( neutron.at(j).first == -1 ) continue;
271 G4ThreeVector p2 = neutron.at(j).second;
272 if ( ! Coalescence( p1, m1, p2, m2, charge
273 return j;
274 }
275 return -1; // no partner found
276 }
277
278
279 G4bool G4CRCoalescence::Coalescence( const G4T
280 const G4ThreeVector &p2, G4double
281 // Returns true if the momenta of the two nu
282 // inside of an sphere of radius p0 (assumin
283 return Coalescence( p1.x(), p1.y(), p1.z(),
284 }
285
286
287 G4bool G4CRCoalescence::Coalescence( G4double
288 G4double
289 G4int charge ) {
290 // Returns true if the momenta of the two nu
291 // inside of a sphere of radius p0 (assuming
292 G4double deltaP = GetPcm( p1x, p1y, p1z, m1,
293 if ( charge > 0 ) return ( deltaP < fP0_d );
294 else return ( deltaP < fP0_dbar
295 }
296
297
298 G4double G4CRCoalescence::GetPcm( const G4Thre
299 const G4ThreeVector& p2, G4double m2
300 // Momentum in the center-of-mass frame of t
301 return GetPcm( p1.x(), p1.y(), p1.z(), m1, p
302 }
303
304
305 G4double G4CRCoalescence::GetPcm( G4double p1x
306 G4double p2x, G4double p2y, G4double
307 // Momentum in the center-of-mass frame of t
308 G4double scm = GetS( p1x, p1y, p1z, m1, p2x,
309 return std::sqrt( (scm - (m1-m2)*(m1-m2))*(s
310 }
311
312
313 G4double G4CRCoalescence::GetS( G4double p1x,
314 G4double p2x, G4double p2y, G4double p
315 // Square of center-of-mass energy of two pa
316 G4double E1 = std::sqrt( p1x*p1x + p1y*p1y +
317 G4double E2 = std::sqrt( p2x*p2x + p2y*p2y +
318 return (E1+E2)*(E1+E2) - (p1x+p2x)*(p1x+p2x)
319 }
320