| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 16.05.17 V. Grichine first implementation 26 // 16.05.17 V. Grichine first implementation
27 27
28 28
29 #include "G4NeutronElectronElXsc.hh" 29 #include "G4NeutronElectronElXsc.hh"
30 #include "G4PhysicalConstants.hh" 30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
32 #include "G4DynamicParticle.hh" 32 #include "G4DynamicParticle.hh"
33 #include "G4ParticleDefinition.hh" 33 #include "G4ParticleDefinition.hh"
34 #include "G4ParticleTable.hh" 34 #include "G4ParticleTable.hh"
35 #include "G4IonTable.hh" 35 #include "G4IonTable.hh"
36 #include "G4HadTmpUtil.hh" 36 #include "G4HadTmpUtil.hh"
37 #include "G4NistManager.hh" 37 #include "G4NistManager.hh"
38 // #include "G4Integrator.hh" 38 // #include "G4Integrator.hh"
39 39
40 #include "G4PhysicsLogVector.hh" 40 #include "G4PhysicsLogVector.hh"
41 #include "G4PhysicsTable.hh" 41 #include "G4PhysicsTable.hh"
42 42
43 #include "G4Neutron.hh" 43 #include "G4Neutron.hh"
44 #include "G4Electron.hh" 44 #include "G4Electron.hh"
45 45
46 46
47 using namespace std; 47 using namespace std;
48 using namespace CLHEP; 48 using namespace CLHEP;
49 49
50 G4NeutronElectronElXsc::G4NeutronElectronElXsc 50 G4NeutronElectronElXsc::G4NeutronElectronElXsc()
51 : G4VCrossSectionDataSet("NuElectronCcXsc") 51 : G4VCrossSectionDataSet("NuElectronCcXsc")
52 { 52 {
53 // neutron magneton squared 53 // neutron magneton squared
54 54
55 fM = neutron_mass_c2; // neutron mass 55 fM = neutron_mass_c2; // neutron mass
56 fM2 = fM*fM; 56 fM2 = fM*fM;
57 fme = electron_mass_c2; 57 fme = electron_mass_c2;
58 fme2 = fme*fme; 58 fme2 = fme*fme;
59 fMv2 = 0.7056*GeV*GeV; 59 fMv2 = 0.7056*GeV*GeV;
60 fee = fme; 60 fee = fme;
61 fee2 = fee*fee; 61 fee2 = fee*fee;
62 fAm = 0.001; 62 fAm = 0.001;
63 63
64 fCofXsc = pi*fine_structure_const*fine_stru 64 fCofXsc = pi*fine_structure_const*fine_structure_const*hbarc*hbarc;
65 fCofXsc *= 3.6481; // neutron Fm^2(0) 65 fCofXsc *= 3.6481; // neutron Fm^2(0)
66 fCofXsc /= fM*fM; 66 fCofXsc /= fM*fM;
67 67
68 // G4cout<<"fCofXsc = "<<fCofXsc*GeV/cm2<<" 68 // G4cout<<"fCofXsc = "<<fCofXsc*GeV/cm2<<" cm2/GeV"<<G4endl;
69 69
70 // G4cout<<"hbarc = "<<hbarc/MeV/fermi<<" Me 70 // G4cout<<"hbarc = "<<hbarc/MeV/fermi<<" MeV*fermi"<<G4endl;
71 71
72 fCutEnergy = 0.; // default value 72 fCutEnergy = 0.; // default value
73 73
74 fEnergyBin = 200; 74 fEnergyBin = 200;
75 fMinEnergy = 1.*MeV; 75 fMinEnergy = 1.*MeV;
76 fMaxEnergy = 10000.*GeV; 76 fMaxEnergy = 10000.*GeV;
77 77
78 fEnergyXscVector = new G4PhysicsLogVector(fM 78 fEnergyXscVector = new G4PhysicsLogVector(fMinEnergy, fMaxEnergy, fEnergyBin);
79 79
80 for( G4int iTkin = 0; iTkin < fEnergyBin; iT 80 for( G4int iTkin = 0; iTkin < fEnergyBin; iTkin++) fEnergyXscVector->PutValue(iTkin, fXscArray[iTkin]*microbarn);
81 81
82 fBiasingFactor = 1.; 82 fBiasingFactor = 1.;
83 83
84 // Initialise(); 84 // Initialise();
85 } 85 }
86 86
87 G4NeutronElectronElXsc::~G4NeutronElectronElXs 87 G4NeutronElectronElXsc::~G4NeutronElectronElXsc()
88 { 88 {
89 if( fEnergyXscVector ) 89 if( fEnergyXscVector )
90 { 90 {
91 delete fEnergyXscVector; 91 delete fEnergyXscVector;
92 fEnergyXscVector = 0; 92 fEnergyXscVector = 0;
93 } 93 }
94 } 94 }
95 95
96 ////////////////////////////////////////////// 96 //////////////////////////////////////////////////////
97 // 97 //
98 // For neutrons in the precalculated energy in 98 // For neutrons in the precalculated energy interval
99 99
100 G4bool 100 G4bool
101 G4NeutronElectronElXsc::IsElementApplicable( c 101 G4NeutronElectronElXsc::IsElementApplicable( const G4DynamicParticle* aPart, G4int, const G4Material*)
102 { 102 {
103 G4bool result = false; 103 G4bool result = false;
104 G4String pName = aPart->GetDefinition()->Get 104 G4String pName = aPart->GetDefinition()->GetParticleName();
105 G4double Tkin = aPart->GetKineticEnergy(); 105 G4double Tkin = aPart->GetKineticEnergy();
106 106
107 if( pName == "neutron" && 107 if( pName == "neutron" &&
108 Tkin >= fMinEnergy && 108 Tkin >= fMinEnergy &&
109 Tkin <= fMaxEnergy ) result = true; 109 Tkin <= fMaxEnergy ) result = true;
110 110
111 return result; 111 return result;
112 } 112 }
113 113
114 ////////////////////////////////////////////// 114 //////////////////////////////////////////////////
115 115
116 void G4NeutronElectronElXsc::Initialise() 116 void G4NeutronElectronElXsc::Initialise()
117 { 117 {
118 G4int iTkin; 118 G4int iTkin;
119 G4double Tkin, rosxsc, xsc, delta, err=1.e-5 119 G4double Tkin, rosxsc, xsc, delta, err=1.e-5;
120 const G4ThreeVector mDir = G4ThreeVector(0., 120 const G4ThreeVector mDir = G4ThreeVector(0.,0.,1.);
121 const G4ParticleDefinition* pD = G4Neutron:: 121 const G4ParticleDefinition* pD = G4Neutron::Neutron();
122 G4Material* mat = G4NistManager::Instance()- 122 G4Material* mat = G4NistManager::Instance()->FindOrBuildMaterial("G4_H");
123 123
124 G4DynamicParticle dP; 124 G4DynamicParticle dP;
125 125
126 for( iTkin = 0; iTkin < fEnergyBin; iTkin++) 126 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
127 { 127 {
128 Tkin = fEnergyXscVector->GetLowEdgeEnergy 128 Tkin = fEnergyXscVector->GetLowEdgeEnergy(iTkin);
129 dP = G4DynamicParticle(pD, mDir, Tkin); 129 dP = G4DynamicParticle(pD, mDir, Tkin);
130 rosxsc = GetRosenbluthXsc(&dP, 1, mat); 130 rosxsc = GetRosenbluthXsc(&dP, 1, mat);
131 fEnergyXscVector->PutValue(iTkin, rosxsc); 131 fEnergyXscVector->PutValue(iTkin, rosxsc); // xscV.PutValue(evt, rosxsc); //
132 xsc= fEnergyXscVector->Value(Tkin); // xsc 132 xsc= fEnergyXscVector->Value(Tkin); // xsc= xscV.GetLowEdgeEnergy(evt); //
133 delta = 0.5*std::abs( (rosxsc-xsc) )/(rosx 133 delta = 0.5*std::abs( (rosxsc-xsc) )/(rosxsc+xsc);
134 if(delta > err) G4cout<<Tkin/GeV<<" GeV, r 134 if(delta > err) G4cout<<Tkin/GeV<<" GeV, rosxsc = "<<rosxsc/microbarn<<"umb, v-xsc = "<<xsc<<" umb"<<G4endl;
135 } 135 }
136 return; 136 return;
137 } 137 }
138 138
139 ////////////////////////////////////////////// 139 ////////////////////////////////////////////////////
140 140
141 G4double G4NeutronElectronElXsc:: 141 G4double G4NeutronElectronElXsc::
142 GetElementCrossSection(const G4DynamicParticle 142 GetElementCrossSection(const G4DynamicParticle* aPart, G4int ZZ,
143 const G4Material*) 143 const G4Material*)
144 { 144 {
145 G4double result = 0., Tkin; 145 G4double result = 0., Tkin;
146 146
147 Tkin = aPart->GetKineticEnergy(); 147 Tkin = aPart->GetKineticEnergy();
148 148
149 result = fEnergyXscVector->Value(Tkin); 149 result = fEnergyXscVector->Value(Tkin);
150 150
151 result *= ZZ; // incoherent sum over all e 151 result *= ZZ; // incoherent sum over all element electrons
152 152
153 result *= fBiasingFactor; 153 result *= fBiasingFactor;
154 154
155 return result; 155 return result;
156 } 156 }
157 157
158 ////////////////////////////////////////////// 158 ////////////////////////////////////////////////////
159 // 159 //
160 // Integration of the Rosenbluth differential 160 // Integration of the Rosenbluth differential xsc
161 161
162 G4double G4NeutronElectronElXsc:: 162 G4double G4NeutronElectronElXsc::
163 GetRosenbluthXsc(const G4DynamicParticle* aPar 163 GetRosenbluthXsc(const G4DynamicParticle* aPart, G4int ZZ,
164 const G4Material*) 164 const G4Material*)
165 { 165 {
166 G4double result = 0., momentum; 166 G4double result = 0., momentum;
167 167
168 fee = aPart->GetTotalEnergy()*fme/fM; 168 fee = aPart->GetTotalEnergy()*fme/fM;
169 fee2 = fee*fee; 169 fee2 = fee*fee;
170 momentum = sqrt( fee2 - fme2 ); 170 momentum = sqrt( fee2 - fme2 );
171 fAm = CalculateAm(momentum); 171 fAm = CalculateAm(momentum);
172 172
173 // G4Integrator<G4NeutronElectronElXsc, G4do 173 // G4Integrator<G4NeutronElectronElXsc, G4double(G4NeutronElectronElXsc::*)(G4double)> integral;
174 174
175 // result = integral.Legendre96( this, &G4Ne 175 // result = integral.Legendre96( this, &G4NeutronElectronElXsc::XscIntegrand, 0., 1. );
176 176
177 result *= fCofXsc; 177 result *= fCofXsc;
178 178
179 result *= ZZ; // incoherent sum over all e 179 result *= ZZ; // incoherent sum over all element electrons
180 180
181 return result; 181 return result;
182 } 182 }
183 183
184 ///////////////////////////////////////// 184 /////////////////////////////////////////
185 // 185 //
186 // Rosenbluth relation in the neutron rest fra 186 // Rosenbluth relation in the neutron rest frame,
187 // x = sin^2(theta/2), theta is the electron s 187 // x = sin^2(theta/2), theta is the electron scattering angle
188 // Magnetic form factor in the dipole approxim 188 // Magnetic form factor in the dipole approximation.
189 189
190 G4double G4NeutronElectronElXsc::XscIntegrand( 190 G4double G4NeutronElectronElXsc::XscIntegrand(G4double x)
191 { 191 {
192 G4double result = 1., q2, znq2, znf, znf2, z 192 G4double result = 1., q2, znq2, znf, znf2, znf4;
193 193
194 znq2 = 1. + 2.*fee*x/fM; 194 znq2 = 1. + 2.*fee*x/fM;
195 195
196 q2 = 4.*fee2*x/znq2; 196 q2 = 4.*fee2*x/znq2;
197 197
198 znf = 1 + q2/fMv2; 198 znf = 1 + q2/fMv2;
199 znf2 = znf*znf; 199 znf2 = znf*znf;
200 znf4 = znf2*znf2; 200 znf4 = znf2*znf2;
201 201
202 result /= ( x + fAm )*znq2*znq2*znf4; 202 result /= ( x + fAm )*znq2*znq2*znf4;
203 203
204 result *= ( 1 - x )/( 1 + q2/4./fM2 ) + 2.*x 204 result *= ( 1 - x )/( 1 + q2/4./fM2 ) + 2.*x;
205 205
206 return result; 206 return result;
207 } 207 }
208 208
209 ////////////////////////////////////////////// 209 //////////////////////////////////////////////////////////
210 // 210 //
211 // Rosenbluth xsc in microbarn from 1*MeV to 1 211 // Rosenbluth xsc in microbarn from 1*MeV to 10*Tev, 200 points
212 212
213 const G4double G4NeutronElectronElXsc::fXscArr 213 const G4double G4NeutronElectronElXsc::fXscArray[200] = {
214 1.52681, 1.54903, 1.57123, 1.59341, 1.61556, 1 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 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 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 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 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 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 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 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 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. 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 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 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. 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 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 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. 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 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. 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 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 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. 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 235 1.34877, 1.30596, 1.26314, 1.22031, 1.17746, 1.13459, 1.0917, 1.04879, 1.00585,
236 0.962892, 0.919908 }; 236 0.962892, 0.919908 };
237 237