| 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 // 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::G4ElNeutrinoNucleusTotXsc()
52 : G4VCrossSectionDataSet("NuElNuclTotXsc")
53 {
54 fCofXsc = 1.e-38*cm2/GeV;
55
56 // G4cout<<"fCofXsc = "<<fCofXsc*GeV/cm2<<" cm2/GeV"<<G4endl;
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 = "<<fCofS<<G4endl;
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::~G4ElNeutrinoNucleusTotXsc()
84 {}
85
86 //////////////////////////////////////////////////////
87
88
89 G4bool
90 G4ElNeutrinoNucleusTotXsc::IsIsoApplicable( const G4DynamicParticle* aPart, G4int, G4int, const G4Element*, const G4Material*)
91 {
92 G4bool result = false;
93 G4String pName = aPart->GetDefinition()->GetParticleName();
94
95 if( pName == "nu_e" || pName == "anti_nu_e" )
96 {
97 result = true;
98 }
99 return result;
100 }
101
102 //////////////////////////////////////
103
104 G4double G4ElNeutrinoNucleusTotXsc::GetElementCrossSection(const G4DynamicParticle* part,
105 G4int Z, const G4Material* mat )
106 {
107 G4int Zi(0);
108 size_t i(0), j(0);
109 const G4ElementVector* theElementVector = mat->GetElementVector();
110
111 for ( i = 0; i < theElementVector->size(); ++i )
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->GetIsotopeVector();
122 const G4double* abundVector = elm->GetRelativeAbundanceVector();
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 && IsIsoApplicable(part, Z, A, elm, mat) )
130 {
131 fact += abundVector[j];
132 xsec += abundVector[j]*GetIsoCrossSection( part, Z, A, iso, elm, mat);
133 }
134 }
135 if( fact > 0.0) { xsec /= fact; }
136 return xsec;
137 }
138
139
140 ////////////////////////////////////////////////////
141 //
142 //
143
144 G4double G4ElNeutrinoNucleusTotXsc::GetIsoCrossSection(const G4DynamicParticle* aPart, G4int, G4int A,
145 const G4Isotope*, const G4Element*, const G4Material* )
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()->GetParticleName();
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/totXsc;
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/totXsc;
187 }
188 else return totXsc;
189
190 totXsc *= fCofXsc;
191 totXsc *= energy;
192 totXsc *= A; // incoherent sum over all isotope nucleons
193
194 totXsc *= fBiasingFactor; // biasing up, if set >1
195
196 fTotXsc = totXsc;
197
198 return totXsc;
199 }
200
201 /////////////////////////////////////////////////////
202 //
203 // Return index of nu/anu energy array corresponding to the neutrino energy
204
205 G4int G4ElNeutrinoNucleusTotXsc::GetEnergyIndex(G4double energy)
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 energy
225
226 G4double G4ElNeutrinoNucleusTotXsc::GetNuElTotCsXsc(G4int index, G4double energy)
227 {
228 G4double xsc(0.);
229
230 if( index <= 0 || energy < theElectron->GetPDGMass() ) xsc = fNuElTotXsc[0];
231 else if (index >= fIndex) xsc = fNuElTotXsc[fIndex-1];
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 energy
252
253 G4double G4ElNeutrinoNucleusTotXsc::GetANuElTotCsXsc(G4int index, G4double energy)
254 {
255 G4double xsc(0.);
256
257 if( index <= 0 || energy < thePositron->GetPDGMass() ) xsc = fANuElTotXsc[0];
258 else if (index >= fIndex) xsc = fANuElTotXsc[fIndex-1];
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 in the array range
279
280 G4double G4ElNeutrinoNucleusTotXsc::GetNuElTotCsArray( G4int index)
281 {
282 if( index >= 0 && index < fIndex) return fNuElTotXsc[index];
283 else
284 {
285 G4cout<<"Improper index of fNuElTotXsc array"<<G4endl;
286 return 0.;
287 }
288 }
289
290 ////////////////////////////////////////////////////////
291 //
292 // return fANuElTotXsc[index] if the index is in the array range
293
294 G4double G4ElNeutrinoNucleusTotXsc::GetANuElTotCsArray( G4int index)
295 {
296 if( index >= 0 && index < fIndex) return fANuElTotXsc[index];
297 else
298 {
299 G4cout<<"Improper index of fANuElTotXsc array"<<G4endl;
300 return 0.;
301 }
302 }
303
304
305 ///////////////////////////////////////////////////////
306 //
307 // E_nu in GeV
308
309 const G4double G4ElNeutrinoNucleusTotXsc::fNuElEnergy[50] =
310 {
311 0.000561138, 0.000735091, 0.000962969, 0.00126149, 0.00165255,
312 0.00216484, 0.00283594, 0.00371508, 0.00486676, 0.00637546,
313 0.00835185, 0.0109409, 0.0143326, 0.0187757, 0.0245962,
314 0.032221, 0.0422095, 0.0552945, 0.0724358, 0.0948908,
315 0.124307, 0.162842, 0.213323, 0.279453, 0.366084,
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::fNuElTotXsc[50] =
328 {
329 0.0026484, 0.00609503, 0.00939421, 0.0132163, 0.0178983,
330 0.0237692, 0.0312066, 0.0406632, 0.0526867, 0.0679357,
331 0.0871913, 0.111359, 0.141458, 0.178584, 0.223838,
332 0.27822, 0.342461, 0.416865, 0.501361, 0.596739,
333 0.713623, 0.905749, 1.20718, 1.52521, 1.75286,
334 1.82072, 1.67119, 1.50074, 1.3077, 1.14923,
335 1.0577, 0.977911, 0.918526, 0.792889, 0.702282,
336 0.678615, 0.687099, 0.725167, 0.706795, 0.678045,
337 0.649791, 0.651328, 0.651934, 0.658062, 0.660659,
338 0.662534, 0.662601, 0.660261, 0.656724, 0.65212
339 };
340
341
342
343 /////////////////////////////////////////////////////////////
344 //
345 // anu_e CC xsc_tot/E_nu, in 10^-38 cm2/GeV
346
347 const G4double G4ElNeutrinoNucleusTotXsc::fANuElTotXsc[50] =
348 {
349 0.00103385, 0.00237807, 0.00366358, 0.00515192, 0.00697434,
350 0.00925859, 0.0121508, 0.0158252, 0.0204908, 0.0263959,
351 0.0338304, 0.0431234, 0.0546346, 0.068735, 0.0857738,
352 0.106025, 0.129614, 0.15643, 0.186063, 0.21784,
353 0.251065, 0.28525, 0.319171, 0.348995, 0.369448,
354 0.378165, 0.377353, 0.371224, 0.363257, 0.355433,
355 0.348618, 0.343082, 0.338825, 0.33574, 0.333684,
356 0.332504, 0.332052, 0.332187, 0.332781, 0.333716,
357 0.33489, 0.336213, 0.337608, 0.339008, 0.340362,
358 0.341606, 0.342706, 0.343628, 0.344305, 0.344675
359 };
360