ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGTwoBody.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGTwoBody.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 <iostream>
29 #include <signal.h>
30 
31 #include "G4RPGTwoBody.hh"
32 #include "G4Log.hh"
33 #include "Randomize.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4Poisson.hh"
38 
40  : G4RPGReaction() {}
41 
42 
44 ReactionStage(const G4HadProjectile* /*originalIncident*/,
45  G4ReactionProduct& modifiedOriginal,
46  G4bool& /*incidentHasChanged*/,
47  const G4DynamicParticle* originalTarget,
48  G4ReactionProduct& targetParticle,
49  G4bool& /*targetHasChanged*/,
50  const G4Nucleus& targetNucleus,
51  G4ReactionProduct& currentParticle,
53  G4int& vecLen,
54  G4bool /*leadFlag*/,
55  G4ReactionProduct& /*leadingStrangeParticle*/)
56 {
57  //
58  // Derived from H. Fesefeldt's original FORTRAN code TWOB
59  //
60  // Generation of momenta for elastic and quasi-elastic 2 body reactions
61  //
62  // The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used.
63  // The b values are parametrizations from experimental data.
64  // Unavailable values are taken from those of similar reactions.
65  //
66 
67  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
68  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
69  const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
70  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
71  G4double currentMass = currentParticle.GetMass()/GeV;
72  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
73 
74  targetMass = targetParticle.GetMass()/GeV;
75  const G4double atomicWeight = targetNucleus.GetA_asInt();
76 
77  G4double etCurrent = currentParticle.GetTotalEnergy()/GeV;
78  G4double pCurrent = currentParticle.GetTotalMomentum()/GeV;
79 
80  G4double cmEnergy = std::sqrt( currentMass*currentMass +
81  targetMass*targetMass +
82  2.0*targetMass*etCurrent ); // in GeV
83 
84  if (cmEnergy < 0.01) { // 2-body scattering not possible
85  targetParticle.SetMass( 0.0 ); // flag that the target particle doesn't exist
86 
87  } else {
88  // Projectile momentum in cm
89 
90  G4double pf = targetMass*pCurrent/cmEnergy;
91 
92  //
93  // Set beam and target in centre of mass system
94  //
95  G4ReactionProduct pseudoParticle[3];
96 
97  if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
98  targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
99 
100  // G4double pf1 = pOriginal*mOriginal/std::sqrt(2.*mOriginal*(mOriginal+etOriginal));
101 
102  pseudoParticle[0].SetMass( targetMass*GeV );
103  pseudoParticle[0].SetTotalEnergy( etOriginal*GeV );
104  pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
105 
106  pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
107  pseudoParticle[1].SetMass( mOriginal*GeV );
108  pseudoParticle[1].SetKineticEnergy( 0.0 );
109 
110  } else {
111  pseudoParticle[0].SetMass( currentMass*GeV );
112  pseudoParticle[0].SetTotalEnergy( etCurrent*GeV );
113  pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*GeV );
114 
115  pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
116  pseudoParticle[1].SetMass( targetMass*GeV );
117  pseudoParticle[1].SetKineticEnergy( 0.0 );
118  }
119  //
120  // Transform into center of mass system
121  //
122  pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
123  pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
124  pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
125  //
126  // Set final state masses and energies in centre of mass system
127  //
128  currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
129  targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
130 
131  //
132  // Calculate slope b for elastic scattering on proton/neutron
133  //
134  const G4double cb = 0.01;
135  const G4double b1 = 4.225;
136  const G4double b2 = 1.795;
137  G4double b = std::max( cb, b1+b2*G4Log(pOriginal) );
138 
139  //
140  // Get cm scattering angle by sampling t from tmin to tmax
141  //
142  G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/GeV;
143  G4double exindt = G4Exp(-btrang) - 1.0;
144  G4double costheta = 1.0 + 2*G4Log( 1.0+G4UniformRand()*exindt ) /btrang;
145  costheta = std::max(-1., std::min(1., costheta) );
146  G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
148 
149  //
150  // Calculate final state momenta in centre of mass system
151  //
152  if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
153  targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
154 
155  currentParticle.SetMomentum( -pf*sintheta*std::cos(phi)*GeV,
156  -pf*sintheta*std::sin(phi)*GeV,
157  -pf*costheta*GeV );
158  } else {
159 
160  currentParticle.SetMomentum( pf*sintheta*std::cos(phi)*GeV,
161  pf*sintheta*std::sin(phi)*GeV,
162  pf*costheta*GeV );
163  }
164  targetParticle.SetMomentum( -currentParticle.GetMomentum() );
165 
166  //
167  // Transform into lab system
168  //
169  currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
170  targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
171 
172  // Rotate final state particle vectors wrt incident momentum
173 
174  Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
175 
176  // Subtract binding energy
177 
178  G4double pp, pp1, ekin;
179  if( atomicWeight >= 1.5 )
180  {
181  const G4double cfa = 0.025*((atomicWeight-1.)/120.)*G4Exp(-(atomicWeight-1.)/120.);
182  pp1 = currentParticle.GetMomentum().mag()/MeV;
183  if( pp1 >= 1.0 )
184  {
185  ekin = currentParticle.GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV;
186  ekin = std::max( 0.0001*GeV, ekin );
187  currentParticle.SetKineticEnergy( ekin*MeV );
188  pp = currentParticle.GetTotalMomentum()/MeV;
189  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
190  }
191  pp1 = targetParticle.GetMomentum().mag()/MeV;
192  if( pp1 >= 1.0 )
193  {
194  ekin = targetParticle.GetKineticEnergy()/MeV - cfa*(1.0+normal()/2.)*GeV;
195  ekin = std::max( 0.0001*GeV, ekin );
196  targetParticle.SetKineticEnergy( ekin*MeV );
197  pp = targetParticle.GetTotalMomentum()/MeV;
198  targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
199  }
200  }
201  }
202 
203  // Get number of final state nucleons and nucleons remaining in
204  // target nucleus
205 
206  std::pair<G4int, G4int> finalStateNucleons =
207  GetFinalStateNucleons(originalTarget, vec, vecLen);
208  G4int protonsInFinalState = finalStateNucleons.first;
209  G4int neutronsInFinalState = finalStateNucleons.second;
210 
211  G4int PinNucleus = std::max(0,
212  G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState);
213  G4int NinNucleus = std::max(0,
214  G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState);
215 
216  if( atomicWeight >= 1.5 )
217  {
218  // Add black track particles
219  // npnb: number of proton/neutron black track particles
220  // ndta: number of deuterons, tritons, and alphas produced
221  // epnb: kinetic energy available for proton/neutron black track
222  // particles
223  // edta: kinetic energy available for deuteron/triton/alpha particles
224 
225  G4double epnb, edta;
226  G4int npnb=0, ndta=0;
227 
228  epnb = targetNucleus.GetPNBlackTrackEnergy(); // was enp1 in fortran code
229  edta = targetNucleus.GetDTABlackTrackEnergy(); // was enp3 in fortran code
230  const G4double pnCutOff = 0.0001; // GeV
231  const G4double dtaCutOff = 0.0001; // GeV
232  // const G4double kineticMinimum = 0.0001;
233  // const G4double kineticFactor = -0.010;
234  // G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
235  if( epnb >= pnCutOff )
236  {
237  npnb = G4Poisson( epnb/0.02 );
238  if( npnb > atomicWeight )npnb = G4int(atomicWeight);
239  if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1;
240  npnb = std::min( npnb, 127-vecLen );
241  }
242  if( edta >= dtaCutOff )
243  {
244  ndta = G4int(2.0 * G4Log(atomicWeight));
245  ndta = std::min( ndta, 127-vecLen );
246  }
247 
248  if (npnb == 0 && ndta == 0) npnb = 1;
249 
250  AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
251  PinNucleus, NinNucleus, targetNucleus,
252  vec, vecLen);
253  }
254 
255  //
256  // calculate time delay for nuclear reactions
257  //
258  if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
259  currentParticle.SetTOF( 1.0-500.0*G4Exp(-ekOriginal/0.04)*G4Log(G4UniformRand()) );
260  else
261  currentParticle.SetTOF( 1.0 );
262 
263  return true;
264 }
265 
266 /* end of file */