| 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
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::G4MuNeutrinoNucleus
48 : G4VCrossSectionDataSet("NuMuNuclTotXsc")
49 {
50 fCofXsc = 1.e-38*cm2/GeV;
51
52 // G4cout<<"fCofXsc = "<<fCofXsc*GeV/cm2<<"
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 = "<<
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::~G4MuNeutrinoNucleu
81 {}
82
83 //////////////////////////////////////////////
84
85 G4bool
86 G4MuNeutrinoNucleusTotXsc::IsIsoApplicable( co
87 {
88 G4bool result = false;
89 G4String pName = aPart->GetDefinition()->Get
90 G4double tKin = aPart->GetKineticEnergy();
91
92 if( ( pName == "nu_mu" || pName == "a
93 {
94 result = true;
95 }
96 return result;
97 }
98
99 //////////////////////////////////////
100
101 G4double G4MuNeutrinoNucleusTotXsc::GetElement
102 G4int Z, const G4Material*
103 {
104 G4int Zi(0);
105 size_t i(0), j(0);
106 const G4ElementVector* theElementVector = ma
107
108 for ( i = 0; i < theElementVector->size(); +
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->GetI
119 const G4double* abundVector = elm->GetRelati
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 && IsIsoApplicabl
127 {
128 fact += abundVector[j];
129 xsec += abundVector[j]*GetIsoCrossSectio
130 }
131 }
132 if( fact > 0.0) { xsec /= fact; }
133 return xsec;
134 }
135
136 //////////////////////////////////////////////
137 //
138 //
139
140 G4double G4MuNeutrinoNucleusTotXsc::GetIsoCros
141 const G4Isotope*, const G4Element*
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()->G
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,
166 ccnuXsc *= fCcFactor;
167 ccanuXsc = GetANuMuTotCsXsc(index, energy, Z
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/tot
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/to
183 }
184 else return totXsc;
185
186 totXsc *= fCofXsc;
187 totXsc *= energy;
188 // totXsc *= A; // incoherent sum over all
189
190 totXsc *= fBiasingFactor; // biasing up, if
191
192 fTotXsc = totXsc;
193
194 return totXsc;
195 }
196
197 //////////////////////////////////////////////
198 //
199 // Return index of nu/anu energy array corresp
200
201 G4int G4MuNeutrinoNucleusTotXsc::GetEnergyInde
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 en
221
222 G4double G4MuNeutrinoNucleusTotXsc::GetNuMuTot
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->
229 if( index <= 0 || energy < fEmc ) xsc = aa*f
230 else if (index >= fIndex) xsc = aa*fNuMuInXs
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] +
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 e
259
260 G4double G4MuNeutrinoNucleusTotXsc::GetANuMuTo
261 {
262 G4double xsc(0.), qexsc(0.), inxsc(0.);
263
264 // if( index <= 0 || energy < theMuonPlus->G
265 if( index <= 0 || energy < fEmc ) xsc = aa*f
266 else if (index >= fIndex) xsc = aa*fANuMuInX
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
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 i
296
297 G4double G4MuNeutrinoNucleusTotXsc::GetNuMuTot
298 {
299 if( index >= 0 && index < fIndex) return fNu
300 else
301 {
302 G4cout<<"Improper index of fNuMuTotXsc arr
303 return 0.;
304 }
305 }
306
307 //////////////////////////////////////////////
308 //
309 // return fANuMuTotXsc[index] if the index is
310
311 G4double G4MuNeutrinoNucleusTotXsc::GetANuMuTo
312 {
313 if( index >= 0 && index < fIndex) return fAN
314 else
315 {
316 G4cout<<"Improper index of fANuMuTotXsc ar
317 return 0.;
318 }
319 }
320
321
322 //////////////////////////////////////////////
323 //
324 // E_nu in GeV, ( Eth = 111.603 MeV, EthW = 33
325
326 const G4double G4MuNeutrinoNucleusTotXsc::fNuM
327 {
328 0.12, 0.141136, 0.165996, 0.195233, 0.229621
329 0.270066, 0.317634, 0.373581, 0.439382, 0.51
330 0.607795, 0.714849, 0.84076, 0.988848, 1.163
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::fNuM
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.479
348 0.524797, 0.562165, 0.592225, 0.61612, 0.63
349 0.649524, 0.660751, 0.669245, 0.675546, 0.68
350 0.683247, 0.685307, 0.686521, 0.687093, 0.68
351 0.686919, 0.686384, 0.685631, 0.684689, 0.68
352 0.682275, 0.680806, 0.67917, 0.677376, 0.675
353 0.673387, 0.671229, 0.668985, 0.666665, 0.66
354 0.661804, 0.65925, 0.656593, 0.65381, 0.6508
355
356 const G4double G4MuNeutrinoNucleusTotXsc::fNuM
357 {
358 0.20787, 0.411055, 0.570762, 0.705379, 0.814
359 0.89543, 0.944299, 0.959743, 0.942906, 0.897
360 0.831331, 0.750948, 0.66443, 0.578191, 0.496
361 0.423071, 0.358103, 0.302016, 0.254241, 0.21
362 0.179971, 0.151527, 0.12769, 0.107706, 0.090
363 0.0768491, 0.0649975, 0.0550143, 0.0465948,
364 0.0334782, 0.0283964, 0.0240945, 0.0204506,
365 0.0147437, 0.0125223, 0.0106374, 0.00903737,
366 0.00652531, 0.00554547, 0.0047131, 0.0040059
367 0.00289436, 0.00246039, 0.00209155, 0.001778
368
369
370
371 //////////////////////////////////////////////
372
373 const G4double G4MuNeutrinoNucleusTotXsc::fANu
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.
378 0.188908, 0.212041, 0.232727, 0.250872, 0.26
379 0.279467, 0.290341, 0.299177, 0.306299, 0.31
380 0.316108, 0.319378, 0.321892, 0.323583, 0.32
381 0.325841, 0.326568, 0.327111, 0.327623, 0.32
382 0.328412, 0.328704, 0.328988, 0.329326, 0.32
383 0.329791, 0.330051, 0.330327, 0.33057, 0.330
384 0.331115, 0.331416, 0.331678, 0.33192, 0.332
385
386 //////////////////////////////////////////////
387
388 const G4double G4MuNeutrinoNucleusTotXsc::fANu
389 {
390 0.0770264, 0.138754, 0.177006, 0.202417, 0.2
391 0.225742, 0.227151, 0.223805, 0.21709, 0.208
392 0.197763, 0.186496, 0.174651, 0.162429, 0.14
393 0.137498, 0.125127, 0.113057, 0.101455, 0.09
394 0.0801914, 0.0707075, 0.0620483, 0.0542192,
395 0.0409571, 0.0354377, 0.0305862, 0.0263422,
396 0.0194358, 0.0166585, 0.0142613, 0.0121968,
397 0.00889912, 0.00759389, 0.00647662, 0.005521
398 0.00400791, 0.00341322, 0.00290607, 0.002473
399 0.00179162, 0.00152441, 0.00129691, 0.001103
400
401 ////////////////////
402
403