| 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 // 16.05.17 V. Grichine first implementation
27
28
29 #include "G4NeutronElectronElXsc.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4DynamicParticle.hh"
33 #include "G4ParticleDefinition.hh"
34 #include "G4ParticleTable.hh"
35 #include "G4IonTable.hh"
36 #include "G4HadTmpUtil.hh"
37 #include "G4NistManager.hh"
38 // #include "G4Integrator.hh"
39
40 #include "G4PhysicsLogVector.hh"
41 #include "G4PhysicsTable.hh"
42
43 #include "G4Neutron.hh"
44 #include "G4Electron.hh"
45
46
47 using namespace std;
48 using namespace CLHEP;
49
50 G4NeutronElectronElXsc::G4NeutronElectronElXsc()
51 : G4VCrossSectionDataSet("NuElectronCcXsc")
52 {
53 // neutron magneton squared
54
55 fM = neutron_mass_c2; // neutron mass
56 fM2 = fM*fM;
57 fme = electron_mass_c2;
58 fme2 = fme*fme;
59 fMv2 = 0.7056*GeV*GeV;
60 fee = fme;
61 fee2 = fee*fee;
62 fAm = 0.001;
63
64 fCofXsc = pi*fine_structure_const*fine_structure_const*hbarc*hbarc;
65 fCofXsc *= 3.6481; // neutron Fm^2(0)
66 fCofXsc /= fM*fM;
67
68 // G4cout<<"fCofXsc = "<<fCofXsc*GeV/cm2<<" cm2/GeV"<<G4endl;
69
70 // G4cout<<"hbarc = "<<hbarc/MeV/fermi<<" MeV*fermi"<<G4endl;
71
72 fCutEnergy = 0.; // default value
73
74 fEnergyBin = 200;
75 fMinEnergy = 1.*MeV;
76 fMaxEnergy = 10000.*GeV;
77
78 fEnergyXscVector = new G4PhysicsLogVector(fMinEnergy, fMaxEnergy, fEnergyBin);
79
80 for( G4int iTkin = 0; iTkin < fEnergyBin; iTkin++) fEnergyXscVector->PutValue(iTkin, fXscArray[iTkin]*microbarn);
81
82 fBiasingFactor = 1.;
83
84 // Initialise();
85 }
86
87 G4NeutronElectronElXsc::~G4NeutronElectronElXsc()
88 {
89 if( fEnergyXscVector )
90 {
91 delete fEnergyXscVector;
92 fEnergyXscVector = 0;
93 }
94 }
95
96 //////////////////////////////////////////////////////
97 //
98 // For neutrons in the precalculated energy interval
99
100 G4bool
101 G4NeutronElectronElXsc::IsElementApplicable( const G4DynamicParticle* aPart, G4int, const G4Material*)
102 {
103 G4bool result = false;
104 G4String pName = aPart->GetDefinition()->GetParticleName();
105 G4double Tkin = aPart->GetKineticEnergy();
106
107 if( pName == "neutron" &&
108 Tkin >= fMinEnergy &&
109 Tkin <= fMaxEnergy ) result = true;
110
111 return result;
112 }
113
114 //////////////////////////////////////////////////
115
116 void G4NeutronElectronElXsc::Initialise()
117 {
118 G4int iTkin;
119 G4double Tkin, rosxsc, xsc, delta, err=1.e-5;
120 const G4ThreeVector mDir = G4ThreeVector(0.,0.,1.);
121 const G4ParticleDefinition* pD = G4Neutron::Neutron();
122 G4Material* mat = G4NistManager::Instance()->FindOrBuildMaterial("G4_H");
123
124 G4DynamicParticle dP;
125
126 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
127 {
128 Tkin = fEnergyXscVector->GetLowEdgeEnergy(iTkin);
129 dP = G4DynamicParticle(pD, mDir, Tkin);
130 rosxsc = GetRosenbluthXsc(&dP, 1, mat);
131 fEnergyXscVector->PutValue(iTkin, rosxsc); // xscV.PutValue(evt, rosxsc); //
132 xsc= fEnergyXscVector->Value(Tkin); // xsc= xscV.GetLowEdgeEnergy(evt); //
133 delta = 0.5*std::abs( (rosxsc-xsc) )/(rosxsc+xsc);
134 if(delta > err) G4cout<<Tkin/GeV<<" GeV, rosxsc = "<<rosxsc/microbarn<<"umb, v-xsc = "<<xsc<<" umb"<<G4endl;
135 }
136 return;
137 }
138
139 ////////////////////////////////////////////////////
140
141 G4double G4NeutronElectronElXsc::
142 GetElementCrossSection(const G4DynamicParticle* aPart, G4int ZZ,
143 const G4Material*)
144 {
145 G4double result = 0., Tkin;
146
147 Tkin = aPart->GetKineticEnergy();
148
149 result = fEnergyXscVector->Value(Tkin);
150
151 result *= ZZ; // incoherent sum over all element electrons
152
153 result *= fBiasingFactor;
154
155 return result;
156 }
157
158 ////////////////////////////////////////////////////
159 //
160 // Integration of the Rosenbluth differential xsc
161
162 G4double G4NeutronElectronElXsc::
163 GetRosenbluthXsc(const G4DynamicParticle* aPart, G4int ZZ,
164 const G4Material*)
165 {
166 G4double result = 0., momentum;
167
168 fee = aPart->GetTotalEnergy()*fme/fM;
169 fee2 = fee*fee;
170 momentum = sqrt( fee2 - fme2 );
171 fAm = CalculateAm(momentum);
172
173 // G4Integrator<G4NeutronElectronElXsc, G4double(G4NeutronElectronElXsc::*)(G4double)> integral;
174
175 // result = integral.Legendre96( this, &G4NeutronElectronElXsc::XscIntegrand, 0., 1. );
176
177 result *= fCofXsc;
178
179 result *= ZZ; // incoherent sum over all element electrons
180
181 return result;
182 }
183
184 /////////////////////////////////////////
185 //
186 // Rosenbluth relation in the neutron rest frame,
187 // x = sin^2(theta/2), theta is the electron scattering angle
188 // Magnetic form factor in the dipole approximation.
189
190 G4double G4NeutronElectronElXsc::XscIntegrand(G4double x)
191 {
192 G4double result = 1., q2, znq2, znf, znf2, znf4;
193
194 znq2 = 1. + 2.*fee*x/fM;
195
196 q2 = 4.*fee2*x/znq2;
197
198 znf = 1 + q2/fMv2;
199 znf2 = znf*znf;
200 znf4 = znf2*znf2;
201
202 result /= ( x + fAm )*znq2*znq2*znf4;
203
204 result *= ( 1 - x )/( 1 + q2/4./fM2 ) + 2.*x;
205
206 return result;
207 }
208
209 //////////////////////////////////////////////////////////
210 //
211 // Rosenbluth xsc in microbarn from 1*MeV to 10*Tev, 200 points
212
213 const G4double G4NeutronElectronElXsc::fXscArray[200] = {
214 1.52681, 1.54903, 1.57123, 1.59341, 1.61556, 1.63769, 1.6598, 1.68189, 1.70396,
215 1.72601, 1.74805, 1.77007, 1.79208, 1.81407, 1.83605, 1.85801, 1.87997, 1.90192,
216 1.92385, 1.94578, 1.96771, 1.98962, 2.01154, 2.03345, 2.05535, 2.07725, 2.09915,
217 2.12105, 2.14295, 2.16485, 2.18675, 2.20865, 2.23055, 2.25244, 2.27433, 2.29621,
218 2.31807, 2.33992, 2.36173, 2.38351, 2.40524, 2.42691, 2.4485, 2.47, 2.49138,
219 2.51262, 2.53369, 2.55457, 2.57524, 2.59565, 2.61577, 2.63559, 2.65505, 2.67414,
220 2.69281, 2.71104, 2.72881, 2.74607, 2.76282, 2.77903, 2.79467, 2.80974, 2.82422,
221 2.83811, 2.85139, 2.86408, 2.87616, 2.88764, 2.89854, 2.90885, 2.91859, 2.92777,
222 2.93641, 2.94453, 2.95213, 2.95924, 2.96588, 2.97207, 2.97782, 2.98316, 2.98811,
223 2.99268, 2.9969, 3.00078, 3.00435, 3.00761, 3.01059, 3.01331, 3.01578, 3.01801,
224 3.02003, 3.02185, 3.02347, 3.02491, 3.02619, 3.02732, 3.0283, 3.02915, 3.02988,
225 3.03049, 3.03099, 3.03139, 3.03169, 3.03191, 3.03203, 3.03208, 3.03205, 3.03195,
226 3.03177, 3.03152, 3.0312, 3.03081, 3.03034, 3.0298, 3.02919, 3.02849, 3.02771,
227 3.02684, 3.02588, 3.02482, 3.02365, 3.02237, 3.02097, 3.01943, 3.01775, 3.0159,
228 3.01389, 3.01169, 3.00929, 3.00666, 3.00379, 3.00065, 2.99722, 2.99347, 2.98936,
229 2.98487, 2.97996, 2.97459, 2.9687, 2.96226, 2.9552, 2.94748, 2.93903, 2.92977,
230 2.91965, 2.90858, 2.89649, 2.88329, 2.86889, 2.85321, 2.83615, 2.81764, 2.7976,
231 2.77594, 2.7526, 2.72754, 2.70071, 2.67209, 2.64171, 2.60957, 2.57575, 2.54031,
232 2.50336, 2.46504, 2.42548, 2.38484, 2.34328, 2.30099, 2.2581, 2.21478, 2.17115,
233 2.12735, 2.08345, 2.03954, 1.99569, 1.95191, 1.90825, 1.86471, 1.82129, 1.77799,
234 1.7348, 1.69171, 1.64869, 1.60575, 1.56286, 1.52, 1.47718, 1.43437, 1.39157,
235 1.34877, 1.30596, 1.26314, 1.22031, 1.17746, 1.13459, 1.0917, 1.04879, 1.00585,
236 0.962892, 0.919908 };
237