source: trunk/source/processes/hadronic/models/rpg/src/G4RPGProtonInelastic.cc@ 1288

Last change on this file since 1288 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 10.0 KB
Line 
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// $Id: G4RPGProtonInelastic.cc,v 1.4 2008/05/05 21:21:55 dennis Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29
30#include "G4RPGProtonInelastic.hh"
31#include "Randomize.hh"
32
33G4HadFinalState*
34G4RPGProtonInelastic::ApplyYourself(const G4HadProjectile& aTrack,
35 G4Nucleus& targetNucleus )
36{
37 theParticleChange.Clear();
38 const G4HadProjectile *originalIncident = &aTrack;
39 if (originalIncident->GetKineticEnergy()<= 0.1)
40 {
41 theParticleChange.SetStatusChange(isAlive);
42 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
43 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
44 return &theParticleChange;
45 }
46
47 //
48 // create the target particle
49 //
50 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
51
52 if (originalIncident->GetKineticEnergy()/GeV < 0.01+2.*G4UniformRand()/9. )
53 {
54 SlowProton( originalIncident, targetNucleus );
55 delete originalTarget;
56 return &theParticleChange;
57 }
58
59 // Fermi motion and evaporation
60 // As of Geant3, the Fermi energy calculation had not been Done
61
62 G4double ek = originalIncident->GetKineticEnergy();
63 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
64 G4ReactionProduct modifiedOriginal;
65 modifiedOriginal = *originalIncident;
66
67 G4double tkin = targetNucleus.Cinema( ek );
68 ek += tkin;
69 modifiedOriginal.SetKineticEnergy(ek);
70 G4double et = ek + amas;
71 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
72 G4double pp = modifiedOriginal.GetMomentum().mag();
73 if (pp > 0.0) {
74 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
75 modifiedOriginal.SetMomentum( momentum * (p/pp) );
76 }
77 //
78 // calculate black track energies
79 //
80 tkin = targetNucleus.EvaporationEffects(ek);
81 ek -= tkin;
82 modifiedOriginal.SetKineticEnergy(ek);
83 et = ek + amas;
84 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
85 pp = modifiedOriginal.GetMomentum().mag();
86 if (pp > 0.0) {
87 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
88 modifiedOriginal.SetMomentum( momentum * (p/pp) );
89 }
90 const G4double cutOff = 0.1;
91 if (modifiedOriginal.GetKineticEnergy() < cutOff) {
92 SlowProton( originalIncident, targetNucleus );
93 delete originalTarget;
94 return &theParticleChange;
95 }
96
97 G4ReactionProduct currentParticle = modifiedOriginal;
98 G4ReactionProduct targetParticle;
99 targetParticle = *originalTarget;
100 currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
101 targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
102 G4bool incidentHasChanged = false;
103 G4bool targetHasChanged = false;
104 G4bool quasiElastic = false;
105 G4FastVector<G4ReactionProduct,256> vec; // vec will contain the sec. particles
106 G4int vecLen = 0;
107 vec.Initialize( 0 );
108
109 InitialCollision(vec, vecLen, currentParticle, targetParticle,
110 incidentHasChanged, targetHasChanged);
111
112 CalculateMomenta(vec, vecLen,
113 originalIncident, originalTarget, modifiedOriginal,
114 targetNucleus, currentParticle, targetParticle,
115 incidentHasChanged, targetHasChanged, quasiElastic);
116
117 SetUpChange( vec, vecLen,
118 currentParticle, targetParticle,
119 incidentHasChanged );
120
121 delete originalTarget;
122 return &theParticleChange;
123}
124
125
126void
127G4RPGProtonInelastic::SlowProton(const G4HadProjectile *originalIncident,
128 G4Nucleus &targetNucleus )
129{
130 const G4double A = targetNucleus.GetN(); // atomic weight
131 const G4double Z = targetNucleus.GetZ(); // atomic number
132// G4double currentKinetic = originalIncident->GetKineticEnergy();
133 //
134 // calculate Q-value of reactions
135 //
136 G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
137 G4double massVec[9];
138 massVec[0] = targetNucleus.AtomicMass( A+1.0, Z+1.0 );
139 massVec[1] = 0.;
140 if (A > Z+1.0)
141 massVec[1] = targetNucleus.AtomicMass( A , Z+1.0 );
142 massVec[2] = theAtomicMass;
143 massVec[3] = 0.;
144 if (A > 1.0 && A-1.0 > Z)
145 massVec[3] = targetNucleus.AtomicMass( A-1.0, Z );
146 massVec[4] = 0.;
147 if (A > 2.0 && A-2.0 > Z)
148 massVec[4] = targetNucleus.AtomicMass( A-2.0, Z );
149 massVec[5] = 0.;
150 if (A > 3.0 && Z > 1.0 && A-3.0 > Z-1.0)
151 massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-1.0 );
152 massVec[6] = 0.;
153 if (A > 1.0 && A-1.0 > Z+1.0)
154 massVec[6] = targetNucleus.AtomicMass( A-1.0, Z+1.0 );
155 massVec[7] = massVec[3];
156 massVec[8] = 0.;
157 if (A > 1.0 && Z > 1.0)
158 massVec[8] = targetNucleus.AtomicMass( A-1.0, Z-1.0 );
159
160 G4FastVector<G4ReactionProduct,4> vec; // vec will contain the secondary particles
161 G4int vecLen = 0;
162 vec.Initialize( 0 );
163
164 twoBody.NuclearReaction( vec, vecLen, originalIncident,
165 targetNucleus, theAtomicMass, massVec );
166
167 theParticleChange.SetStatusChange( stopAndKill );
168 theParticleChange.SetEnergyChange( 0.0 );
169
170 G4DynamicParticle *pd;
171 for( G4int i=0; i<vecLen; ++i )
172 {
173 pd = new G4DynamicParticle();
174 pd->SetDefinition( vec[i]->GetDefinition() );
175 pd->SetMomentum( vec[i]->GetMomentum() );
176 theParticleChange.AddSecondary( pd );
177 delete vec[i];
178 }
179}
180
181
182// Initial Collision
183// selects the particle types arising from the initial collision of
184// the proton and target nucleon. Secondaries are assigned to forward
185// and backward reaction hemispheres, but final state energies and
186// momenta are not calculated here.
187
188void
189G4RPGProtonInelastic::InitialCollision(G4FastVector<G4ReactionProduct,256>& vec,
190 G4int& vecLen,
191 G4ReactionProduct& currentParticle,
192 G4ReactionProduct& targetParticle,
193 G4bool& incidentHasChanged,
194 G4bool& targetHasChanged)
195{
196 G4double KE = currentParticle.GetKineticEnergy()/GeV;
197
198 G4int mult;
199 G4int partType;
200 std::vector<G4int> fsTypes;
201 G4int part1;
202 G4int part2;
203
204 G4double testCharge;
205 G4double testBaryon;
206 G4double testStrange;
207
208 // Get particle types according to incident and target types
209
210 if (targetParticle.GetDefinition() == particleDef[pro]) {
211 mult = GetMultiplicityT1(KE);
212 fsTypes = GetFSPartTypesForPP(mult, KE);
213
214 part1 = fsTypes[0];
215 part2 = fsTypes[1];
216 currentParticle.SetDefinition(particleDef[part1]);
217 targetParticle.SetDefinition(particleDef[part2]);
218 if (part1 == pro) {
219 if (part2 == neu) {
220 if (G4UniformRand() > 0.5) {
221 incidentHasChanged = true;
222 targetParticle.SetDefinition(particleDef[part1]);
223 currentParticle.SetDefinition(particleDef[part2]);
224 } else {
225 targetHasChanged = true;
226 }
227 } else if (part2 > neu && part2 < xi0) {
228 targetHasChanged = true;
229 }
230
231 } else { // neutron
232 targetHasChanged = true;
233 incidentHasChanged = true;
234 }
235
236 testCharge = 2.0;
237 testBaryon = 2.0;
238 testStrange = 0.0;
239
240 } else { // target was a neutron
241 mult = GetMultiplicityT0(KE);
242 fsTypes = GetFSPartTypesForPN(mult, KE);
243
244 part1 = fsTypes[0];
245 part2 = fsTypes[1];
246 currentParticle.SetDefinition(particleDef[part1]);
247 targetParticle.SetDefinition(particleDef[part2]);
248 if (part1 == pro) {
249 if (part2 == pro) {
250 targetHasChanged = true;
251 } else if (part2 == neu) {
252 if (G4UniformRand() > 0.5) {
253 incidentHasChanged = true;
254 targetHasChanged = true;
255 targetParticle.SetDefinition(particleDef[part1]);
256 currentParticle.SetDefinition(particleDef[part2]);
257 }
258 } else { // hyperon
259 targetHasChanged = true;
260 }
261
262 } else { // neutron
263 incidentHasChanged = true;
264 if (part2 > neu && part2 < xi0) targetHasChanged = true;
265 }
266
267 testCharge = 1.0;
268 testBaryon = 2.0;
269 testStrange = 0.0;
270 }
271
272 // Remove incident and target from fsTypes
273
274 fsTypes.erase(fsTypes.begin());
275 fsTypes.erase(fsTypes.begin());
276
277 // Remaining particles are secondaries. Put them into vec.
278
279 G4ReactionProduct* rp(0);
280 for(G4int i=0; i < mult-2; ++i ) {
281 partType = fsTypes[i];
282 rp = new G4ReactionProduct();
283 rp->SetDefinition(particleDef[partType]);
284 (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
285 vec.SetElement(vecLen++, rp);
286 }
287
288 // Check conservation of charge, strangeness, baryon number
289
290 CheckQnums(vec, vecLen, currentParticle, targetParticle,
291 testCharge, testBaryon, testStrange);
292
293 return;
294}
Note: See TracBrowser for help on using the repository browser.