ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuBetheBlochModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuBetheBlochModel.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 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 //
32 // File name: G4MuBetheBlochModel
33 //
34 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
35 //
36 // Creation date: 09.08.2002
37 //
38 // Modifications:
39 //
40 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
42 // 27-01-03 Make models region aware (V.Ivanchenko)
43 // 13-02-03 Add name (V.Ivanchenko)
44 // 10-02-04 Calculation of radiative corrections using R.Kokoulin model (V.I)
45 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
46 // 12-04-05 Add usage of G4EmCorrections (V.Ivanchenko)
47 // 13-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
48 //
49 
50 //
51 // -------------------------------------------------------------------
52 //
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 #include "G4MuBetheBlochModel.hh"
58 #include "G4PhysicalConstants.hh"
59 #include "G4SystemOfUnits.hh"
60 #include "Randomize.hh"
61 #include "G4Electron.hh"
62 #include "G4LossTableManager.hh"
63 #include "G4EmCorrections.hh"
64 #include "G4ParticleChangeForLoss.hh"
65 #include "G4Log.hh"
66 #include "G4Exp.hh"
67 
68 G4double G4MuBetheBlochModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083, 0.5917,
69  0.7628, 0.8983, 0.9801 };
70 
71 G4double G4MuBetheBlochModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813, 0.1813,
72  0.1569, 0.1112, 0.0506 };
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 
76 using namespace std;
77 
78 G4MuBetheBlochModel::G4MuBetheBlochModel(const G4ParticleDefinition* p,
79  const G4String& nam)
80  : G4VEmModel(nam),
81  particle(nullptr),
82  limitKinEnergy(100.*keV),
83  logLimitKinEnergy(G4Log(limitKinEnergy)),
84  twoln10(2.0*G4Log(10.0)),
85  //bg2lim(0.0169),
86  //taulim(8.4146e-3),
87  alphaprime(fine_structure_const/twopi)
88 {
89  theElectron = G4Electron::Electron();
90  corr = G4LossTableManager::Instance()->EmCorrections();
91  fParticleChange = nullptr;
92 
93  // initial initialisation of memeber should be overwritten
94  // by SetParticle
95  mass = massSquare = ratio = 1.0;
96 
97  if(p) { SetParticle(p); }
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 
102 G4double G4MuBetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
103  const G4MaterialCutsCouple* couple)
104 {
105  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
109 
110 G4double G4MuBetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
111  G4double kinEnergy)
112 {
113  G4double tau = kinEnergy/mass;
114  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
115  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
116  return tmax;
117 }
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 
121 void G4MuBetheBlochModel::Initialise(const G4ParticleDefinition* p,
122  const G4DataVector&)
123 {
124  if(p) { SetParticle(p); }
125  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 
130 G4double G4MuBetheBlochModel::ComputeCrossSectionPerElectron(
131  const G4ParticleDefinition* p,
132  G4double kineticEnergy,
133  G4double cutEnergy,
134  G4double maxKinEnergy)
135 {
136  G4double cross = 0.0;
137  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
138  G4double maxEnergy = min(tmax,maxKinEnergy);
139  if(cutEnergy < maxEnergy) {
140 
141  G4double totEnergy = kineticEnergy + mass;
142  G4double energy2 = totEnergy*totEnergy;
143  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
144 
145  cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*G4Log(maxEnergy/cutEnergy)/tmax
146  + 0.5*(maxEnergy - cutEnergy)/energy2;
147 
148  // radiative corrections of R. Kokoulin
149  if (maxEnergy > limitKinEnergy) {
150 
151  G4double logtmax = G4Log(maxEnergy);
152  G4double logtmin = G4Log(max(cutEnergy,limitKinEnergy));
153  G4double logstep = logtmax - logtmin;
154  G4double dcross = 0.0;
155 
156  for (G4int ll=0; ll<8; ll++)
157  {
158  G4double ep = G4Exp(logtmin + xgi[ll]*logstep);
159  G4double a1 = G4Log(1.0 + 2.0*ep/electron_mass_c2);
160  G4double a3 = G4Log(4.0*totEnergy*(totEnergy - ep)/massSquare);
161  dcross += wgi[ll]*(1.0/ep - beta2/tmax + 0.5*ep/energy2)*a1*(a3 - a1);
162  }
163 
164  cross += dcross*logstep*alphaprime;
165  }
166 
167  cross *= twopi_mc2_rcl2/beta2;
168 
169  }
170 
171  // G4cout << "tmin= " << cutEnergy << " tmax= " << tmax
172  // << " cross= " << cross << G4endl;
173 
174  return cross;
175 }
176 
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178 
179 G4double G4MuBetheBlochModel::ComputeCrossSectionPerAtom(
180  const G4ParticleDefinition* p,
181  G4double kineticEnergy,
182  G4double Z, G4double,
183  G4double cutEnergy,
184  G4double maxEnergy)
185 {
186  G4double cross = Z*ComputeCrossSectionPerElectron
187  (p,kineticEnergy,cutEnergy,maxEnergy);
188  return cross;
189 }
190 
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
192 
193 G4double G4MuBetheBlochModel::CrossSectionPerVolume(
194  const G4Material* material,
195  const G4ParticleDefinition* p,
196  G4double kineticEnergy,
197  G4double cutEnergy,
198  G4double maxEnergy)
199 {
200  G4double eDensity = material->GetElectronDensity();
201  G4double cross = eDensity*ComputeCrossSectionPerElectron
202  (p,kineticEnergy,cutEnergy,maxEnergy);
203  return cross;
204 }
205 
206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
207 
208 G4double G4MuBetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
209  const G4ParticleDefinition* p,
210  G4double kineticEnergy,
211  G4double cut)
212 {
213  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
214  G4double tau = kineticEnergy/mass;
215  G4double cutEnergy = min(cut,tmax);
216  G4double gam = tau + 1.0;
217  G4double bg2 = tau * (tau+2.0);
218  G4double beta2 = bg2/(gam*gam);
219 
220  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
221  G4double eexc2 = eexc*eexc;
222 
223  G4double eDensity = material->GetElectronDensity();
224 
225  G4double dedx = G4Log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
226  -(1.0 + cutEnergy/tmax)*beta2;
227 
228  G4double totEnergy = kineticEnergy + mass;
229  G4double del = 0.5*cutEnergy/totEnergy;
230  dedx += del*del;
231 
232  // density correction
233  G4double x = G4Log(bg2)/twoln10;
234  dedx -= material->GetIonisation()->DensityCorrection(x);
235 
236  // shell correction
237  dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
238 
239  // now compute the total ionization loss
240 
241  if (dedx < 0.0) dedx = 0.0 ;
242 
243  // radiative corrections of R. Kokoulin
244  if (cutEnergy > limitKinEnergy) {
245 
246  G4double logtmax = G4Log(cutEnergy);
247  G4double logstep = logtmax - logLimitKinEnergy;
248  G4double dloss = 0.0;
249  G4double ftot2= 0.5/(totEnergy*totEnergy);
250 
251  for (G4int ll=0; ll<8; ll++)
252  {
253  G4double ep = G4Exp(logLimitKinEnergy + xgi[ll]*logstep);
254  G4double a1 = G4Log(1.0 + 2.0*ep/electron_mass_c2);
255  G4double a3 = G4Log(4.0*totEnergy*(totEnergy - ep)/massSquare);
256  dloss += wgi[ll]*(1.0 - beta2*ep/tmax + ep*ep*ftot2)*a1*(a3 - a1);
257  }
258  dedx += dloss*logstep*alphaprime;
259  }
260 
261  dedx *= twopi_mc2_rcl2*eDensity/beta2;
262 
263  //High order corrections
264  dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
265 
266  return dedx;
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 
271 void G4MuBetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
272  const G4MaterialCutsCouple*,
273  const G4DynamicParticle* dp,
274  G4double minKinEnergy,
275  G4double maxEnergy)
276 {
277  G4double tmax = MaxSecondaryKinEnergy(dp);
278  G4double maxKinEnergy = min(maxEnergy,tmax);
279  if(minKinEnergy >= maxKinEnergy) { return; }
280 
281  G4double kineticEnergy = dp->GetKineticEnergy();
282  G4double totEnergy = kineticEnergy + mass;
283  G4double etot2 = totEnergy*totEnergy;
284  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
285 
286  G4double grej = 1.;
287  if(tmax > limitKinEnergy) {
288  G4double a0 = G4Log(2.*totEnergy/mass);
289  grej += alphaprime*a0*a0;
290  }
291 
292  G4double deltaKinEnergy, f;
293 
294  // sampling follows ...
295  do {
296  G4double q = G4UniformRand();
297  deltaKinEnergy = minKinEnergy*maxKinEnergy
298  /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
299 
300 
301  f = 1.0 - beta2*deltaKinEnergy/tmax
302  + 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
303 
304  if(deltaKinEnergy > limitKinEnergy) {
305  G4double a1 = G4Log(1.0 + 2.0*deltaKinEnergy/electron_mass_c2);
306  G4double a3 = G4Log(4.0*totEnergy*(totEnergy - deltaKinEnergy)/massSquare);
307  f *= (1. + alphaprime*a1*(a3 - a1));
308  }
309 
310  if(f > grej) {
311  G4cout << "G4MuBetheBlochModel::SampleSecondary Warning! "
312  << "Majorant " << grej << " < "
313  << f << " for edelta= " << deltaKinEnergy
314  << " tmin= " << minKinEnergy << " max= " << maxKinEnergy
315  << G4endl;
316  }
317  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
318  } while( grej*G4UniformRand() > f );
319 
320  G4double deltaMomentum =
321  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
322  G4double totalMomentum = totEnergy*sqrt(beta2);
323  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
324  (deltaMomentum * totalMomentum);
325 
326  G4double sint = sqrt(1.0 - cost*cost);
327 
328  G4double phi = twopi * G4UniformRand() ;
329 
330  G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
331  G4ThreeVector direction = dp->GetMomentumDirection();
332  deltaDirection.rotateUz(direction);
333 
334  // primary change
335  kineticEnergy -= deltaKinEnergy;
336  G4ThreeVector dir = totalMomentum*direction - deltaMomentum*deltaDirection;
337  direction = dir.unit();
338  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
339  fParticleChange->SetProposedMomentumDirection(direction);
340 
341  // create G4DynamicParticle object for delta ray
342  G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
343  deltaDirection,deltaKinEnergy);
344  vdp->push_back(delta);
345 }
346 
347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......