| Geant4 Cross Reference |
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 //
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::G4NeutronElectronElModel(const G4String& name)
50 : G4HadronElastic(name)
51 {
52 secID = G4PhysicsModelCatalog::GetModelID( "model_" + name );
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(fMinEnergy, fMaxEnergy, fEnergyBin, false);
73
74 fAngleBin = 500;
75 fAngleTable = 0;
76
77 fCutEnergy = 0.; // default value
78
79 Initialise();
80 }
81
82 ////////////////////////////////////////////////
83
84 G4NeutronElectronElModel::~G4NeutronElectronElModel()
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::ModelDescription(std::ostream& outFile) const
102 {
103 outFile << "G4NeutronElectronElModel is a neutrino-electron (neutral current) elastic scattering\n"
104 << "model which uses the standard model \n"
105 << "transfer parameterization. The model is fully relativistic\n";
106 }
107
108 /////////////////////////////////////////////////////////
109
110 G4bool G4NeutronElectronElModel::IsApplicable(const G4HadProjectile & aTrack, G4Nucleus&)
111 {
112 G4String pName = aTrack.GetDefinition()->GetParticleName();
113 G4double energy = aTrack.GetTotalEnergy();
114
115 return (pName == "neutron" && energy >= fMinEnergy && energy <= fMaxEnergy);
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, G4double(G4NeutronElectronElModel::*)(G4double)> integral;
125
126 fAngleTable = new G4PhysicsTable(fEnergyBin);
127
128 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
129 {
130 Tkin = fEnergyVector->GetLowEdgeEnergy(iTkin);
131 fAm = CalculateAm(Tkin);
132 dt = 1./fAngleBin;
133
134 G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fAngleBin);
135
136 sum = 0.;
137
138 for( jTransfer = 0; jTransfer < fAngleBin; jTransfer++)
139 {
140 t1 = dt*jTransfer;
141 t2 = t1 + dt;
142
143 result = integral.Legendre96( this, &G4NeutronElectronElModel::XscIntegrand, t1, t2 );
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::SampleSin2HalfTheta(G4double Tkin)
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) ) break;
167 }
168 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
169 if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
170
171 position = (*(*fAngleTable)(iTkin))(fAngleBin-1)*G4UniformRand();
172
173 // G4cout<<"position = "<<position<<G4endl;
174
175 for( iTransfer = 0; iTransfer < fAngleBin; iTransfer++)
176 {
177 if( position <= (*(*fAngleTable)(iTkin))(iTransfer) ) break;
178 }
179 if (iTransfer >= fAngleBin-1) iTransfer = fAngleBin-1;
180
181 // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
182
183 result = GetTransfer(iTkin, iTransfer, position);
184
185 // G4cout<<"t = "<<t<<G4endl;
186
187
188 return result;
189 }
190
191 /////////////////////////////////////////////////
192
193 G4double
194 G4NeutronElectronElModel:: GetTransfer( G4int iTkin, G4int iTransfer, G4double position )
195 {
196 G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6;
197
198 if( iTransfer == 0 || iTransfer == fAngleBin-1 )
199 {
200 randTransfer = (*fAngleTable)(iTkin)->Energy(iTransfer);
201 }
202 else
203 {
204 if ( iTransfer >= G4int((*fAngleTable)(iTkin)->GetVectorLength()) )
205 {
206 iTransfer = G4int((*fAngleTable)(iTkin)->GetVectorLength() - 1);
207 }
208 y1 = (*(*fAngleTable)(iTkin))(iTransfer-1);
209 y2 = (*(*fAngleTable)(iTkin))(iTransfer);
210
211 x1 = (*fAngleTable)(iTkin)->Energy(iTransfer-1);
212 x2 = (*fAngleTable)(iTkin)->Energy(iTransfer);
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 )*G4UniformRand();
223 }
224 else
225 {
226 randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 );
227 }
228 }
229 }
230 return randTransfer;
231 }
232
233 //////////////////////////////////////////////////////////////
234 //
235 // Rosenbluth relation (ultra-relativistic!) in the neutron rest frame,
236 // x = sin^2(theta/2), theta is the electron scattering angle
237 // Magnetic form factor in the dipole approximation.
238
239 G4double G4NeutronElectronElModel::XscIntegrand(G4double x)
240 {
241 G4double result = 1., q2, znq2, znf, znf2, znf4;
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::ApplyYourself(
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.Get4Momentum().vect().unit());
276 return &theParticleChange;
277 }
278 // sample e-scattering angle and make final state in lab frame
279
280 G4double sin2ht = SampleSin2HalfTheta( Tkin); // in n-rrest frame
281
282 // G4cout<<"sin2ht = "<<sin2ht<<G4endl;
283
284 G4double eTkin = fee; // fM;
285
286 eTkin /= 1.+2.*fee*sin2ht/fM; // fme/En + 2*sin2ht;
287
288 eTkin -= fme;
289
290 // G4cout<<"eTkin = "<<eTkin<<G4endl;
291
292 if( eTkin > fCutEnergy )
293 {
294 G4double ePlab = sqrt( eTkin*(eTkin + 2.*fme) );
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.0 + cost) );
304 G4double phi = G4UniformRand()*CLHEP::twopi;
305
306 G4ThreeVector eP( sint*std::cos(phi), sint*std::sin(phi), cost );
307 eP *= ePlab;
308 G4LorentzVector lvt2( eP, eTkin + electron_mass_c2 ); // recoil e- in n-rest frame
309
310 G4LorentzVector lvp1 = aParticle->Get4Momentum();
311 G4LorentzVector lvt1(0.,0.,0.,electron_mass_c2);
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 G4DynamicParticle( theElectron, lvt2 );
320 theParticleChange.AddSecondary( aSec, secID );
321
322 G4LorentzVector lvp2 = lvsum-lvt2;
323
324 // G4cout<<"lvp2 = "<<lvp2<<G4endl;
325
326 G4double Tkin2 = lvp2.e()-aParticle->GetDefinition()->GetPDGMass();
327 theParticleChange.SetEnergyChange(Tkin2);
328 theParticleChange.SetMomentumChange(lvp2.vect().unit());
329 }
330 else if( eTkin > 0.0 )
331 {
332 theParticleChange.SetLocalEnergyDeposit( eTkin );
333 Tkin -= eTkin;
334
335 if( Tkin > 0. )
336 {
337 theParticleChange.SetEnergyChange( Tkin );
338 theParticleChange.SetMomentumChange( aTrack.Get4Momentum().vect().unit() );
339 }
340 }
341 else
342 {
343 theParticleChange.SetEnergyChange( Tkin );
344 theParticleChange.SetMomentumChange( aTrack.Get4Momentum().vect().unit() );
345 }
346 return &theParticleChange;
347 }
348
349 //
350 //
351 ///////////////////////////
352