source: trunk/source/processes/hadronic/models/rpg/src/G4RPGTwoBody.cc@ 1347

Last change on this file since 1347 was 1340, checked in by garnier, 15 years ago

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