ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuBremsstrahlungModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuBremsstrahlungModel.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 file
30 //
31 //
32 // File name: G4MuBremsstrahlungModel
33 //
34 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
35 //
36 // Creation date: 24.06.2002
37 //
38 // Modifications:
39 //
40 // 04-12-02 Change G4DynamicParticle constructor in PostStepDoIt (V.Ivanchenko)
41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
42 // 24-01-03 Fix for compounds (V.Ivanchenko)
43 // 27-01-03 Make models region aware (V.Ivanchenko)
44 // 13-02-03 Add name (V.Ivanchenko)
45 // 10-02-04 Add lowestKinEnergy (V.Ivanchenko)
46 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
47 // 03-08-05 Angular correlations according to PRM (V.Ivanchenko)
48 // 13-02-06 add ComputeCrossSectionPerAtom (mma)
49 // 21-03-06 Fix problem of initialisation in case when cuts are not defined (VI)
50 // 07-11-07 Improve sampling of final state (A.Bogdanov)
51 // 28-02-08 Use precomputed Z^1/3 and Log(A) (V.Ivanchenko)
52 // 31-05-13 Use element selectors instead of local data structure (V.Ivanchenko)
53 //
54 
55 //
56 // Class Description:
57 //
58 //
59 // -------------------------------------------------------------------
60 //
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 
64 #include "G4MuBremsstrahlungModel.hh"
65 #include "G4PhysicalConstants.hh"
66 #include "G4SystemOfUnits.hh"
67 #include "G4Gamma.hh"
68 #include "G4MuonMinus.hh"
69 #include "G4MuonPlus.hh"
70 #include "Randomize.hh"
71 #include "G4Material.hh"
72 #include "G4Element.hh"
73 #include "G4ElementVector.hh"
74 #include "G4ProductionCutsTable.hh"
75 #include "G4ParticleChangeForLoss.hh"
76 #include "G4Log.hh"
77 #include "G4Exp.hh"
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81 
82 using namespace std;
83 
84 const G4double G4MuBremsstrahlungModel::xgi[] =
85  {0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
86 const G4double G4MuBremsstrahlungModel::wgi[] =
87  {0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
88 G4double G4MuBremsstrahlungModel::fDN[] = {0.0};
89 
90 G4MuBremsstrahlungModel::G4MuBremsstrahlungModel(const G4ParticleDefinition* p,
91  const G4String& nam)
92  : G4VEmModel(nam),
93  particle(nullptr),
94  sqrte(sqrt(G4Exp(1.))),
95  bh(202.4),
96  bh1(446.),
97  btf(183.),
98  btf1(1429.),
99  fParticleChange(nullptr),
100  lowestKinEnergy(1.0*GeV),
101  minThreshold(0.9*keV)
102 {
103  theGamma = G4Gamma::Gamma();
104  nist = G4NistManager::Instance();
105 
106  lowestKinEnergy = 1.*GeV;
107 
108  mass = rmass = cc = coeff = 1.0;
109 
110  if(0.0 == fDN[1]) {
111  for(G4int i=1; i<93; ++i) {
112  G4double dn = 1.54*nist->GetA27(i);
113  fDN[i] = dn;
114  if(1 < i) {
115  fDN[i] /= std::pow(dn, 1./G4double(i));
116  }
117  }
118  }
119 
120  if(p) { SetParticle(p); }
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 
125 G4double G4MuBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
126  const G4MaterialCutsCouple*)
127 {
128  return minThreshold;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 G4double G4MuBremsstrahlungModel::MinPrimaryEnergy(const G4Material*,
134  const G4ParticleDefinition*,
135  G4double cut)
136 {
137  return std::max(lowestKinEnergy,cut);
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 
142 void G4MuBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
143  const G4DataVector& cuts)
144 {
145  if(p) { SetParticle(p); }
146 
147  // define pointer to G4ParticleChange
148  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
149 
150  if(IsMaster() && p == particle && lowestKinEnergy < HighEnergyLimit()) {
151  InitialiseElementSelectors(p, cuts);
152  }
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156 
157 void G4MuBremsstrahlungModel::InitialiseLocal(const G4ParticleDefinition* p,
158  G4VEmModel* masterModel)
159 {
160  if(p == particle && lowestKinEnergy < HighEnergyLimit()) {
161  SetElementSelectors(masterModel->GetElementSelectors());
162  }
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166 
167 G4double G4MuBremsstrahlungModel::ComputeDEDXPerVolume(
168  const G4Material* material,
169  const G4ParticleDefinition*,
170  G4double kineticEnergy,
171  G4double cutEnergy)
172 {
173  G4double dedx = 0.0;
174  if (kineticEnergy <= lowestKinEnergy) { return dedx; }
175 
176  G4double tmax = kineticEnergy;
177  G4double cut = std::min(cutEnergy,tmax);
178  if(cut < minThreshold) { cut = minThreshold; }
179 
180  const G4ElementVector* theElementVector = material->GetElementVector();
181  const G4double* theAtomicNumDensityVector =
182  material->GetAtomicNumDensityVector();
183 
184  // loop for elements in the material
185  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
186 
187  G4double loss =
188  ComputMuBremLoss((*theElementVector)[i]->GetZ(), kineticEnergy, cut);
189 
190  dedx += loss*theAtomicNumDensityVector[i];
191  }
192  // G4cout << "BR e= " << kineticEnergy << " dedx= " << dedx << G4endl;
193  if(dedx < 0.) dedx = 0.;
194  return dedx;
195 }
196 
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
198 
199 G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z,
200  G4double tkin, G4double cut)
201 {
202  G4double totalEnergy = mass + tkin;
203  static const G4double ak1 = 0.05;
204  static const G4int k2=5;
205  G4double loss = 0.;
206 
207  G4double vcut = cut/totalEnergy;
208  G4double vmax = tkin/totalEnergy;
209 
210  G4double aaa = 0.;
211  G4double bbb = vcut;
212  if(vcut>vmax) { bbb = vmax; }
213  G4int kkk = (G4int)((bbb-aaa)/ak1)+k2;
214  if(kkk < 1) { kkk = 1; }
215 
216  G4double hhh=(bbb-aaa)/G4double(kkk);
217 
218  G4double aa = aaa;
219  for(G4int l=0; l<kkk; l++)
220  {
221  for(G4int i=0; i<6; i++)
222  {
223  G4double ep = (aa + xgi[i]*hhh)*totalEnergy;
224  loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
225  }
226  aa += hhh;
227  }
228 
229  loss *=hhh*totalEnergy ;
230 
231  return loss;
232 }
233 
234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 
236 G4double G4MuBremsstrahlungModel::ComputeMicroscopicCrossSection(
237  G4double tkin,
238  G4double Z,
239  G4double cut)
240 {
241  G4double totalEnergy = tkin + mass;
242  static const G4double ak1 = 2.3;
243  static const G4int k2 = 4;
244  G4double cross = 0.;
245 
246  if(cut >= tkin) return cross;
247 
248  G4double vcut = cut/totalEnergy;
249  G4double vmax = tkin/totalEnergy;
250 
251  G4double aaa = G4Log(vcut);
252  G4double bbb = G4Log(vmax);
253  G4int kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
254  if(kkk < 1) { kkk = 1; }
255 
256  G4double hhh = (bbb-aaa)/G4double(kkk);
257 
258  G4double aa = aaa;
259 
260  for(G4int l=0; l<kkk; l++)
261  {
262  for(G4int i=0; i<6; i++)
263  {
264  G4double ep = G4Exp(aa + xgi[i]*hhh)*totalEnergy;
265  cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
266  }
267  aa += hhh;
268  }
269 
270  cross *=hhh;
271 
272  //G4cout << "BR e= " << tkin<< " cross= " << cross/barn << G4endl;
273 
274  return cross;
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278 
279 G4double G4MuBremsstrahlungModel::ComputeDMicroscopicCrossSection(
280  G4double tkin,
281  G4double Z,
282  G4double gammaEnergy)
283 // differential cross section
284 {
285  G4double dxsection = 0.;
286 
287  if(gammaEnergy > tkin) { return dxsection; }
288 
289  G4double E = tkin + mass ;
290  G4double v = gammaEnergy/E ;
291  G4double delta = 0.5*mass*mass*v/(E-gammaEnergy) ;
292  G4double rab0 = delta*sqrte ;
293 
294  G4int iz = G4lrint(Z);
295  if(iz < 1) { iz = 1; }
296  else if(iz > 92) { iz = 92; }
297 
298  G4double z13 = 1.0/nist->GetZ13(iz);
299  G4double dnstar = fDN[iz];
300 
301  G4double b,b1;
302 
303  if(1 == iz) {
304  b = bh;
305  b1 = bh1;
306  } else {
307  b = btf;
308  b1 = btf1;
309  }
310 
311  // nucleus contribution logarithm
312  G4double rab1=b*z13;
313  G4double fn=G4Log(rab1/(dnstar*(electron_mass_c2+rab0*rab1))*
314  (mass+delta*(dnstar*sqrte-2.))) ;
315  if(fn <0.) { fn = 0.; }
316  // electron contribution logarithm
317  G4double epmax1=E/(1.+0.5*mass*rmass/E) ;
318  G4double fe=0.;
319  if(gammaEnergy<epmax1)
320  {
321  G4double rab2=b1*z13*z13 ;
322  fe=G4Log(rab2*mass/((1.+delta*rmass/(electron_mass_c2*sqrte))*
323  (electron_mass_c2+rab0*rab2))) ;
324  if(fe<0.) { fe=0.; }
325  }
326 
327  dxsection = coeff*(1.-v*(1. - 0.75*v))*Z*(fn*Z + fe)/gammaEnergy;
328 
329  return dxsection;
330 }
331 
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
333 
334 G4double G4MuBremsstrahlungModel::ComputeCrossSectionPerAtom(
335  const G4ParticleDefinition*,
336  G4double kineticEnergy,
337  G4double Z, G4double,
338  G4double cutEnergy,
339  G4double maxEnergy)
340 {
341  G4double cross = 0.0;
342  if (kineticEnergy <= lowestKinEnergy) return cross;
343  G4double tmax = std::min(maxEnergy, kineticEnergy);
344  G4double cut = std::min(cutEnergy, kineticEnergy);
345  if(cut < minThreshold) cut = minThreshold;
346  if (cut >= tmax) return cross;
347 
348  cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
349  if(tmax < kineticEnergy) {
350  cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
351  }
352  return cross;
353 }
354 
355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
356 
357 void G4MuBremsstrahlungModel::SampleSecondaries(
358  std::vector<G4DynamicParticle*>* vdp,
359  const G4MaterialCutsCouple* couple,
360  const G4DynamicParticle* dp,
361  G4double minEnergy,
362  G4double maxEnergy)
363 {
364  G4double kineticEnergy = dp->GetKineticEnergy();
365  // check against insufficient energy
366  G4double tmax = std::min(kineticEnergy, maxEnergy);
367  G4double tmin = std::min(kineticEnergy, minEnergy);
368  if(tmin < minThreshold) tmin = minThreshold;
369  if(tmin >= tmax) return;
370 
371  // ===== sampling of energy transfer ======
372 
373  G4ParticleMomentum partDirection = dp->GetMomentumDirection();
374 
375  // select randomly one element constituing the material
376  const G4Element* anElement = SelectRandomAtom(couple,particle,kineticEnergy);
377  G4double Z = anElement->GetZ();
378 
379  G4double totalEnergy = kineticEnergy + mass;
380  G4double totalMomentum = sqrt(kineticEnergy*(kineticEnergy + 2.0*mass));
381 
382  G4double func1 = tmin*
383  ComputeDMicroscopicCrossSection(kineticEnergy,Z,tmin);
384 
385  G4double lnepksi, epksi;
386  G4double func2;
387 
388  G4double xmin = G4Log(tmin/MeV);
389  G4double xmax = G4Log(kineticEnergy/tmin);
390 
391  do {
392  lnepksi = xmin + G4UniformRand()*xmax;
393  epksi = MeV*G4Exp(lnepksi);
394  func2 = epksi*ComputeDMicroscopicCrossSection(kineticEnergy,Z,epksi);
395 
396  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
397  } while(func2 < func1*G4UniformRand());
398 
399  G4double gEnergy = epksi;
400 
401  // ===== sample angle =====
402 
403  G4double gam = totalEnergy/mass;
404  G4double rmax = gam*std::min(1.0, totalEnergy/gEnergy - 1.0);
405  G4double rmax2= rmax*rmax;
406  G4double x = G4UniformRand()*rmax2/(1.0 + rmax2);
407 
408  G4double theta = sqrt(x/(1.0 - x))/gam;
409  G4double sint = sin(theta);
410  G4double phi = twopi * G4UniformRand() ;
411  G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
412 
413  G4ThreeVector gDirection(dirx, diry, dirz);
414  gDirection.rotateUz(partDirection);
415 
416  partDirection *= totalMomentum;
417  partDirection -= gEnergy*gDirection;
418  partDirection = partDirection.unit();
419 
420  // primary change
421  kineticEnergy -= gEnergy;
422  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
423  fParticleChange->SetProposedMomentumDirection(partDirection);
424 
425  // save secondary
426  G4DynamicParticle* aGamma =
427  new G4DynamicParticle(theGamma,gDirection,gEnergy);
428  vdp->push_back(aGamma);
429 }
430 
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......