ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KokoulinMuonNuclearXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4KokoulinMuonNuclearXS.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 // Author: D.H. Wright
28 // Date: 1 February 2011
29 //
30 // Modified:
31 //
32 // 19 Aug 2011, V.Ivanchenko move to new design and make x-section per element
33 
34 // Description: use Kokoulin's parameterized calculation of virtual
35 // photon production cross section and conversion to
36 // real photons.
37 
39 
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4PhysicsLogVector.hh"
43 #include "G4PhysicsVector.hh"
44 #include "G4MuonMinus.hh"
45 #include "G4MuonPlus.hh"
46 #include "G4NucleiProperties.hh"
47 #include "G4NistManager.hh"
48 #include "G4Log.hh"
49 #include "G4Exp.hh"
50 
51 // factory
52 #include "G4CrossSectionFactory.hh"
53 //
55 
57 
59  :G4VCrossSectionDataSet(Default_Name()),
60  LowestKineticEnergy(1*GeV), HighestKineticEnergy(1*PeV),
61  TotBin(60), CutFixed(0.2*GeV), isInitialized(false), isMaster(false)
62 {}
63 
65 {
66  if (isMaster) {
67  for(G4int i=0; i<MAXZMUN; ++i) {
68  delete theCrossSection[i];
69  theCrossSection[i] = 0;
70  }
71  }
72 }
73 
74 
75 void
77 {
78  outFile << "G4KokoulinMuonNuclearXS provides the total inelastic\n"
79  << "cross section for mu- and mu+ interactions with nuclei.\n"
80  << "R. Kokoulin's approximation of the Borog and Petrukhin double\n"
81  << "differential cross section at high energy and low Q**2 is integrated\n"
82  << "over the muon energy loss to get the total cross section as a\n"
83  << "function of muon kinetic energy\n" ;
84 }
85 
86 
87 G4bool
89  G4int, const G4Material*)
90 {
91  return true;
92 }
93 
94 void
96 {
97  if(!isInitialized) {
98  isInitialized = true;
99  for(G4int i=0; i<MAXZMUN; ++i) {
100  if(theCrossSection[i]) { return; }
101  }
102  isMaster = true;
103  }
105 }
106 
108 {
109  G4double energy, A, Value;
110  G4int Z;
111 
113  const G4ElementTable* theElementTable = G4Element::GetElementTable();
114  G4NistManager* nistManager = G4NistManager::Instance();
115 
116  for (G4int j = 0; j < nEl; j++) {
117  Z = G4lrint((*theElementTable)[j]->GetZ());
118 
119  //AR-24Apr2018 Switch to treat transuranic elements as uranium
120  const G4bool isHeavyElementAllowed = true; if ( isHeavyElementAllowed && Z>92 ) Z=92;
121 
122  A = nistManager->GetAtomicMassAmu(Z);
123  if(Z < MAXZMUN && !theCrossSection[Z]) {
126  TotBin);
127  for (G4int i = 0; i <= TotBin; ++i) {
128  energy = theCrossSection[Z]->Energy(i);
129  Value = ComputeMicroscopicCrossSection(energy, A);
130  theCrossSection[Z]->PutValue(i,Value);
131  }
132  }
133  }
134 }
135 
138 {
139  // Calculate cross section (differential in muon incident kinetic energy) by
140  // integrating the double differential cross section over the energy loss
141 
142  static const G4double xgi[] =
143  {0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801};
144  static const G4double wgi[] =
145  {0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506};
146  static const G4double ak1 = 6.9;
147  static const G4double ak2 = 1.0;
148 
150 
151  G4double CrossSection = 0.0;
152  if (KineticEnergy <= CutFixed) return CrossSection;
153 
154  G4double epmin = CutFixed;
155  G4double epmax = KineticEnergy + Mass - 0.5*proton_mass_c2;
156  if (epmax <= epmin) return CrossSection; // NaN bug correction
157 
158  G4double aaa = G4Log(epmin);
159  G4double bbb = G4Log(epmax);
160  G4int kkk = std::max(1,G4int((bbb-aaa)/ak1 +ak2));
161  G4double hhh = (bbb-aaa)/kkk ;
162  G4double epln;
163  G4double ep;
164  G4double x;
165 
166  for (G4int l = 0; l < kkk; ++l) {
167  x = aaa + hhh*l;
168  for (G4int ll = 0; ll < 8; ++ll) {
169  epln=x+xgi[ll]*hhh;
170  ep = G4Exp(epln);
171  CrossSection +=
172  ep*wgi[ll]*ComputeDDMicroscopicCrossSection(KineticEnergy, 0, A, ep);
173  }
174  }
175 
176  CrossSection *= hhh ;
177  if (CrossSection < 0.) { CrossSection = 0.; }
178  return CrossSection;
179 }
180 
184 {
185  // Calculate the double-differential microscopic cross section (in muon
186  // incident kinetic energy and energy loss) using the cross section formula
187  // of R.P. Kokoulin (18/01/98)
188 
189  static const G4double alam2 = 0.400*GeV*GeV;
190  static const G4double alam = 0.632456*GeV;
191  static const G4double coeffn = fine_structure_const/pi;
192 
193  G4double ParticleMass = G4MuonMinus::MuonMinus()->GetPDGMass();
194  G4double TotalEnergy = KineticEnergy + ParticleMass;
195 
196  G4double DCrossSection = 0.;
197 
198  if ((epsilon >= TotalEnergy - 0.5*proton_mass_c2) ||
199  (epsilon <= CutFixed) ) { return DCrossSection; }
200 
201  G4double ep = epsilon/GeV;
202  G4double aeff = 0.22*A+0.78*G4Exp(0.89*G4Log(A)); //shadowing
203  G4double sigph = (49.2+11.1*G4Log(ep)+151.8/std::sqrt(ep))*microbarn;
204 
205  G4double v = epsilon/TotalEnergy;
206  G4double v1 = 1.-v;
207  G4double v2 = v*v;
208  G4double mass2 = ParticleMass*ParticleMass;
209 
210  G4double up = TotalEnergy*TotalEnergy*v1/mass2*(1.+mass2*v2/(alam2*v1));
211  G4double down = 1.+epsilon/alam*(1.+alam/(2.*proton_mass_c2)+epsilon/alam);
212 
213  DCrossSection = coeffn*aeff*sigph/epsilon*
214  (-v1+(v1+0.5*v2*(1.+2.*mass2/alam2))*G4Log(up/down));
215 
216  if (DCrossSection < 0.) { DCrossSection = 0.; }
217  return DCrossSection;
218 }
219 
222  G4int Z, const G4Material*)
223 {
224  //AR-24Apr2018 Switch to treat transuranic elements as uranium
225  const G4bool isHeavyElementAllowed = true; if ( isHeavyElementAllowed && Z>92 ) Z=92;
226 
227  return theCrossSection[Z]->Value(aPart->GetKineticEnergy());
228 }
229