ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGNeutronInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGNeutronInelastic.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 "G4RPGNeutronInelastic.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "Randomize.hh"
32 
35  G4Nucleus& targetNucleus)
36 {
38  const G4HadProjectile* originalIncident = &aTrack;
39 
40  // create the target particle
41  G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
42 
43  G4ReactionProduct modifiedOriginal;
44  modifiedOriginal = *originalIncident;
45  G4ReactionProduct targetParticle;
46  targetParticle = *originalTarget;
47  if( originalIncident->GetKineticEnergy()/GeV < 0.01 + 2.*G4UniformRand()/9. )
48  {
49  SlowNeutron(originalIncident,modifiedOriginal,targetParticle,targetNucleus );
50  delete originalTarget;
51  return &theParticleChange;
52  }
53 
54  // Fermi motion and evaporation
55  // As of Geant3, the Fermi energy calculation had not been Done
56  G4double ek = originalIncident->GetKineticEnergy()/MeV;
57  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
58 
59  G4double tkin = targetNucleus.Cinema( ek );
60  ek += tkin;
61  modifiedOriginal.SetKineticEnergy( ek*MeV );
62  G4double et = ek + amas;
63  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
64  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
65  if( pp > 0.0 )
66  {
67  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
68  modifiedOriginal.SetMomentum( momentum * (p/pp) );
69  }
70  //
71  // calculate black track energies
72  //
73  tkin = targetNucleus.EvaporationEffects( ek );
74  ek -= tkin;
75  modifiedOriginal.SetKineticEnergy(ek);
76  et = ek + amas;
77  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
78  pp = modifiedOriginal.GetMomentum().mag();
79  if( pp > 0.0 )
80  {
81  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
82  modifiedOriginal.SetMomentum( momentum * (p/pp) );
83  }
84  const G4double cutOff = 0.1;
85  if( modifiedOriginal.GetKineticEnergy()/MeV <= cutOff )
86  {
87  SlowNeutron( originalIncident, modifiedOriginal, targetParticle, targetNucleus );
88  delete originalTarget;
89  return &theParticleChange;
90  }
91 
92  G4ReactionProduct currentParticle = modifiedOriginal;
93  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
94  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
95  G4bool incidentHasChanged = false;
96  G4bool targetHasChanged = false;
97  G4bool quasiElastic = false;
98  G4FastVector<G4ReactionProduct,256> vec; // vec will contain sec. particles
99  G4int vecLen = 0;
100  vec.Initialize( 0 );
101 
102  InitialCollision(vec, vecLen, currentParticle, targetParticle,
103  incidentHasChanged, targetHasChanged);
104 
105  CalculateMomenta(vec, vecLen,
106  originalIncident, originalTarget, modifiedOriginal,
107  targetNucleus, currentParticle, targetParticle,
108  incidentHasChanged, targetHasChanged, quasiElastic);
109 
110  SetUpChange(vec, vecLen,
111  currentParticle, targetParticle,
112  incidentHasChanged);
113 
114  delete originalTarget;
115  return &theParticleChange;
116 }
117 
118 void
120  G4ReactionProduct& modifiedOriginal,
121  G4ReactionProduct& targetParticle,
122  G4Nucleus& targetNucleus)
123 {
124  const G4double A = targetNucleus.GetA_asInt(); // atomic weight
125  const G4double Z = targetNucleus.GetZ_asInt(); // atomic number
126 
127  G4double currentKinetic = modifiedOriginal.GetKineticEnergy()/MeV;
128  G4double currentMass = modifiedOriginal.GetMass()/MeV;
129  if( A < 1.5 ) // Hydrogen
130  {
131  //
132  // very simple simulation of scattering angle and energy
133  // nonrelativistic approximation with isotropic angular
134  // distribution in the cms system
135  //
136  G4double cost1, eka = 0.0;
137  while (eka <= 0.0) { /* Loop checking, 01.09.2015, D.Wright */
138  cost1 = -1.0 + 2.0*G4UniformRand();
139  eka = 1.0 + 2.0*cost1*A + A*A;
140  }
141  G4double cost = std::min( 1.0, std::max( -1.0, (A*cost1+1.0)/std::sqrt(eka) ) );
142  eka /= (1.0+A)*(1.0+A);
143  G4double ek = currentKinetic*MeV/GeV;
144  G4double amas = currentMass*MeV/GeV;
145  ek *= eka;
146  G4double en = ek + amas;
147  G4double p = std::sqrt(std::abs(en*en-amas*amas));
148  G4double sint = std::sqrt(std::abs(1.0-cost*cost));
150  G4double px = sint*std::sin(phi);
151  G4double py = sint*std::cos(phi);
152  G4double pz = cost;
153  targetParticle.SetMomentum( px*GeV, py*GeV, pz*GeV );
154  G4double pxO = originalIncident->Get4Momentum().x()/GeV;
155  G4double pyO = originalIncident->Get4Momentum().y()/GeV;
156  G4double pzO = originalIncident->Get4Momentum().z()/GeV;
157  G4double ptO = pxO*pxO + pyO+pyO;
158  if( ptO > 0.0 )
159  {
160  G4double pO = std::sqrt(pxO*pxO+pyO*pyO+pzO*pzO);
161  cost = pzO/pO;
162  sint = 0.5*(std::sqrt(std::abs((1.0-cost)*(1.0+cost)))+std::sqrt(ptO)/pO);
163  G4double ph = pi/2.0;
164  if( pyO < 0.0 )ph = ph*1.5;
165  if( std::abs(pxO) > 0.000001 )ph = std::atan2(pyO,pxO);
166  G4double cosp = std::cos(ph);
167  G4double sinp = std::sin(ph);
168  px = cost*cosp*px - sinp*py+sint*cosp*pz;
169  py = cost*sinp*px + cosp*py+sint*sinp*pz;
170  pz = -sint*px + cost*pz;
171  }
172  else
173  {
174  if( pz < 0.0 )pz *= -1.0;
175  }
176  G4double pu = std::sqrt(px*px+py*py+pz*pz);
177  modifiedOriginal.SetMomentum( targetParticle.GetMomentum() * (p/pu) );
178  modifiedOriginal.SetKineticEnergy( ek*GeV );
179 
180  targetParticle.SetMomentum(
181  originalIncident->Get4Momentum().vect() - modifiedOriginal.GetMomentum() );
182  G4double pp = targetParticle.GetMomentum().mag();
183  G4double tarmas = targetParticle.GetMass();
184  targetParticle.SetTotalEnergy( std::sqrt( pp*pp + tarmas*tarmas ) );
185 
186  theParticleChange.SetEnergyChange( modifiedOriginal.GetKineticEnergy() );
188  pd->SetDefinition( targetParticle.GetDefinition() );
189  pd->SetMomentum( targetParticle.GetMomentum() );
191  return;
192  }
193 
194  G4FastVector<G4ReactionProduct,4> vec; // vec will contain the secondary particles
195  G4int vecLen = 0;
196  vec.Initialize( 0 );
197 
198  G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
199  G4double massVec[9];
200  massVec[0] = targetNucleus.AtomicMass( A+1.0, Z );
201  massVec[1] = theAtomicMass;
202  massVec[2] = 0.;
203  if (Z > 1.0) massVec[2] = targetNucleus.AtomicMass(A, Z-1.0);
204  massVec[3] = 0.;
205  if (Z > 1.0 && A > 1.0) massVec[3] = targetNucleus.AtomicMass(A-1.0, Z-1.0 );
206  massVec[4] = 0.;
207  if (Z > 1.0 && A > 2.0 && A-2.0 > Z-1.0)
208  massVec[4] = targetNucleus.AtomicMass( A-2.0, Z-1.0 );
209  massVec[5] = 0.;
210  if (Z > 2.0 && A > 3.0 && A-3.0 > Z-2.0)
211  massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-2.0 );
212  massVec[6] = 0.;
213  if (A > 1.0 && A-1.0 > Z) massVec[6] = targetNucleus.AtomicMass(A-1.0, Z);
214  massVec[7] = massVec[3];
215  massVec[8] = 0.;
216  if (Z > 2.0 && A > 1.0) massVec[8] = targetNucleus.AtomicMass( A-1.0,Z-2.0 );
217 
218  twoBody.NuclearReaction(vec, vecLen, originalIncident,
219  targetNucleus, theAtomicMass, massVec );
220 
223 
224  G4DynamicParticle* pd;
225  for( G4int i=0; i<vecLen; ++i ) {
226  pd = new G4DynamicParticle();
227  pd->SetDefinition( vec[i]->GetDefinition() );
228  pd->SetMomentum( vec[i]->GetMomentum() );
230  delete vec[i];
231  }
232 }
233 
234 
235 // Initial Collision
236 // selects the particle types arising from the initial collision of
237 // the neutron and target nucleon. Secondaries are assigned to
238 // forward and backward reaction hemispheres, but final state energies
239 // and momenta are not calculated here.
240 
241 void
243  G4int& vecLen,
244  G4ReactionProduct& currentParticle,
245  G4ReactionProduct& targetParticle,
246  G4bool& incidentHasChanged,
247  G4bool& targetHasChanged)
248 {
249  G4double KE = currentParticle.GetKineticEnergy()/GeV;
250 
251  G4int mult;
252  G4int partType;
253  std::vector<G4int> fsTypes;
254  G4int part1;
255  G4int part2;
256 
257  G4double testCharge;
258  G4double testBaryon;
259  G4double testStrange;
260 
261  // Get particle types according to incident and target types
262 
263  if (targetParticle.GetDefinition() == particleDef[neu]) {
264  mult = GetMultiplicityT1(KE);
265  fsTypes = GetFSPartTypesForNN(mult, KE);
266 
267  part1 = fsTypes[0];
268  part2 = fsTypes[1];
269  currentParticle.SetDefinition(particleDef[part1]);
270  targetParticle.SetDefinition(particleDef[part2]);
271  if (part1 == pro) {
272  if (part2 == neu) {
273  if (G4UniformRand() > 0.5) {
274  incidentHasChanged = true;
275  } else {
276  targetHasChanged = true;
277  currentParticle.SetDefinition(particleDef[part2]);
278  targetParticle.SetDefinition(particleDef[part1]);
279  }
280  } else {
281  targetHasChanged = true;
282  incidentHasChanged = true;
283  }
284 
285  } else { // neutron
286  if (part2 > neu && part2 < xi0) targetHasChanged = true;
287  }
288 
289  testCharge = 0.0;
290  testBaryon = 2.0;
291  testStrange = 0.0;
292 
293  } else { // target was a proton
294  mult = GetMultiplicityT0(KE);
295  fsTypes = GetFSPartTypesForNP(mult, KE);
296 
297  part1 = fsTypes[0];
298  part2 = fsTypes[1];
299  currentParticle.SetDefinition(particleDef[part1]);
300  targetParticle.SetDefinition(particleDef[part2]);
301  if (part1 == pro) {
302  if (part2 == pro) {
303  incidentHasChanged = true;
304  } else if (part2 == neu) {
305  if (G4UniformRand() > 0.5) {
306  incidentHasChanged = true;
307  targetHasChanged = true;
308  } else {
309  currentParticle.SetDefinition(particleDef[part2]);
310  targetParticle.SetDefinition(particleDef[part1]);
311  }
312 
313  } else if (part2 > neu && part2 < xi0) {
314  incidentHasChanged = true;
315  targetHasChanged = true;
316  }
317 
318  } else { // neutron
319  targetHasChanged = true;
320  }
321 
322  testCharge = 1.0;
323  testBaryon = 2.0;
324  testStrange = 0.0;
325  }
326 
327  // if (mult == 2 && !incidentHasChanged && !targetHasChanged)
328  // quasiElastic = true;
329 
330  // Remove incident and target from fsTypes
331 
332  fsTypes.erase(fsTypes.begin());
333  fsTypes.erase(fsTypes.begin());
334 
335  // Remaining particles are secondaries. Put them into vec.
336 
337  G4ReactionProduct* rp(0);
338  for(G4int i=0; i < mult-2; ++i ) {
339  partType = fsTypes[i];
340  rp = new G4ReactionProduct();
341  rp->SetDefinition(particleDef[partType]);
342  (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
343  vec.SetElement(vecLen++, rp);
344  }
345 
346  // Check conservation of charge, strangeness, baryon number
347 
348  CheckQnums(vec, vecLen, currentParticle, targetParticle,
349  testCharge, testBaryon, testStrange);
350 
351  return;
352 }