| 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
28 #include "G4MuNeutrinoNucleusTotXsc.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4DynamicParticle.hh"
32 #include "G4ParticleTable.hh"
33 #include "G4IonTable.hh"
34 #include "G4HadTmpUtil.hh"
35 #include "G4NistManager.hh"
36 #include "G4Material.hh"
37 #include "G4Element.hh"
38 #include "G4Isotope.hh"
39 #include "G4ElementVector.hh"
40
41 #include "G4MuonMinus.hh"
42 #include "G4MuonPlus.hh"
43
44 using namespace std;
45 using namespace CLHEP;
46
47 G4MuNeutrinoNucleusTotXsc::G4MuNeutrinoNucleusTotXsc()
48 : G4VCrossSectionDataSet("NuMuNuclTotXsc")
49 {
50 fCofXsc = 1.e-38*cm2/GeV;
51
52 // G4cout<<"fCofXsc = "<<fCofXsc*GeV/cm2<<" cm2/GeV"<<G4endl;
53
54 // PDG2016: sin^2 theta Weinberg
55
56 fSin2tW = 0.23129; // 0.2312;
57
58 // 9 <-> 6, 5/9 or 5/6 ?
59
60 fCofS = 5.*fSin2tW*fSin2tW/9.;
61 fCofL = 1. - fSin2tW + fCofS;
62
63 // G4cout<<"fCosL = "<<fCofL<<", fCofS = "<<fCofS<<G4endl;
64
65 fCutEnergy = 0.; // default value
66
67 fBiasingFactor = 1.; // default as physics
68 fEmc = 0.2*GeV;
69 fIndex = 50;
70
71 fTotXsc = 0.;
72 fCcTotRatio = 0.75; // from nc/cc~0.33 ratio
73 fCcFactor = fNcFactor = 1.;
74 fQEratio = 0.5; // mean in the 1 GeV range
75
76 // theMuonMinus = G4MuonMinus::MuonMinus();
77 // theMuonPlus = G4MuonPlus::MuonPlus();
78 }
79
80 G4MuNeutrinoNucleusTotXsc::~G4MuNeutrinoNucleusTotXsc()
81 {}
82
83 //////////////////////////////////////////////////////
84
85 G4bool
86 G4MuNeutrinoNucleusTotXsc::IsIsoApplicable( const G4DynamicParticle* aPart, G4int, G4int, const G4Element*, const G4Material*)
87 {
88 G4bool result = false;
89 G4String pName = aPart->GetDefinition()->GetParticleName();
90 G4double tKin = aPart->GetKineticEnergy();
91
92 if( ( pName == "nu_mu" || pName == "anti_nu_mu") && tKin >= fEmc )
93 {
94 result = true;
95 }
96 return result;
97 }
98
99 //////////////////////////////////////
100
101 G4double G4MuNeutrinoNucleusTotXsc::GetElementCrossSection(const G4DynamicParticle* part,
102 G4int Z, const G4Material* mat )
103 {
104 G4int Zi(0);
105 size_t i(0), j(0);
106 const G4ElementVector* theElementVector = mat->GetElementVector();
107
108 for ( i = 0; i < theElementVector->size(); ++i )
109 {
110 Zi = (*theElementVector)[i]->GetZasInt();
111 if( Zi == Z ) break;
112 }
113 const G4Element* elm = (*theElementVector)[i];
114 size_t nIso = elm->GetNumberOfIsotopes();
115 G4double fact = 0.0;
116 G4double xsec = 0.0;
117 const G4Isotope* iso = nullptr;
118 const G4IsotopeVector* isoVector = elm->GetIsotopeVector();
119 const G4double* abundVector = elm->GetRelativeAbundanceVector();
120
121 for (j = 0; j<nIso; ++j)
122 {
123 iso = (*isoVector)[j];
124 G4int A = iso->GetN();
125
126 if( abundVector[j] > 0.0 && IsIsoApplicable(part, Z, A, elm, mat) )
127 {
128 fact += abundVector[j];
129 xsec += abundVector[j]*GetIsoCrossSection( part, Z, A, iso, elm, mat);
130 }
131 }
132 if( fact > 0.0) { xsec /= fact; }
133 return xsec;
134 }
135
136 ////////////////////////////////////////////////////
137 //
138 //
139
140 G4double G4MuNeutrinoNucleusTotXsc::GetIsoCrossSection(const G4DynamicParticle* aPart, G4int Z, G4int A,
141 const G4Isotope*, const G4Element*, const G4Material* )
142 {
143 fCcFactor = fNcFactor = 1.;
144 fCcTotRatio = 0.25;
145
146 G4double ccnuXsc, ccanuXsc, ncXsc, totXsc(0.);
147
148 G4double energy = aPart->GetTotalEnergy();
149 G4String pName = aPart->GetDefinition()->GetParticleName();
150
151 G4int index = GetEnergyIndex(energy);
152
153 if( index >= fIndex )
154 {
155 G4double pm = proton_mass_c2;
156 G4double s2 = 2.*energy*pm+pm*pm;
157 G4double aa = 1.;
158 G4double bb = 1.085;
159 G4double mw = 80.385*GeV;
160 fCcFactor = bb/(1.+ aa*s2/mw/mw);
161
162 G4double mz = 91.1876*GeV;
163 fNcFactor = bb/(1.+ aa*s2/mz/mz);
164 }
165 ccnuXsc = GetNuMuTotCsXsc(index, energy, Z, A);
166 ccnuXsc *= fCcFactor;
167 ccanuXsc = GetANuMuTotCsXsc(index, energy, Z, A);
168 ccanuXsc *= fCcFactor;
169
170 if( pName == "nu_mu" )
171 {
172 ncXsc = fCofL*ccnuXsc + fCofS*ccanuXsc;
173 ncXsc *= fNcFactor/fCcFactor;
174 totXsc = ccnuXsc + ncXsc;
175 if( totXsc > 0.) fCcTotRatio = ccnuXsc/totXsc;
176 }
177 else if( pName == "anti_nu_mu" )
178 {
179 ncXsc = fCofL*ccanuXsc + fCofS*ccnuXsc;
180 ncXsc *= fNcFactor/fCcFactor;
181 totXsc = ccanuXsc + ncXsc;
182 if( totXsc > 0.) fCcTotRatio = ccanuXsc/totXsc;
183 }
184 else return totXsc;
185
186 totXsc *= fCofXsc;
187 totXsc *= energy;
188 // totXsc *= A; // incoherent sum over all isotope nucleons
189
190 totXsc *= fBiasingFactor; // biasing up, if set >1
191
192 fTotXsc = totXsc;
193
194 return totXsc;
195 }
196
197 /////////////////////////////////////////////////////
198 //
199 // Return index of nu/anu energy array corresponding to the neutrino energy
200
201 G4int G4MuNeutrinoNucleusTotXsc::GetEnergyIndex(G4double energy)
202 {
203 G4int i, eIndex = 0;
204
205 for( i = 0; i < fIndex; i++)
206 {
207 if( energy <= fNuMuEnergy[i]*GeV )
208 {
209 eIndex = i;
210 break;
211 }
212 }
213 if( i >= fIndex ) eIndex = i;
214 // G4cout<<"eIndex = "<<eIndex<<G4endl;
215 return eIndex;
216 }
217
218 /////////////////////////////////////////////////////
219 //
220 // nu_mu xsc for index-1, index linear over energy
221
222 G4double G4MuNeutrinoNucleusTotXsc::GetNuMuTotCsXsc(G4int index, G4double energy, G4int zz, G4int aa)
223 {
224 G4double xsc(0.), qexsc(0.), inxsc(0.);
225 G4int nn = aa - zz;
226 if(nn < 1) nn = 0;
227
228 // if( index <= 0 || energy < theMuonMinus->GetPDGMass() ) xsc = aa*fNuMuInXsc[0] + nn*fNuMuQeXsc[0];
229 if( index <= 0 || energy < fEmc ) xsc = aa*fNuMuInXsc[0] + nn*fNuMuQeXsc[0];
230 else if (index >= fIndex) xsc = aa*fNuMuInXsc[fIndex-1] + nn*fNuMuQeXsc[fIndex-1];
231 else
232 {
233 G4double x1 = fNuMuEnergy[index-1]*GeV;
234 G4double x2 = fNuMuEnergy[index]*GeV;
235 G4double y1 = fNuMuInXsc[index-1];
236 G4double y2 = fNuMuInXsc[index];
237 G4double z1 = fNuMuQeXsc[index-1];
238 G4double z2 = fNuMuQeXsc[index];
239
240 if(x1 >= x2) return aa*fNuMuInXsc[index] + nn*fNuMuQeXsc[index];
241 else
242 {
243 G4double angle = (y2-y1)/(x2-x1);
244 inxsc = y1 + (energy-x1)*angle;
245 angle = (z2-z1)/(x2-x1);
246 qexsc = z1 + (energy-x1)*angle;
247 qexsc *= nn;
248 xsc = inxsc*aa + qexsc;
249
250 if( xsc > 0.) fQEratio = qexsc/xsc;
251 }
252 }
253 return xsc;
254 }
255
256 /////////////////////////////////////////////////////
257 //
258 // anu_mu xsc for index-1, index linear over energy
259
260 G4double G4MuNeutrinoNucleusTotXsc::GetANuMuTotCsXsc(G4int index, G4double energy, G4int zz, G4int aa)
261 {
262 G4double xsc(0.), qexsc(0.), inxsc(0.);
263
264 // if( index <= 0 || energy < theMuonPlus->GetPDGMass() ) xsc = aa*fANuMuInXsc[0] + zz*fANuMuQeXsc[0];
265 if( index <= 0 || energy < fEmc ) xsc = aa*fANuMuInXsc[0] + zz*fANuMuQeXsc[0];
266 else if (index >= fIndex) xsc = aa*fANuMuInXsc[fIndex-1] + zz*fANuMuQeXsc[fIndex-1];
267 else
268 {
269 G4double x1 = fNuMuEnergy[index-1]*GeV;
270 G4double x2 = fNuMuEnergy[index]*GeV;
271 G4double y1 = fANuMuInXsc[index-1];
272 G4double y2 = fANuMuInXsc[index];
273 G4double z1 = fANuMuQeXsc[index-1];
274 G4double z2 = fANuMuQeXsc[index];
275
276 if( x1 >= x2 ) return aa*fANuMuInXsc[index] + zz*fANuMuQeXsc[index];
277 else
278 {
279 G4double angle = (y2-y1)/(x2-x1);
280 inxsc = y1 + (energy-x1)*angle;
281
282 angle = (z2-z1)/(x2-x1);
283 qexsc = z1 + (energy-x1)*angle;
284 qexsc *= zz;
285 xsc = inxsc*aa + qexsc;
286
287 if( xsc > 0.) fQEratio = qexsc/xsc;
288 }
289 }
290 return xsc;
291 }
292
293 ////////////////////////////////////////////////////////
294 //
295 // return fNuMuTotXsc[index] if the index is in the array range
296
297 G4double G4MuNeutrinoNucleusTotXsc::GetNuMuTotCsArray( G4int index)
298 {
299 if( index >= 0 && index < fIndex) return fNuMuInXsc[index] + fNuMuQeXsc[index];
300 else
301 {
302 G4cout<<"Improper index of fNuMuTotXsc array"<<G4endl;
303 return 0.;
304 }
305 }
306
307 ////////////////////////////////////////////////////////
308 //
309 // return fANuMuTotXsc[index] if the index is in the array range
310
311 G4double G4MuNeutrinoNucleusTotXsc::GetANuMuTotCsArray( G4int index)
312 {
313 if( index >= 0 && index < fIndex) return fANuMuInXsc[index] + fANuMuQeXsc[index];
314 else
315 {
316 G4cout<<"Improper index of fANuMuTotXsc array"<<G4endl;
317 return 0.;
318 }
319 }
320
321
322 ///////////////////////////////////////////////////////
323 //
324 // E_nu in GeV, ( Eth = 111.603 MeV, EthW = 330.994 MeV)
325
326 const G4double G4MuNeutrinoNucleusTotXsc::fNuMuEnergy[50] =
327 {
328 0.12, 0.141136, 0.165996, 0.195233, 0.229621,
329 0.270066, 0.317634, 0.373581, 0.439382, 0.516773,
330 0.607795, 0.714849, 0.84076, 0.988848, 1.16302,
331 1.36787, 1.6088, 1.89217, 2.22545, 2.61743,
332 3.07845, 3.62068, 4.25841, 5.00847, 5.89065,
333 6.9282, 8.14851, 9.58376, 11.2718, 13.2572,
334 15.5922, 18.3386, 21.5687, 25.3677, 29.8359,
335 35.0911, 41.2719, 48.5413, 57.0912, 67.147,
336 78.974, 92.8842, 109.244, 128.486, 151.117,
337 177.735, 209.04, 245.86, 289.164, 340.097 };
338
339 ////////////////////////////////////////////////////
340 //
341 // XS/E arrays in 10^-38cm2/GeV
342
343 const G4double G4MuNeutrinoNucleusTotXsc::fNuMuInXsc[50] =
344 {
345 0, 0, 0, 0, 0,
346 0, 0, 0.0166853, 0.0649693, 0.132346,
347 0.209102, 0.286795, 0.3595, 0.423961, 0.479009,
348 0.524797, 0.562165, 0.592225, 0.61612, 0.63491,
349 0.649524, 0.660751, 0.669245, 0.675546, 0.680092,
350 0.683247, 0.685307, 0.686521, 0.687093, 0.687184,
351 0.686919, 0.686384, 0.685631, 0.684689, 0.68357,
352 0.682275, 0.680806, 0.67917, 0.677376, 0.675442,
353 0.673387, 0.671229, 0.668985, 0.666665, 0.664272,
354 0.661804, 0.65925, 0.656593, 0.65381, 0.650871 };
355
356 const G4double G4MuNeutrinoNucleusTotXsc::fNuMuQeXsc[50] =
357 {
358 0.20787, 0.411055, 0.570762, 0.705379, 0.814702,
359 0.89543, 0.944299, 0.959743, 0.942906, 0.897917,
360 0.831331, 0.750948, 0.66443, 0.578191, 0.496828,
361 0.423071, 0.358103, 0.302016, 0.254241, 0.213889,
362 0.179971, 0.151527, 0.12769, 0.107706, 0.0909373,
363 0.0768491, 0.0649975, 0.0550143, 0.0465948, 0.0394861,
364 0.0334782, 0.0283964, 0.0240945, 0.0204506, 0.0173623,
365 0.0147437, 0.0125223, 0.0106374, 0.00903737, 0.00767892,
366 0.00652531, 0.00554547, 0.0047131, 0.0040059, 0.003405,
367 0.00289436, 0.00246039, 0.00209155, 0.00177804, 0.00151152 };
368
369
370
371 //////////////////////////////////////////////////////////////////
372
373 const G4double G4MuNeutrinoNucleusTotXsc::fANuMuInXsc[50] =
374 {
375 0, 0, 0, 0, 0,
376 0, 0, 0.00437363, 0.0161485, 0.0333162,
377 0.0557621, 0.0814548, 0.108838, 0.136598, 0.163526,
378 0.188908, 0.212041, 0.232727, 0.250872, 0.26631,
379 0.279467, 0.290341, 0.299177, 0.306299, 0.311864,
380 0.316108, 0.319378, 0.321892, 0.323583, 0.324909,
381 0.325841, 0.326568, 0.327111, 0.327623, 0.32798,
382 0.328412, 0.328704, 0.328988, 0.329326, 0.329559,
383 0.329791, 0.330051, 0.330327, 0.33057, 0.330834,
384 0.331115, 0.331416, 0.331678, 0.33192, 0.332124 };
385
386 //////////////////////////////////////////////////////////////////
387
388 const G4double G4MuNeutrinoNucleusTotXsc::fANuMuQeXsc[50] =
389 {
390 0.0770264, 0.138754, 0.177006, 0.202417, 0.21804,
391 0.225742, 0.227151, 0.223805, 0.21709, 0.208137,
392 0.197763, 0.186496, 0.174651, 0.162429, 0.14999,
393 0.137498, 0.125127, 0.113057, 0.101455, 0.0904642,
394 0.0801914, 0.0707075, 0.0620483, 0.0542192, 0.0472011,
395 0.0409571, 0.0354377, 0.0305862, 0.0263422, 0.0226451,
396 0.0194358, 0.0166585, 0.0142613, 0.0121968, 0.0104221,
397 0.00889912, 0.00759389, 0.00647662, 0.00552119, 0.00470487,
398 0.00400791, 0.00341322, 0.00290607, 0.00247377, 0.0021054,
399 0.00179162, 0.00152441, 0.00129691, 0.00110323, 0.000938345 };
400
401 ////////////////////
402
403