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