| 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 // 24.04.20 V. Grichine
27 //
28 // (nu_e,anti_nu_e)-nucleus xsc
29
30
31
32 #include "G4ElNeutrinoNucleusTotXsc.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4DynamicParticle.hh"
36 #include "G4ParticleTable.hh"
37 #include "G4IonTable.hh"
38 #include "G4HadTmpUtil.hh"
39 #include "G4NistManager.hh"
40 #include "G4Material.hh"
41 #include "G4Element.hh"
42 #include "G4Isotope.hh"
43 #include "G4ElementVector.hh"
44
45 #include "G4Electron.hh"
46 #include "G4Positron.hh"
47
48 using namespace std;
49 using namespace CLHEP;
50
51 G4ElNeutrinoNucleusTotXsc::G4ElNeutrinoNucleus
52 : G4VCrossSectionDataSet("NuElNuclTotXsc")
53 {
54 fCofXsc = 1.e-38*cm2/GeV;
55
56 // G4cout<<"fCofXsc = "<<fCofXsc*GeV/cm2<<"
57
58 // PDG2016: sin^2 theta Weinberg
59
60 fSin2tW = 0.23129; // 0.2312;
61
62 // 9 <-> 6, 5/9 or 5/6 ?
63
64 fCofS = 5.*fSin2tW*fSin2tW/9.;
65 fCofL = 1. - fSin2tW + fCofS;
66
67 // G4cout<<"fCosL = "<<fCofL<<", fCofS = "<<
68
69 fCutEnergy = 0.; // default value
70
71 fBiasingFactor = 1.; // default as physics
72
73 fIndex = 50;
74
75 fTotXsc = 0.;
76 fCcTotRatio = 0.75; // from nc/cc~0.33 ratio
77 fCcFactor = fNcFactor = 1.;
78
79 theElectron = G4Electron::Electron();
80 thePositron = G4Positron::Positron();
81 }
82
83 G4ElNeutrinoNucleusTotXsc::~G4ElNeutrinoNucleu
84 {}
85
86 //////////////////////////////////////////////
87
88
89 G4bool
90 G4ElNeutrinoNucleusTotXsc::IsIsoApplicable( co
91 {
92 G4bool result = false;
93 G4String pName = aPart->GetDefinition()->Get
94
95 if( pName == "nu_e" || pName == "anti
96 {
97 result = true;
98 }
99 return result;
100 }
101
102 //////////////////////////////////////
103
104 G4double G4ElNeutrinoNucleusTotXsc::GetElement
105 G4int Z, const G4Material*
106 {
107 G4int Zi(0);
108 size_t i(0), j(0);
109 const G4ElementVector* theElementVector = ma
110
111 for ( i = 0; i < theElementVector->size(); +
112 {
113 Zi = (*theElementVector)[i]->GetZasInt();
114 if( Zi == Z ) break;
115 }
116 const G4Element* elm = (*theElementVector)[i
117 size_t nIso = elm->GetNumberOfIsotopes();
118 G4double fact = 0.0;
119 G4double xsec = 0.0;
120 const G4Isotope* iso = nullptr;
121 const G4IsotopeVector* isoVector = elm->GetI
122 const G4double* abundVector = elm->GetRelati
123
124 for (j = 0; j<nIso; ++j)
125 {
126 iso = (*isoVector)[j];
127 G4int A = iso->GetN();
128
129 if( abundVector[j] > 0.0 && IsIsoApplicabl
130 {
131 fact += abundVector[j];
132 xsec += abundVector[j]*GetIsoCrossSectio
133 }
134 }
135 if( fact > 0.0) { xsec /= fact; }
136 return xsec;
137 }
138
139
140 //////////////////////////////////////////////
141 //
142 //
143
144 G4double G4ElNeutrinoNucleusTotXsc::GetIsoCros
145 const G4Isotope*, const G4Element*
146 {
147 fCcFactor = fNcFactor = 1.;
148 fCcTotRatio = 0.25;
149
150 G4double ccnuXsc, ccanuXsc, ncXsc, totXsc(0.
151
152 G4double energy = aPart->GetTotalEnergy();
153 G4String pName = aPart->GetDefinition()->G
154
155 G4int index = GetEnergyIndex(energy);
156
157 if( index >= fIndex )
158 {
159 G4double pm = proton_mass_c2;
160 G4double s2 = 2.*energy*pm+pm*pm;
161 G4double aa = 1.;
162 G4double bb = 1.085;
163 G4double mw = 80.385*GeV;
164 fCcFactor = bb/(1.+ aa*s2/mw/mw);
165
166 G4double mz = 91.1876*GeV;
167 fNcFactor = bb/(1.+ aa*s2/mz/mz);
168 }
169 ccnuXsc = GetNuElTotCsXsc(index, energy);
170 ccnuXsc *= fCcFactor;
171 ccanuXsc = GetANuElTotCsXsc(index, energy);
172 ccanuXsc *= fCcFactor;
173
174 if( pName == "nu_e")
175 {
176 ncXsc = fCofL*ccnuXsc + fCofS*ccanuXsc;
177 ncXsc *= fNcFactor/fCcFactor;
178 totXsc = ccnuXsc + ncXsc;
179 if( totXsc > 0.) fCcTotRatio = ccnuXsc/tot
180 }
181 else if( pName == "anti_nu_e")
182 {
183 ncXsc = fCofL*ccanuXsc + fCofS*ccnuXsc;
184 ncXsc *= fNcFactor/fCcFactor;
185 totXsc = ccanuXsc + ncXsc;
186 if( totXsc > 0.) fCcTotRatio = ccanuXsc/to
187 }
188 else return totXsc;
189
190 totXsc *= fCofXsc;
191 totXsc *= energy;
192 totXsc *= A; // incoherent sum over all is
193
194 totXsc *= fBiasingFactor; // biasing up, if
195
196 fTotXsc = totXsc;
197
198 return totXsc;
199 }
200
201 //////////////////////////////////////////////
202 //
203 // Return index of nu/anu energy array corresp
204
205 G4int G4ElNeutrinoNucleusTotXsc::GetEnergyInde
206 {
207 G4int i, eIndex = 0;
208
209 for( i = 0; i < fIndex; i++)
210 {
211 if( energy <= fNuElEnergy[i]*GeV )
212 {
213 eIndex = i;
214 break;
215 }
216 }
217 if( i >= fIndex-1 ) eIndex = fIndex-1;
218 // G4cout<<"eIndex = "<<eIndex<<G4endl;
219 return eIndex;
220 }
221
222 //////////////////////////////////////////////
223 //
224 // nu_e xsc for index-1, index linear over ene
225
226 G4double G4ElNeutrinoNucleusTotXsc::GetNuElTot
227 {
228 G4double xsc(0.);
229
230 if( index <= 0 || energy < theElectron->GetP
231 else if (index >= fIndex) xsc = fNuElTotXsc[
232 else
233 {
234 G4double x1 = fNuElEnergy[index-1]*GeV;
235 G4double x2 = fNuElEnergy[index]*GeV;
236 G4double y1 = fNuElTotXsc[index-1];
237 G4double y2 = fNuElTotXsc[index];
238
239 if(x1 >= x2) return fNuElTotXsc[index];
240 else
241 {
242 G4double angle = (y2-y1)/(x2-x1);
243 xsc = y1 + (energy-x1)*angle;
244 }
245 }
246 return xsc;
247 }
248
249 //////////////////////////////////////////////
250 //
251 // anu_e xsc for index-1, index linear over en
252
253 G4double G4ElNeutrinoNucleusTotXsc::GetANuElTo
254 {
255 G4double xsc(0.);
256
257 if( index <= 0 || energy < thePositron->GetP
258 else if (index >= fIndex) xsc = fANuElTotXsc
259 else
260 {
261 G4double x1 = fNuElEnergy[index-1]*GeV;
262 G4double x2 = fNuElEnergy[index]*GeV;
263 G4double y1 = fANuElTotXsc[index-1];
264 G4double y2 = fANuElTotXsc[index];
265
266 if( x1 >= x2 ) return fANuElTotXsc[index];
267 else
268 {
269 G4double angle = (y2-y1)/(x2-x1);
270 xsc = y1 + (energy-x1)*angle;
271 }
272 }
273 return xsc;
274 }
275
276 //////////////////////////////////////////////
277 //
278 // return fNuElTotXsc[index] if the index is i
279
280 G4double G4ElNeutrinoNucleusTotXsc::GetNuElTot
281 {
282 if( index >= 0 && index < fIndex) return fNu
283 else
284 {
285 G4cout<<"Improper index of fNuElTotXsc arr
286 return 0.;
287 }
288 }
289
290 //////////////////////////////////////////////
291 //
292 // return fANuElTotXsc[index] if the index is
293
294 G4double G4ElNeutrinoNucleusTotXsc::GetANuElTo
295 {
296 if( index >= 0 && index < fIndex) return fAN
297 else
298 {
299 G4cout<<"Improper index of fANuElTotXsc ar
300 return 0.;
301 }
302 }
303
304
305 //////////////////////////////////////////////
306 //
307 // E_nu in GeV
308
309 const G4double G4ElNeutrinoNucleusTotXsc::fNuE
310 {
311 0.000561138, 0.000735091, 0.000962969, 0.001
312 0.00216484, 0.00283594, 0.00371508, 0.004866
313 0.00835185, 0.0109409, 0.0143326, 0.0187757,
314 0.032221, 0.0422095, 0.0552945, 0.0724358, 0
315 0.124307, 0.162842, 0.213323, 0.279453, 0.36
316 0.47957, 0.628237, 0.82299, 1.07812, 1.41233
317 1.85016, 2.42371, 3.17505, 4.15932, 5.44871,
318 7.13781, 9.35053, 12.2492, 16.0464, 21.0208,
319 27.5373, 36.0739, 47.2568, 61.9064, 81.0973,
320 106.238, 139.171, 182.314, 238.832, 312.869
321 };
322
323 //////////////////////////////////////////////
324 //
325 // nu_e CC xsc_tot/E_nu, in 10^-38 cm2/GeV
326
327 const G4double G4ElNeutrinoNucleusTotXsc::fNuE
328 {
329 0.0026484, 0.00609503, 0.00939421, 0.0132163
330 0.0237692, 0.0312066, 0.0406632, 0.0526867,
331 0.0871913, 0.111359, 0.141458, 0.178584, 0.2
332 0.27822, 0.342461, 0.416865, 0.501361, 0.596
333 0.713623, 0.905749, 1.20718, 1.52521, 1.7528
334 1.82072, 1.67119, 1.50074, 1.3077, 1.14923,
335 1.0577, 0.977911, 0.918526, 0.792889, 0.7022
336 0.678615, 0.687099, 0.725167, 0.706795, 0.67
337 0.649791, 0.651328, 0.651934, 0.658062, 0.66
338 0.662534, 0.662601, 0.660261, 0.656724, 0.65
339 };
340
341
342
343 //////////////////////////////////////////////
344 //
345 // anu_e CC xsc_tot/E_nu, in 10^-38 cm2/GeV
346
347 const G4double G4ElNeutrinoNucleusTotXsc::fANu
348 {
349 0.00103385, 0.00237807, 0.00366358, 0.005151
350 0.00925859, 0.0121508, 0.0158252, 0.0204908,
351 0.0338304, 0.0431234, 0.0546346, 0.068735, 0
352 0.106025, 0.129614, 0.15643, 0.186063, 0.217
353 0.251065, 0.28525, 0.319171, 0.348995, 0.369
354 0.378165, 0.377353, 0.371224, 0.363257, 0.35
355 0.348618, 0.343082, 0.338825, 0.33574, 0.333
356 0.332504, 0.332052, 0.332187, 0.332781, 0.33
357 0.33489, 0.336213, 0.337608, 0.339008, 0.340
358 0.341606, 0.342706, 0.343628, 0.344305, 0.34
359 };
360