ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuNeutrinoNucleusTotXsc.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuNeutrinoNucleusTotXsc.cc
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 
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 
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 
69  fIndex = 50;
70 
71  fTotXsc = 0.;
72  fCcTotRatio = 0.75; // from nc/cc~0.33 ratio
73  fCcFactor = fNcFactor = 1.;
74 
77 }
78 
80 {}
81 
83 
84 G4bool
86 {
87  G4bool result = false;
88  G4String pName = aPart->GetDefinition()->GetParticleName();
89 
90  if( pName == "nu_mu" || pName == "anti_nu_mu" )
91  {
92  result = true;
93  }
94  return result;
95 }
96 
98 
100  G4int Z, const G4Material* mat )
101 {
102  G4int Zi(0);
103  size_t i(0), j(0);
104  const G4ElementVector* theElementVector = mat->GetElementVector();
105 
106  for ( i = 0; i < theElementVector->size(); ++i )
107  {
108  Zi = (*theElementVector)[i]->GetZasInt();
109  if( Zi == Z ) break;
110  }
111  const G4Element* elm = (*theElementVector)[i];
112  size_t nIso = elm->GetNumberOfIsotopes();
113  G4double fact = 0.0;
114  G4double xsec = 0.0;
115  const G4Isotope* iso = nullptr;
116  const G4IsotopeVector* isoVector = elm->GetIsotopeVector();
117  const G4double* abundVector = elm->GetRelativeAbundanceVector();
118 
119  for (j = 0; j<nIso; ++j)
120  {
121  iso = (*isoVector)[j];
122  G4int A = iso->GetN();
123 
124  if( abundVector[j] > 0.0 && IsIsoApplicable(part, Z, A, elm, mat) )
125  {
126  fact += abundVector[j];
127  xsec += abundVector[j]*GetIsoCrossSection( part, Z, A, iso, elm, mat);
128  }
129  }
130  if( fact > 0.0) { xsec /= fact; }
131  return xsec;
132 }
133 
135 //
136 //
137 
139  const G4Isotope*, const G4Element*, const G4Material* )
140 {
141  fCcFactor = fNcFactor = 1.;
142  fCcTotRatio = 0.25;
143 
144  G4double ccnuXsc, ccanuXsc, ncXsc, totXsc(0.);
145 
146  G4double energy = aPart->GetTotalEnergy();
147  G4String pName = aPart->GetDefinition()->GetParticleName();
148 
149  G4int index = GetEnergyIndex(energy);
150 
151  if( index >= fIndex )
152  {
154  G4double s2 = 2.*energy*pm+pm*pm;
155  G4double aa = 1.;
156  G4double bb = 1.085;
157  G4double mw = 80.385*GeV;
158  fCcFactor = bb/(1.+ aa*s2/mw/mw);
159 
160  G4double mz = 91.1876*GeV;
161  fNcFactor = bb/(1.+ aa*s2/mz/mz);
162  }
163  ccnuXsc = GetNuMuTotCsXsc(index, energy);
164  ccnuXsc *= fCcFactor;
165  ccanuXsc = GetANuMuTotCsXsc(index, energy);
166  ccanuXsc *= fCcFactor;
167 
168  if( pName == "nu_mu")
169  {
170  ncXsc = fCofL*ccnuXsc + fCofS*ccanuXsc;
171  ncXsc *= fNcFactor/fCcFactor;
172  totXsc = ccnuXsc + ncXsc;
173  if( totXsc > 0.) fCcTotRatio = ccnuXsc/totXsc;
174  }
175  else if( pName == "anti_nu_mu")
176  {
177  ncXsc = fCofL*ccanuXsc + fCofS*ccnuXsc;
178  ncXsc *= fNcFactor/fCcFactor;
179  totXsc = ccanuXsc + ncXsc;
180  if( totXsc > 0.) fCcTotRatio = ccanuXsc/totXsc;
181  }
182  else return totXsc;
183 
184  totXsc *= fCofXsc;
185  totXsc *= energy;
186  totXsc *= A; // incoherent sum over all isotope nucleons
187 
188  totXsc *= fBiasingFactor; // biasing up, if set >1
189 
190  fTotXsc = totXsc;
191 
192  return totXsc;
193 }
194 
196 //
197 // Return index of nu/anu energy array corresponding to the neutrino energy
198 
200 {
201  G4int i, eIndex = 0;
202 
203  for( i = 0; i < fIndex; i++)
204  {
205  if( energy <= fNuMuEnergy[i]*GeV )
206  {
207  eIndex = i;
208  break;
209  }
210  }
211  if( i >= fIndex ) eIndex = i;
212  // G4cout<<"eIndex = "<<eIndex<<G4endl;
213  return eIndex;
214 }
215 
217 //
218 // nu_mu xsc for index-1, index linear over energy
219 
221 {
222  G4double xsc(0.);
223 
224  if( index <= 0 || energy < theMuonMinus->GetPDGMass() ) xsc = fNuMuTotXsc[0];
225  else if (index >= fIndex) xsc = fNuMuTotXsc[fIndex-1];
226  else
227  {
228  G4double x1 = fNuMuEnergy[index-1]*GeV;
229  G4double x2 = fNuMuEnergy[index]*GeV;
230  G4double y1 = fNuMuTotXsc[index-1];
231  G4double y2 = fNuMuTotXsc[index];
232 
233  if(x1 >= x2) return fNuMuTotXsc[index];
234  else
235  {
236  G4double angle = (y2-y1)/(x2-x1);
237  xsc = y1 + (energy-x1)*angle;
238  }
239  }
240  return xsc;
241 }
242 
244 //
245 // anu_mu xsc for index-1, index linear over energy
246 
248 {
249  G4double xsc(0.);
250 
251  if( index <= 0 || energy < theMuonPlus->GetPDGMass() ) xsc = fANuMuTotXsc[0];
252  else if (index >= fIndex) xsc = fANuMuTotXsc[fIndex-1];
253  else
254  {
255  G4double x1 = fNuMuEnergy[index-1]*GeV;
256  G4double x2 = fNuMuEnergy[index]*GeV;
257  G4double y1 = fANuMuTotXsc[index-1];
258  G4double y2 = fANuMuTotXsc[index];
259 
260  if( x1 >= x2 ) return fANuMuTotXsc[index];
261  else
262  {
263  G4double angle = (y2-y1)/(x2-x1);
264  xsc = y1 + (energy-x1)*angle;
265  }
266  }
267  return xsc;
268 }
269 
271 //
272 // return fNuMuTotXsc[index] if the index is in the array range
273 
275 {
276  if( index >= 0 && index < fIndex) return fNuMuTotXsc[index];
277  else
278  {
279  G4cout<<"Inproper index of fNuMuTotXsc array"<<G4endl;
280  return 0.;
281  }
282 }
283 
285 //
286 // return fANuMuTotXsc[index] if the index is in the array range
287 
289 {
290  if( index >= 0 && index < fIndex) return fANuMuTotXsc[index];
291  else
292  {
293  G4cout<<"Inproper index of fANuMuTotXsc array"<<G4endl;
294  return 0.;
295  }
296 }
297 
298 
300 //
301 // E_nu in GeV
302 
304 {
305  0.112103, 0.117359, 0.123119, 0.129443, 0.136404,
306  0.144084, 0.152576, 0.161991, 0.172458, 0.184126,
307  0.197171, 0.211801, 0.228261, 0.24684, 0.267887,
308  0.291816, 0.319125, 0.350417, 0.386422, 0.428032,
309  0.47634, 0.532692, 0.598756, 0.676612, 0.768868,
310  0.878812, 1.01062, 1.16963, 1.36271, 1.59876,
311  1.88943, 2.25002, 2.70086, 3.26916, 3.99166,
312  4.91843, 6.11836, 7.6872, 9.75942, 12.5259,
313  16.2605, 21.3615, 28.4141, 38.2903, 52.3062,
314  72.4763, 101.93, 145.6, 211.39, 312.172};
315 
317 //
318 // nu_mu CC xsc_tot/E_nu, in 10^-38 cm2/GeV
319 
321 {
322  0.0716001, 0.241013, 0.337793, 0.416033, 0.484616,
323  0.546945, 0.604661, 0.658623, 0.709277, 0.756815,
324  0.801263, 0.842519, 0.880396, 0.914642, 0.944968,
325  0.971069, 0.992655, 1.00947, 1.02133, 1.02812,
326  1.02987, 1.02671, 1.01892, 1.00689, 0.991167,
327  0.972618, 0.951518, 0.928806, 0.90521, 0.881404,
328  0.857978, 0.835424, 0.814112, 0.794314, 0.776204,
329  0.759884, 0.745394, 0.732719, 0.721809, 0.712164,
330  0.704299, 0.697804, 0.692491, 0.688137, 0.68448,
331  0.681232, 0.676128, 0.674154, 0.670553, 0.666034 };
332 
333 
334 
336 //
337 // anu_mu CC xsc_tot/E_nu, in 10^-38 cm2/GeV
338 
340 {
341  0.0291812, 0.0979725, 0.136884, 0.16794, 0.194698,
342  0.218468, 0.23992, 0.259241, 0.27665, 0.292251,
343  0.30612, 0.318314, 0.328886, 0.337885, 0.345464,
344  0.351495, 0.356131, 0.359448, 0.361531, 0.362474,
345  0.362382, 0.361365, 0.359538, 0.357024, 0.353943,
346  0.350422, 0.346685, 0.342662, 0.338567, 0.334514,
347  0.330612, 0.326966, 0.323668, 0.320805, 0.318451,
348  0.316671, 0.315514, 0.315013, 0.315187, 0.316036,
349  0.317541, 0.319667, 0.322362, 0.325556, 0.329159,
350  0.332577, 0.337133, 0.341214, 0.345128, 0.347657 };