| 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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4DynamicParticleFluctuation
33 //
34 // Author: V. Ivanchenko
35 //
36 // Creation date: 23.08.2024
37 //
38 // -------------------------------------------------------------------
39
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42
43 #include "G4DynamicParticleFluctuation.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "Randomize.hh"
47 #include "G4Poisson.hh"
48 #include "G4Material.hh"
49 #include "G4MaterialCutsCouple.hh"
50 #include "G4DynamicParticle.hh"
51 #include "G4Log.hh"
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54
55 G4DynamicParticleFluctuation::G4DynamicParticleFluctuation(const G4String& nam)
56 : G4UniversalFluctuation(nam)
57 {}
58
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60
61 void G4DynamicParticleFluctuation::InitialiseLocal(const G4DynamicParticle* part)
62 {
63 particleMass = part->GetMass();
64 const G4double q = part->GetCharge()/CLHEP::eplus;
65
66 // Derived quantities
67 m_Inv_particleMass = 1.0 / particleMass;
68 m_massrate = CLHEP::electron_mass_c2 * m_Inv_particleMass;
69 chargeSquare = q*q;
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74 G4double G4DynamicParticleFluctuation::SampleFluctuations(
75 const G4MaterialCutsCouple* couple,
76 const G4DynamicParticle* dp,
77 const G4double tcut,
78 const G4double tmax,
79 const G4double length,
80 const G4double averageLoss)
81 {
82 // Calculate actual loss from the mean loss.
83 // The model used to get the fluctuations is essentially the same
84 // as in Glandz in Geant3 (Cern program library W5013, phys332).
85 // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
86
87 // shortcut for very small loss or from a step nearly equal to the range
88 // (out of validity of the model)
89 //
90 if (averageLoss < minLoss) { return averageLoss; }
91 meanLoss = averageLoss;
92 const G4double tkin = dp->GetKineticEnergy();
93 //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
94
95 CLHEP::HepRandomEngine* rndmEngineF = G4Random::getTheEngine();
96
97 InitialiseLocal(dp);
98 const G4double gam = tkin * m_Inv_particleMass + 1.0;
99 const G4double gam2 = gam*gam;
100 const G4double beta = dp->GetBeta();
101 const G4double beta2 = beta*beta;
102
103 G4double loss(0.), siga(0.);
104
105 const G4Material* material = couple->GetMaterial();
106
107 // Gaussian regime
108 // for heavy particles only and conditions
109 // for Gauusian fluct. has been changed
110 //
111 if (particleMass > CLHEP::electron_mass_c2 &&
112 meanLoss >= minNumberInteractionsBohr*tcut && tmax <= 2.*tcut) {
113
114 siga = std::sqrt((tmax/beta2 - 0.5*tcut)*CLHEP::twopi_mc2_rcl2*
115 length*chargeSquare*material->GetElectronDensity());
116 const G4double sn = meanLoss/siga;
117
118 // thick target case
119 if (sn >= 2.0) {
120
121 const G4double twomeanLoss = meanLoss + meanLoss;
122 do {
123 loss = G4RandGauss::shoot(rndmEngineF, meanLoss, siga);
124 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
125 } while (0.0 > loss || twomeanLoss < loss);
126
127 // Gamma distribution
128 } else {
129
130 const G4double neff = sn*sn;
131 loss = meanLoss*G4RandGamma::shoot(rndmEngineF, neff, 1.0)/neff;
132 }
133 //G4cout << "Gauss: " << loss << G4endl;
134 return loss;
135 }
136
137 auto ioni = material->GetIonisation();
138 e0 = ioni->GetEnergy0fluct();
139
140 // very small step or low-density material
141 if(tcut <= e0) { return meanLoss; }
142
143 ipotFluct = ioni->GetMeanExcitationEnergy();
144 ipotLogFluct = ioni->GetLogMeanExcEnergy();
145
146 // width correction for small cuts
147 const G4double scaling = std::min(1.+0.5*CLHEP::keV/tcut, 1.50);
148 meanLoss /= scaling;
149
150 w2 = (tcut > ipotFluct) ?
151 G4Log(2.*CLHEP::electron_mass_c2*beta2*gam2) - beta2 : 0.0;
152 return SampleGlandz(rndmEngineF, material, tcut)*scaling;
153 }
154
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156
157
158 G4double G4DynamicParticleFluctuation::Dispersion(
159 const G4Material* material,
160 const G4DynamicParticle* dp,
161 const G4double tcut,
162 const G4double tmax,
163 const G4double length)
164 {
165 InitialiseLocal(dp);
166 const G4double beta = dp->GetBeta();
167 return (tmax/(beta*beta) - 0.5*tcut) * CLHEP::twopi_mc2_rcl2 * length
168 * material->GetElectronDensity() * chargeSquare;
169 }
170
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172