Geant4  9.6.p02
 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$
27 //
28 
29 #include <iostream>
30 #include <signal.h>
31 
32 #include "G4RPGTwoBody.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*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 }
265 
266 /* end of file */