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