| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Geant4 Header : G4NeutronElectronElModel
28 //
29 // 16.5.17: V.Grichine
30 //
31
32 #include "G4NeutronElectronElModel.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4ParticleTable.hh"
35 #include "G4ParticleDefinition.hh"
36 #include "G4IonTable.hh"
37 #include "Randomize.hh"
38 #include "G4Integrator.hh"
39 #include "G4Electron.hh"
40 #include "G4PhysicsTable.hh"
41 #include "G4PhysicsLogVector.hh"
42 #include "G4PhysicsFreeVector.hh"
43 #include "G4PhysicsModelCatalog.hh"
44
45
46 using namespace std;
47 using namespace CLHEP;
48
49 G4NeutronElectronElModel::G4NeutronElectronElM
50 : G4HadronElastic(name)
51 {
52 secID = G4PhysicsModelCatalog::GetModelID( "
53
54 // neutron magneton squared
55
56 fM = neutron_mass_c2; // neutron mass
57 fM2 = fM*fM;
58 fme = electron_mass_c2;
59 fme2 = fme*fme;
60 fMv2 = 0.7056*GeV*GeV;
61
62 SetMinEnergy( 0.001*GeV );
63 SetMaxEnergy( 10.*TeV );
64 SetLowestEnergyLimit(1.e-6*eV);
65
66 theElectron = G4Electron::Electron();
67 // PDG2016: sin^2 theta Weinberg
68
69 fEnergyBin = 200;
70 fMinEnergy = 1.*MeV;
71 fMaxEnergy = 10000.*GeV;
72 fEnergyVector = new G4PhysicsLogVector(fMinE
73
74 fAngleBin = 500;
75 fAngleTable = 0;
76
77 fCutEnergy = 0.; // default value
78
79 Initialise();
80 }
81
82 //////////////////////////////////////////////
83
84 G4NeutronElectronElModel::~G4NeutronElectronEl
85 {
86 if( fEnergyVector )
87 {
88 delete fEnergyVector;
89 fEnergyVector = nullptr;
90 }
91 if( fAngleTable )
92 {
93 fAngleTable->clearAndDestroy();
94 delete fAngleTable;
95 fAngleTable = nullptr;
96 }
97 }
98
99 /////////////////////////////////////////
100
101 void G4NeutronElectronElModel::ModelDescriptio
102 {
103 outFile << "G4NeutronElectronElModel is a ne
104 << "model which uses the standard model \n
105 << "transfer parameterization. The model
106 }
107
108 //////////////////////////////////////////////
109
110 G4bool G4NeutronElectronElModel::IsApplicable(
111 {
112 G4String pName = aTrack.GetDefinition()->Get
113 G4double energy = aTrack.GetTotalEnergy();
114
115 return (pName == "neutron" && energy >= fMin
116 }
117
118 //////////////////////////////////////////////
119
120 void G4NeutronElectronElModel::Initialise()
121 {
122 G4double result = 0., sum, Tkin, dt, t1, t2;
123 G4int iTkin, jTransfer;
124 G4Integrator<G4NeutronElectronElModel, G4dou
125
126 fAngleTable = new G4PhysicsTable(fEnergyBin)
127
128 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
129 {
130 Tkin = fEnergyVector->GetLowEdgeEnergy(iT
131 fAm = CalculateAm(Tkin);
132 dt = 1./fAngleBin;
133
134 G4PhysicsFreeVector* vectorT = new G4Physi
135
136 sum = 0.;
137
138 for( jTransfer = 0; jTransfer < fAngleBin;
139 {
140 t1 = dt*jTransfer;
141 t2 = t1 + dt;
142
143 result = integral.Legendre96( this, &G4N
144
145 sum += result;
146 // G4cout<<sum<<", ";
147 vectorT->PutValue(jTransfer, t1, sum);
148 }
149 // G4cout<<G4endl;
150 fAngleTable->insertAt(iTkin,vectorT);
151 }
152 return;
153 }
154
155 //////////////////////////////////////////////
156 //
157 // sample recoil electron energy in lab frame
158
159 G4double G4NeutronElectronElModel::SampleSin2H
160 {
161 G4double result = 0., position;
162 G4int iTkin, iTransfer;
163
164 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
165 {
166 if( Tkin < fEnergyVector->Energy(iTkin) )
167 }
168 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBi
169 if ( iTkin < 0 ) iTkin = 0; // aga
170
171 position = (*(*fAngleTable)(iTkin))(fAngle
172
173 // G4cout<<"position = "<<position<<G4endl
174
175 for( iTransfer = 0; iTransfer < fAngleBin;
176 {
177 if( position <= (*(*fAngleTable)(iTkin))
178 }
179 if (iTransfer >= fAngleBin-1) iTransfer =
180
181 // G4cout<<"iTransfer = "<<iTransfer<<G4en
182
183 result = GetTransfer(iTkin, iTransfer, pos
184
185 // G4cout<<"t = "<<t<<G4endl;
186
187
188 return result;
189 }
190
191 //////////////////////////////////////////////
192
193 G4double
194 G4NeutronElectronElModel:: GetTransfer( G4int
195 {
196 G4double x1, x2, y1, y2, randTransfer, delta
197
198 if( iTransfer == 0 || iTransfer == fAngleBi
199 {
200 randTransfer = (*fAngleTable)(iTkin)->Ener
201 }
202 else
203 {
204 if ( iTransfer >= G4int((*fAngleTable)(iTk
205 {
206 iTransfer = G4int((*fAngleTable)(iTkin)-
207 }
208 y1 = (*(*fAngleTable)(iTkin))(iTransfer-1)
209 y2 = (*(*fAngleTable)(iTkin))(iTransfer);
210
211 x1 = (*fAngleTable)(iTkin)->Energy(iTransf
212 x2 = (*fAngleTable)(iTkin)->Energy(iTransf
213
214 delta = y2 - y1;
215 mean = y2 + y1;
216
217 if ( x1 == x2 ) randTransfer = x2;
218 else
219 {
220 if ( delta < epsilon*mean )
221 {
222 randTransfer = x1 + ( x2 - x1 )*G4Unif
223 }
224 else
225 {
226 randTransfer = x1 + ( position - y1 )*
227 }
228 }
229 }
230 return randTransfer;
231 }
232
233 //////////////////////////////////////////////
234 //
235 // Rosenbluth relation (ultra-relativistic!) i
236 // x = sin^2(theta/2), theta is the electron s
237 // Magnetic form factor in the dipole approxim
238
239 G4double G4NeutronElectronElModel::XscIntegran
240 {
241 G4double result = 1., q2, znq2, znf, znf2, z
242
243 znq2 = 1. + 2.*fee*x/fM;
244
245 q2 = 4.*fee2*x/znq2;
246
247 znf = 1 + q2/fMv2;
248 znf2 = znf*znf;
249 znf4 = znf2*znf2;
250
251 result /= ( x + fAm )*znq2*znq2*znf4;
252
253 result *= ( 1 - x )/( 1 + q2/4./fM2 ) + 2.*x
254
255 return result;
256 }
257
258 //////////////////////////////////////////////
259 //
260 //
261
262 G4HadFinalState* G4NeutronElectronElModel::App
263 const G4HadProjectile& aTrack, G4Nucleus&
264 {
265 theParticleChange.Clear();
266
267 const G4HadProjectile* aParticle = &aTrack;
268 G4double Tkin = aParticle->GetKineticEnergy(
269 fAm = CalculateAm( Tkin);
270 // G4double En = aParticle->GetTotalEnergy
271
272 if( Tkin <= LowestEnergyLimit() )
273 {
274 theParticleChange.SetEnergyChange(Tkin);
275 theParticleChange.SetMomentumChange(aTrack
276 return &theParticleChange;
277 }
278 // sample e-scattering angle and make final
279
280 G4double sin2ht = SampleSin2HalfTheta( Tkin)
281
282 // G4cout<<"sin2ht = "<<sin2ht<<G4endl;
283
284 G4double eTkin = fee; // fM;
285
286 eTkin /= 1.+2.*fee*sin2ht/fM; // fme/En + 2*
287
288 eTkin -= fme;
289
290 // G4cout<<"eTkin = "<<eTkin<<G4endl;
291
292 if( eTkin > fCutEnergy )
293 {
294 G4double ePlab = sqrt( eTkin*(eTkin + 2.*f
295
296 // G4cout<<"ePlab = "<<ePlab<<G4endl;
297
298 G4double cost = 1. - 2*sin2ht;
299
300 if( cost > 1. ) cost = 1.;
301 if( cost < -1. ) cost = -1.;
302
303 G4double sint = std::sqrt( (1.0 - cost)*(1
304 G4double phi = G4UniformRand()*CLHEP::two
305
306 G4ThreeVector eP( sint*std::cos(phi), sint
307 eP *= ePlab;
308 G4LorentzVector lvt2( eP, eTkin + electron
309
310 G4LorentzVector lvp1 = aParticle->Get4Mome
311 G4LorentzVector lvt1(0.,0.,0.,electron_mas
312 G4LorentzVector lvsum = lvp1+lvt1;
313
314 G4ThreeVector bst = lvp1.boostVector();
315 lvt2.boost(bst);
316
317 // G4cout<<"lvt2 = "<<lvt2<<G4endl;
318
319 G4DynamicParticle * aSec = new G4DynamicPa
320 theParticleChange.AddSecondary( aSec, secI
321
322 G4LorentzVector lvp2 = lvsum-lvt2;
323
324 // G4cout<<"lvp2 = "<<lvp2<<G4endl;
325
326 G4double Tkin2 = lvp2.e()-aParticle->GetDe
327 theParticleChange.SetEnergyChange(Tkin2);
328 theParticleChange.SetMomentumChange(lvp2.v
329 }
330 else if( eTkin > 0.0 )
331 {
332 theParticleChange.SetLocalEnergyDeposit( e
333 Tkin -= eTkin;
334
335 if( Tkin > 0. )
336 {
337 theParticleChange.SetEnergyChange( Tkin
338 theParticleChange.SetMomentumChange( aTr
339 }
340 }
341 else
342 {
343 theParticleChange.SetEnergyChange( Tkin );
344 theParticleChange.SetMomentumChange( aTrac
345 }
346 return &theParticleChange;
347 }
348
349 //
350 //
351 ///////////////////////////
352