| 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 "G4TauNeutrinoNucleusTotXsc.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 G4TauNeutrinoNucleusTotXsc::G4TauNeutrinoNucle
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 G4double mt = 1.77686*GeV;
70 G4double mnp = 0.5*(proton_mass_c2+neutron_m
71 fEtc = mt + 0.5*mt*mt/mnp;
72 fDtc = fEtc - fEmc;
73 fIndex = 50;
74
75 fTotXsc = 0.;
76 fCcTotRatio = 0.75; // from nc/cc~0.33 ratio
77 fCcFactor = fNcFactor = 1.;
78 fQEratio = 0.5; // mean in the 1 GeV range
79
80 // theMuonMinus = G4MuonMinus::MuonMinus();
81 // theMuonPlus = G4MuonPlus::MuonPlus();
82 }
83
84 G4TauNeutrinoNucleusTotXsc::~G4TauNeutrinoNucl
85 {}
86
87 //////////////////////////////////////////////
88
89 G4bool
90 G4TauNeutrinoNucleusTotXsc::IsIsoApplicable( c
91 {
92 G4bool result = false;
93 G4String pName = aPart->GetDefinition()->Get
94 G4double tKin = aPart->GetKineticEnergy();
95
96 if( ( pName == "nu_tau" || pName == "
97 {
98 result = true;
99 }
100 return result;
101 }
102
103 //////////////////////////////////////
104
105 G4double G4TauNeutrinoNucleusTotXsc::GetElemen
106 G4int Z, const G4Material*
107 {
108 G4int Zi(0);
109 size_t i(0), j(0);
110 const G4ElementVector* theElementVector = ma
111
112 for ( i = 0; i < theElementVector->size(); +
113 {
114 Zi = (*theElementVector)[i]->GetZasInt();
115 if( Zi == Z ) break;
116 }
117 const G4Element* elm = (*theElementVector)[i
118 size_t nIso = elm->GetNumberOfIsotopes();
119 G4double fact = 0.0;
120 G4double xsec = 0.0;
121 const G4Isotope* iso = nullptr;
122 const G4IsotopeVector* isoVector = elm->GetI
123 const G4double* abundVector = elm->GetRelati
124
125 for (j = 0; j<nIso; ++j)
126 {
127 iso = (*isoVector)[j];
128 G4int A = iso->GetN();
129
130 if( abundVector[j] > 0.0 && IsIsoApplicabl
131 {
132 fact += abundVector[j];
133 xsec += abundVector[j]*GetIsoCrossSectio
134 }
135 }
136 if( fact > 0.0) { xsec /= fact; }
137 return xsec;
138 }
139
140 //////////////////////////////////////////////
141 //
142 //
143
144 G4double G4TauNeutrinoNucleusTotXsc::GetIsoCro
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 if( pName == "nu_tau" || pName == "ant_nu_ta
156
157 G4int index = GetEnergyIndex(energy);
158
159 if( index >= fIndex )
160 {
161 G4double pm = proton_mass_c2;
162 G4double s2 = 2.*energy*pm+pm*pm;
163 G4double aa = 1.;
164 G4double bb = 1.085;
165 G4double mw = 80.385*GeV;
166 fCcFactor = bb/(1.+ aa*s2/mw/mw);
167
168 G4double mz = 91.1876*GeV;
169 fNcFactor = bb/(1.+ aa*s2/mz/mz);
170 }
171 ccnuXsc = GetNuMuTotCsXsc(index, energy, Z,
172 ccnuXsc *= fCcFactor;
173 ccanuXsc = GetANuMuTotCsXsc(index, energy, Z
174 ccanuXsc *= fCcFactor;
175
176 if( pName == "nu_tau")
177 {
178 ncXsc = fCofL*ccnuXsc + fCofS*ccanuXsc;
179 ncXsc *= fNcFactor/fCcFactor;
180 totXsc = ccnuXsc + ncXsc;
181 if( totXsc > 0.) fCcTotRatio = ccnuXsc/tot
182 }
183 else if( pName == "anti_nu_tau")
184 {
185 ncXsc = fCofL*ccanuXsc + fCofS*ccnuXsc;
186 ncXsc *= fNcFactor/fCcFactor;
187 totXsc = ccanuXsc + ncXsc;
188 if( totXsc > 0.) fCcTotRatio = ccanuXsc/to
189 }
190 else return totXsc;
191
192 totXsc *= fCofXsc;
193 totXsc *= energy;
194 // totXsc *= A; // incoherent sum over all
195
196 totXsc *= fBiasingFactor; // biasing up, if
197
198 fTotXsc = totXsc;
199
200 return totXsc;
201 }
202
203 //////////////////////////////////////////////
204 //
205 // Return index of nu/anu energy array corresp
206
207 G4int G4TauNeutrinoNucleusTotXsc::GetEnergyInd
208 {
209 G4int i, eIndex = 0;
210
211 for( i = 0; i < fIndex; i++)
212 {
213 if( energy <= fNuMuEnergy[i]*GeV )
214 {
215 eIndex = i;
216 break;
217 }
218 }
219 if( i >= fIndex ) eIndex = i;
220 // G4cout<<"eIndex = "<<eIndex<<G4endl;
221 return eIndex;
222 }
223
224 //////////////////////////////////////////////
225 //
226 // nu_mu xsc for index-1, index linear over en
227
228 G4double G4TauNeutrinoNucleusTotXsc::GetNuMuTo
229 {
230 G4double xsc(0.), qexsc(0.), inxsc(0.);
231 G4int nn = aa - zz;
232 if(nn < 1) nn = 0;
233
234 // if( index <= 0 || energy < theMuonMinus->
235 if( index <= 0 || energy < fEmc ) xsc = aa*f
236 else if (index >= fIndex) xsc = aa*fNuMuInXs
237 else
238 {
239 G4double x1 = fNuMuEnergy[index-1]*GeV;
240 G4double x2 = fNuMuEnergy[index]*GeV;
241 G4double y1 = fNuMuInXsc[index-1];
242 G4double y2 = fNuMuInXsc[index];
243 G4double z1 = fNuMuQeXsc[index-1];
244 G4double z2 = fNuMuQeXsc[index];
245
246 if(x1 >= x2) return aa*fNuMuInXsc[index] +
247 else
248 {
249 G4double angle = (y2-y1)/(x2-x1);
250 inxsc = y1 + (energy-x1)*angle;
251 angle = (z2-z1)/(x2-x1);
252 qexsc = z1 + (energy-x1)*angle;
253 qexsc *= nn;
254 xsc = inxsc*aa + qexsc;
255
256 if( xsc > 0.) fQEratio = qexsc/xsc;
257 }
258 }
259 return xsc;
260 }
261
262 //////////////////////////////////////////////
263 //
264 // anu_mu xsc for index-1, index linear over e
265
266 G4double G4TauNeutrinoNucleusTotXsc::GetANuMuT
267 {
268 G4double xsc(0.), qexsc(0.), inxsc(0.);
269
270 // if( index <= 0 || energy < theMuonPlus->G
271 if( index <= 0 || energy < fEmc ) xsc = aa*f
272 else if (index >= fIndex) xsc = aa*fANuMuInX
273 else
274 {
275 G4double x1 = fNuMuEnergy[index-1]*GeV;
276 G4double x2 = fNuMuEnergy[index]*GeV;
277 G4double y1 = fANuMuInXsc[index-1];
278 G4double y2 = fANuMuInXsc[index];
279 G4double z1 = fANuMuQeXsc[index-1];
280 G4double z2 = fANuMuQeXsc[index];
281
282 if( x1 >= x2 ) return aa*fANuMuInXsc[index
283 else
284 {
285 G4double angle = (y2-y1)/(x2-x1);
286 inxsc = y1 + (energy-x1)*angle;
287
288 angle = (z2-z1)/(x2-x1);
289 qexsc = z1 + (energy-x1)*angle;
290 qexsc *= zz;
291 xsc = inxsc*aa + qexsc;
292
293 if( xsc > 0.) fQEratio = qexsc/xsc;
294 }
295 }
296 return xsc;
297 }
298
299 //////////////////////////////////////////////
300 //
301 // return fNuMuTotXsc[index] if the index is i
302
303 G4double G4TauNeutrinoNucleusTotXsc::GetNuMuTo
304 {
305 if( index >= 0 && index < fIndex) return fNu
306 else
307 {
308 G4cout<<"Improper index of fNuMuTotXsc arr
309 return 0.;
310 }
311 }
312
313 //////////////////////////////////////////////
314 //
315 // return fANuMuTotXsc[index] if the index is
316
317 G4double G4TauNeutrinoNucleusTotXsc::GetANuMuT
318 {
319 if( index >= 0 && index < fIndex) return fAN
320 else
321 {
322 G4cout<<"Improper index of fANuMuTotXsc ar
323 return 0.;
324 }
325 }
326
327
328 //////////////////////////////////////////////
329 //
330 // E_nu in GeV, ( Eth = 111.603 MeV, EthW = 33
331
332 const G4double G4TauNeutrinoNucleusTotXsc::fNu
333 {
334 0.12, 0.141136, 0.165996, 0.195233, 0.229621
335 0.270066, 0.317634, 0.373581, 0.439382, 0.51
336 0.607795, 0.714849, 0.84076, 0.988848, 1.163
337 1.36787, 1.6088, 1.89217, 2.22545, 2.61743,
338 3.07845, 3.62068, 4.25841, 5.00847, 5.89065,
339 6.9282, 8.14851, 9.58376, 11.2718, 13.2572,
340 15.5922, 18.3386, 21.5687, 25.3677, 29.8359,
341 35.0911, 41.2719, 48.5413, 57.0912, 67.147,
342 78.974, 92.8842, 109.244, 128.486, 151.117,
343 177.735, 209.04, 245.86, 289.164, 340.097 };
344
345 //////////////////////////////////////////////
346 //
347 // XS/E arrays in 10^-38cm2/GeV
348
349 const G4double G4TauNeutrinoNucleusTotXsc::fNu
350 {
351 0, 0, 0, 0, 0,
352 0, 0, 0.0166853, 0.0649693, 0.132346,
353 0.209102, 0.286795, 0.3595, 0.423961, 0.479
354 0.524797, 0.562165, 0.592225, 0.61612, 0.63
355 0.649524, 0.660751, 0.669245, 0.675546, 0.68
356 0.683247, 0.685307, 0.686521, 0.687093, 0.68
357 0.686919, 0.686384, 0.685631, 0.684689, 0.68
358 0.682275, 0.680806, 0.67917, 0.677376, 0.675
359 0.673387, 0.671229, 0.668985, 0.666665, 0.66
360 0.661804, 0.65925, 0.656593, 0.65381, 0.6508
361
362 const G4double G4TauNeutrinoNucleusTotXsc::fNu
363 {
364 0.20787, 0.411055, 0.570762, 0.705379, 0.814
365 0.89543, 0.944299, 0.959743, 0.942906, 0.897
366 0.831331, 0.750948, 0.66443, 0.578191, 0.496
367 0.423071, 0.358103, 0.302016, 0.254241, 0.21
368 0.179971, 0.151527, 0.12769, 0.107706, 0.090
369 0.0768491, 0.0649975, 0.0550143, 0.0465948,
370 0.0334782, 0.0283964, 0.0240945, 0.0204506,
371 0.0147437, 0.0125223, 0.0106374, 0.00903737,
372 0.00652531, 0.00554547, 0.0047131, 0.0040059
373 0.00289436, 0.00246039, 0.00209155, 0.001778
374
375
376
377 //////////////////////////////////////////////
378
379 const G4double G4TauNeutrinoNucleusTotXsc::fAN
380 {
381 0, 0, 0, 0, 0,
382 0, 0, 0.00437363, 0.0161485, 0.0333162,
383 0.0557621, 0.0814548, 0.108838, 0.136598, 0.
384 0.188908, 0.212041, 0.232727, 0.250872, 0.26
385 0.279467, 0.290341, 0.299177, 0.306299, 0.31
386 0.316108, 0.319378, 0.321892, 0.323583, 0.32
387 0.325841, 0.326568, 0.327111, 0.327623, 0.32
388 0.328412, 0.328704, 0.328988, 0.329326, 0.32
389 0.329791, 0.330051, 0.330327, 0.33057, 0.330
390 0.331115, 0.331416, 0.331678, 0.33192, 0.332
391
392 //////////////////////////////////////////////
393
394 const G4double G4TauNeutrinoNucleusTotXsc::fAN
395 {
396 0.0770264, 0.138754, 0.177006, 0.202417, 0.2
397 0.225742, 0.227151, 0.223805, 0.21709, 0.208
398 0.197763, 0.186496, 0.174651, 0.162429, 0.14
399 0.137498, 0.125127, 0.113057, 0.101455, 0.09
400 0.0801914, 0.0707075, 0.0620483, 0.0542192,
401 0.0409571, 0.0354377, 0.0305862, 0.0263422,
402 0.0194358, 0.0166585, 0.0142613, 0.0121968,
403 0.00889912, 0.00759389, 0.00647662, 0.005521
404 0.00400791, 0.00341322, 0.00290607, 0.002473
405 0.00179162, 0.00152441, 0.00129691, 0.001103
406
407 ////////////////////
408
409