| 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 "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::G4TauNeutrinoNucleusTotXsc()
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 G4double mt = 1.77686*GeV;
70 G4double mnp = 0.5*(proton_mass_c2+neutron_mass_c2);
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::~G4TauNeutrinoNucleusTotXsc()
85 {}
86
87 //////////////////////////////////////////////////////
88
89 G4bool
90 G4TauNeutrinoNucleusTotXsc::IsIsoApplicable( const G4DynamicParticle* aPart, G4int, G4int, const G4Element*, const G4Material*)
91 {
92 G4bool result = false;
93 G4String pName = aPart->GetDefinition()->GetParticleName();
94 G4double tKin = aPart->GetKineticEnergy();
95
96 if( ( pName == "nu_tau" || pName == "anti_nu_tau") && tKin >= fEtc )
97 {
98 result = true;
99 }
100 return result;
101 }
102
103 //////////////////////////////////////
104
105 G4double G4TauNeutrinoNucleusTotXsc::GetElementCrossSection(const G4DynamicParticle* part,
106 G4int Z, const G4Material* mat )
107 {
108 G4int Zi(0);
109 size_t i(0), j(0);
110 const G4ElementVector* theElementVector = mat->GetElementVector();
111
112 for ( i = 0; i < theElementVector->size(); ++i )
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->GetIsotopeVector();
123 const G4double* abundVector = elm->GetRelativeAbundanceVector();
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 && IsIsoApplicable(part, Z, A, elm, mat) )
131 {
132 fact += abundVector[j];
133 xsec += abundVector[j]*GetIsoCrossSection( part, Z, A, iso, elm, mat);
134 }
135 }
136 if( fact > 0.0) { xsec /= fact; }
137 return xsec;
138 }
139
140 ////////////////////////////////////////////////////
141 //
142 //
143
144 G4double G4TauNeutrinoNucleusTotXsc::GetIsoCrossSection(const G4DynamicParticle* aPart, G4int Z, 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 if( pName == "nu_tau" || pName == "ant_nu_tau" ) energy -= fDtc; // scaling energy for tau-neutrinos
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, A);
172 ccnuXsc *= fCcFactor;
173 ccanuXsc = GetANuMuTotCsXsc(index, energy, Z, A);
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/totXsc;
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/totXsc;
189 }
190 else return totXsc;
191
192 totXsc *= fCofXsc;
193 totXsc *= energy;
194 // totXsc *= A; // incoherent sum over all isotope nucleons
195
196 totXsc *= fBiasingFactor; // biasing up, if set >1
197
198 fTotXsc = totXsc;
199
200 return totXsc;
201 }
202
203 /////////////////////////////////////////////////////
204 //
205 // Return index of nu/anu energy array corresponding to the neutrino energy
206
207 G4int G4TauNeutrinoNucleusTotXsc::GetEnergyIndex(G4double energy)
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 energy
227
228 G4double G4TauNeutrinoNucleusTotXsc::GetNuMuTotCsXsc(G4int index, G4double energy, G4int zz, G4int aa)
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->GetPDGMass() ) xsc = aa*fNuMuInXsc[0] + nn*fNuMuQeXsc[0];
235 if( index <= 0 || energy < fEmc ) xsc = aa*fNuMuInXsc[0] + nn*fNuMuQeXsc[0];
236 else if (index >= fIndex) xsc = aa*fNuMuInXsc[fIndex-1] + nn*fNuMuQeXsc[fIndex-1];
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] + nn*fNuMuQeXsc[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 energy
265
266 G4double G4TauNeutrinoNucleusTotXsc::GetANuMuTotCsXsc(G4int index, G4double energy, G4int zz, G4int aa)
267 {
268 G4double xsc(0.), qexsc(0.), inxsc(0.);
269
270 // if( index <= 0 || energy < theMuonPlus->GetPDGMass() ) xsc = aa*fANuMuInXsc[0] + zz*fANuMuQeXsc[0];
271 if( index <= 0 || energy < fEmc ) xsc = aa*fANuMuInXsc[0] + zz*fANuMuQeXsc[0];
272 else if (index >= fIndex) xsc = aa*fANuMuInXsc[fIndex-1] + zz*fANuMuQeXsc[fIndex-1];
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] + zz*fANuMuQeXsc[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 in the array range
302
303 G4double G4TauNeutrinoNucleusTotXsc::GetNuMuTotCsArray( G4int index)
304 {
305 if( index >= 0 && index < fIndex) return fNuMuInXsc[index] + fNuMuQeXsc[index];
306 else
307 {
308 G4cout<<"Improper index of fNuMuTotXsc array"<<G4endl;
309 return 0.;
310 }
311 }
312
313 ////////////////////////////////////////////////////////
314 //
315 // return fANuMuTotXsc[index] if the index is in the array range
316
317 G4double G4TauNeutrinoNucleusTotXsc::GetANuMuTotCsArray( G4int index)
318 {
319 if( index >= 0 && index < fIndex) return fANuMuInXsc[index] + fANuMuQeXsc[index];
320 else
321 {
322 G4cout<<"Improper index of fANuMuTotXsc array"<<G4endl;
323 return 0.;
324 }
325 }
326
327
328 ///////////////////////////////////////////////////////
329 //
330 // E_nu in GeV, ( Eth = 111.603 MeV, EthW = 330.994 MeV)
331
332 const G4double G4TauNeutrinoNucleusTotXsc::fNuMuEnergy[50] =
333 {
334 0.12, 0.141136, 0.165996, 0.195233, 0.229621,
335 0.270066, 0.317634, 0.373581, 0.439382, 0.516773,
336 0.607795, 0.714849, 0.84076, 0.988848, 1.16302,
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::fNuMuInXsc[50] =
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.479009,
354 0.524797, 0.562165, 0.592225, 0.61612, 0.63491,
355 0.649524, 0.660751, 0.669245, 0.675546, 0.680092,
356 0.683247, 0.685307, 0.686521, 0.687093, 0.687184,
357 0.686919, 0.686384, 0.685631, 0.684689, 0.68357,
358 0.682275, 0.680806, 0.67917, 0.677376, 0.675442,
359 0.673387, 0.671229, 0.668985, 0.666665, 0.664272,
360 0.661804, 0.65925, 0.656593, 0.65381, 0.650871 };
361
362 const G4double G4TauNeutrinoNucleusTotXsc::fNuMuQeXsc[50] =
363 {
364 0.20787, 0.411055, 0.570762, 0.705379, 0.814702,
365 0.89543, 0.944299, 0.959743, 0.942906, 0.897917,
366 0.831331, 0.750948, 0.66443, 0.578191, 0.496828,
367 0.423071, 0.358103, 0.302016, 0.254241, 0.213889,
368 0.179971, 0.151527, 0.12769, 0.107706, 0.0909373,
369 0.0768491, 0.0649975, 0.0550143, 0.0465948, 0.0394861,
370 0.0334782, 0.0283964, 0.0240945, 0.0204506, 0.0173623,
371 0.0147437, 0.0125223, 0.0106374, 0.00903737, 0.00767892,
372 0.00652531, 0.00554547, 0.0047131, 0.0040059, 0.003405,
373 0.00289436, 0.00246039, 0.00209155, 0.00177804, 0.00151152 };
374
375
376
377 //////////////////////////////////////////////////////////////////
378
379 const G4double G4TauNeutrinoNucleusTotXsc::fANuMuInXsc[50] =
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.163526,
384 0.188908, 0.212041, 0.232727, 0.250872, 0.26631,
385 0.279467, 0.290341, 0.299177, 0.306299, 0.311864,
386 0.316108, 0.319378, 0.321892, 0.323583, 0.324909,
387 0.325841, 0.326568, 0.327111, 0.327623, 0.32798,
388 0.328412, 0.328704, 0.328988, 0.329326, 0.329559,
389 0.329791, 0.330051, 0.330327, 0.33057, 0.330834,
390 0.331115, 0.331416, 0.331678, 0.33192, 0.332124 };
391
392 //////////////////////////////////////////////////////////////////
393
394 const G4double G4TauNeutrinoNucleusTotXsc::fANuMuQeXsc[50] =
395 {
396 0.0770264, 0.138754, 0.177006, 0.202417, 0.21804,
397 0.225742, 0.227151, 0.223805, 0.21709, 0.208137,
398 0.197763, 0.186496, 0.174651, 0.162429, 0.14999,
399 0.137498, 0.125127, 0.113057, 0.101455, 0.0904642,
400 0.0801914, 0.0707075, 0.0620483, 0.0542192, 0.0472011,
401 0.0409571, 0.0354377, 0.0305862, 0.0263422, 0.0226451,
402 0.0194358, 0.0166585, 0.0142613, 0.0121968, 0.0104221,
403 0.00889912, 0.00759389, 0.00647662, 0.00552119, 0.00470487,
404 0.00400791, 0.00341322, 0.00290607, 0.00247377, 0.0021054,
405 0.00179162, 0.00152441, 0.00129691, 0.00110323, 0.000938345 };
406
407 ////////////////////
408
409