Geant4.10
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4RPGTwoBody Class Reference

#include <G4RPGTwoBody.hh>

Inheritance diagram for G4RPGTwoBody:
G4RPGReaction

Public Member Functions

 G4RPGTwoBody ()
 
G4bool ReactionStage (const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
 
- Public Member Functions inherited from G4RPGReaction
 G4RPGReaction ()
 
virtual ~G4RPGReaction ()
 
G4bool ReactionStage (const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
 
void AddBlackTrackParticles (const G4double, const G4int, const G4double, const G4int, const G4ReactionProduct &, G4int, G4int, const G4Nucleus &, G4FastVector< G4ReactionProduct, 256 > &, G4int &)
 
G4double GenerateNBodyEvent (const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
G4double GenerateNBodyEventT (const G4double totalEnergy, const G4bool constantCrossSection, std::vector< G4ReactionProduct * > &list)
 
void NuclearReaction (G4FastVector< G4ReactionProduct, 4 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4Nucleus &aNucleus, const G4double theAtomicMass, const G4double *massVec)
 

Additional Inherited Members

- Protected Member Functions inherited from G4RPGReaction
void Rotate (const G4double numberofFinalStateNucleons, const G4ThreeVector &temp, const G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, const G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
void Defs1 (const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
std::pair< G4int, G4int > GetFinalStateNucleons (const G4DynamicParticle *originalTarget, const G4FastVector< G4ReactionProduct, 256 > &vec, const G4int &vecLen)
 
void MomentumCheck (const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
G4double normal ()
 
G4ThreeVector Isotropic (const G4double &)
 

Detailed Description

Definition at line 50 of file G4RPGTwoBody.hh.

Constructor & Destructor Documentation

G4RPGTwoBody::G4RPGTwoBody ( )

Definition at line 39 of file G4RPGTwoBody.cc.

40  : G4RPGReaction() {}

Member Function Documentation

G4bool G4RPGTwoBody::ReactionStage ( const G4HadProjectile *  ,
G4ReactionProduct &  modifiedOriginal,
G4bool &  ,
const G4DynamicParticle *  originalTarget,
G4ReactionProduct &  targetParticle,
G4bool &  ,
const G4Nucleus &  targetNucleus,
G4ReactionProduct &  currentParticle,
G4FastVector< G4ReactionProduct, 256 > &  vec,
G4int &  vecLen,
G4bool  ,
G4ReactionProduct &   
)

Definition at line 44 of file G4RPGTwoBody.cc.

References G4RPGReaction::AddBlackTrackParticles(), test::b, G4RPGReaction::Defs1(), G4Poisson(), G4UniformRand, G4Nucleus::GetA_asInt(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4RPGReaction::GetFinalStateNucleons(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), python.hepunit::GeV, G4ReactionProduct::Lorentz(), CLHEP::Hep3Vector::mag(), G4INCL::Math::max(), python.hepunit::MeV, G4INCL::Math::min(), G4RPGReaction::normal(), G4InuclParticleNames::pp, G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTOF(), G4ReactionProduct::SetTotalEnergy(), and python.hepunit::twopi.

Referenced by G4RPGInelastic::CalculateMomenta().

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*std::log(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 = std::exp(-btrang) - 1.0;
144  G4double costheta = 1.0 + 2*std::log( 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));
147  G4double phi = twopi * G4UniformRand();
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.)*std::exp(-(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 * std::log(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*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
260  else
261  currentParticle.SetTOF( 1.0 );
262 
263  return true;
264 }
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:199
G4RPGReaction::Defs1
void Defs1(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
Definition: G4RPGReaction.cc:693
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:132
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:165
G4Nucleus::GetDTABlackTrackEnergy
G4double GetDTABlackTrackEnergy() const
Definition: G4Nucleus.hh:153
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:181
G4int
int G4int
Definition: G4Types.hh:78
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:107
G4RPGReaction::GetFinalStateNucleons
std::pair< G4int, G4int > GetFinalStateNucleons(const G4DynamicParticle *originalTarget, const G4FastVector< G4ReactionProduct, 256 > &vec, const G4int &vecLen)
Definition: G4RPGReaction.cc:1011
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:147
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:138
G4RPGReaction::normal
G4double normal()
Definition: G4RPGReaction.cc:685
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
G4RPGReaction::AddBlackTrackParticles
void AddBlackTrackParticles(const G4double, const G4int, const G4double, const G4int, const G4ReactionProduct &, G4int, G4int, const G4Nucleus &, G4FastVector< G4ReactionProduct, 256 > &, G4int &)
Definition: G4RPGReaction.cc:57
G4ReactionProduct::SetTOF
void SetTOF(const G4double t)
Definition: G4ReactionProduct.hh:153
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4Nucleus::GetPNBlackTrackEnergy
G4double GetPNBlackTrackEnergy() const
Definition: G4Nucleus.hh:150
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:57:01 for Geant4.10 by   doxygen 1.8.7