| 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: G4DynamicParticleIonisation
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 17.08.2024
37 //
38 // -------------------------------------------------------------------
39 //
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42
43 #include "G4DynamicParticleIonisation.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4DynamicParticleFluctuation.hh"
47 #include "G4EmSecondaryParticleType.hh"
48 #include "G4Electron.hh"
49 #include "G4EmParameters.hh"
50 #include "G4EmProcessSubType.hh"
51 #include "G4LossTableManager.hh"
52 #include "G4MaterialCutsCouple.hh"
53 #include "G4ProductionCutsTable.hh"
54 #include "G4Material.hh"
55 #include "G4Step.hh"
56 #include "G4Track.hh"
57 #include "G4Log.hh"
58
59 namespace
60 {
61 constexpr G4double ekinLimit = 0.2*CLHEP::MeV;
62 const G4double twoln10 = 2*G4Log(10.0);
63 }
64
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
67 G4DynamicParticleIonisation::G4DynamicParticleIonisation()
68 : G4VContinuousDiscreteProcess("dynPartIoni")
69 {
70 SetVerboseLevel(1);
71 SetProcessSubType(fDynamicIonisation);
72 theElectron = G4Electron::Electron();
73
74 lManager = G4LossTableManager::Instance();
75 lManager->Register(this);
76
77 fUrban = new G4DynamicParticleFluctuation();
78
79 // define these flags only once
80 auto param = G4EmParameters::Instance();
81 fFluct = param->LossFluctuation();
82 fLinLimit = 5*param->LinearLossLimit();
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87 G4DynamicParticleIonisation::~G4DynamicParticleIonisation()
88 {
89 lManager->DeRegister(this);
90 }
91
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93
94 void
95 G4DynamicParticleIonisation::BuildPhysicsTable(const G4ParticleDefinition&)
96 {
97 auto theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
98 fCuts = theCoupleTable->GetEnergyCutsVector(1);
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102
103 void G4DynamicParticleIonisation::PreStepInitialisation(const G4Track& track)
104 {
105 fCouple = track.GetMaterialCutsCouple();
106 fMaterial = fCouple->GetMaterial();
107 auto dpart = track.GetDynamicParticle();
108 fEkinPreStep = dpart->GetKineticEnergy();
109 fMass = std::max(dpart->GetMass(), CLHEP::electron_mass_c2);
110 fCharge = dpart->GetCharge()/CLHEP::eplus;
111 fRatio = fMass/CLHEP::proton_mass_c2;
112 fLowestEkin = ekinLimit*fRatio;
113 G4double tau = fEkinPreStep/fMass;
114 G4double ratio = CLHEP::electron_mass_c2/fMass;
115 fTmax = 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.) /
116 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
117 fCut = (*fCuts)[fCouple->GetIndex()];
118 fCut = std::max(fCut, fMaterial->GetIonisation()->GetMeanExcitationEnergy());
119 fCut = std::min(fCut, fTmax);
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123
124 G4double G4DynamicParticleIonisation::AlongStepGetPhysicalInteractionLength(
125 const G4Track&, G4double, G4double, G4double&,
126 G4GPILSelection* selection)
127 {
128 *selection = CandidateForSelection;
129
130 // no step limit for the time being
131 return DBL_MAX;
132 }
133
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135
136 G4double G4DynamicParticleIonisation::PostStepGetPhysicalInteractionLength(
137 const G4Track& track, G4double previousStepSize,
138 G4ForceCondition* condition)
139 {
140 *condition = NotForced;
141 G4double x = DBL_MAX;
142 G4double xsec = 0.0;
143
144 PreStepInitialisation(track);
145
146 if (fCharge != 0.0) {
147 xsec = ComputeCrossSection(fEkinPreStep);
148 }
149
150 if (xsec <= 0.0) {
151 theNumberOfInteractionLengthLeft = -1.0;
152 currentInteractionLength = DBL_MAX;
153
154 } else {
155 if (theNumberOfInteractionLengthLeft < 0.0) {
156
157 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
158 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
159
160 } else if(currentInteractionLength < DBL_MAX) {
161
162 // subtract NumberOfInteractionLengthLeft using previous step
163 theNumberOfInteractionLengthLeft -=
164 previousStepSize/currentInteractionLength;
165
166 theNumberOfInteractionLengthLeft =
167 std::max(theNumberOfInteractionLengthLeft, 0.0);
168 }
169 currentInteractionLength = 1.0/xsec;
170 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
171 }
172 #ifdef G4VERBOSE
173 if (verboseLevel>2) {
174 G4cout << "G4DynamicParticleIonisation::PostStepGetPhysicalInteractionLength ";
175 G4cout << " Process: " << GetProcessName()
176 << " for unknown particle Mass(GeV)=" << fMass/CLHEP::GeV
177 << " charge=" << fCharge
178 << " Material " << fMaterial->GetName()
179 << " Ekin(MeV)=" << fEkinPreStep/CLHEP::MeV
180 << " MFP(cm)=" << currentInteractionLength/CLHEP::cm
181 << " ProposedLength(cm)=" << x/CLHEP::cm <<G4endl;
182 }
183 #endif
184 return x;
185 }
186
187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188
189 G4VParticleChange*
190 G4DynamicParticleIonisation::AlongStepDoIt(const G4Track& track,
191 const G4Step& step)
192 {
193 fParticleChange.InitializeForAlongStep(track);
194
195 // no energy loss
196 if (fCharge == 0.0) { return &fParticleChange; }
197
198 // stop low-energy object
199 if (fEkinPreStep <= fLowestEkin) {
200 fParticleChange.SetProposedKineticEnergy(0.0);
201 fParticleChange.ProposeLocalEnergyDeposit(fEkinPreStep);
202 return &fParticleChange;
203 }
204
205 G4double length = step.GetStepLength();
206 G4double dedxPre = ComputeDEDX(fEkinPreStep);
207 G4double eloss = dedxPre*length;
208 G4double ekinPostStep = fEkinPreStep - eloss;
209
210 // correction for large step if it is not the last step
211 if (fEkinPreStep*fLinLimit < eloss && ekinPostStep > fLowestEkin) {
212 G4double dedxPost = ComputeDEDX(ekinPostStep);
213 eloss = (eloss + dedxPost*length)*0.5;
214 }
215
216 // do not sample fluctuations at the last step
217 if (fFluct && fEkinPreStep > eloss) {
218 eloss = fUrban->SampleFluctuations(fCouple, track.GetDynamicParticle(),
219 fCut, fTmax, length, eloss);
220 }
221
222 ekinPostStep = fEkinPreStep - eloss;
223
224 // stop low-energy object
225 if (ekinPostStep <= fLowestEkin) {
226 fParticleChange.SetProposedKineticEnergy(0.0);
227 fParticleChange.ProposeLocalEnergyDeposit(fEkinPreStep);
228 } else {
229 fParticleChange.SetProposedKineticEnergy(ekinPostStep);
230 fParticleChange.ProposeLocalEnergyDeposit(eloss);
231 }
232 return &fParticleChange;
233 }
234
235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236
237 G4VParticleChange*
238 G4DynamicParticleIonisation::PostStepDoIt(const G4Track& track, const G4Step&)
239 {
240 theNumberOfInteractionLengthLeft = -1.0;
241 fParticleChange.InitializeForPostStep(track);
242
243 auto dp = track.GetDynamicParticle();
244 G4double kinEnergy = dp->GetKineticEnergy();
245 const G4double totEnergy = kinEnergy + fMass;
246 const G4double beta2 = kinEnergy*(kinEnergy + 2.0*fMass)/(totEnergy*totEnergy);
247
248 G4double deltaKinEnergy, f;
249
250 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
251 G4double rndm[2];
252
253 // sampling without nuclear size effect
254 do {
255 rndmEngineMod->flatArray(2, rndm);
256 deltaKinEnergy = fCut*fTmax/(fCut*(1.0 - rndm[0]) + fTmax*rndm[0]);
257 f = 1.0 - beta2*deltaKinEnergy/fTmax;
258 // Loop checking, 14-Aug-2024, Vladimir Ivanchenko
259 } while( rndm[1] > f);
260
261 G4double deltaMomentum =
262 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*CLHEP::electron_mass_c2));
263 G4double cost = deltaKinEnergy * (totEnergy + CLHEP::electron_mass_c2) /
264 (deltaMomentum * dp->GetTotalMomentum());
265 cost = std::min(cost, 1.0);
266 const G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
267 const G4double phi = CLHEP::twopi*rndmEngineMod->flat();
268
269 G4ThreeVector deltaDirection(sint*std::cos(phi), sint*std::sin(phi), cost);
270 deltaDirection.rotateUz(dp->GetMomentumDirection());
271
272 // create G4DynamicParticle object for delta ray
273 auto delta = new G4DynamicParticle(theElectron, deltaDirection, deltaKinEnergy);
274 auto t = new G4Track(delta, track.GetGlobalTime(), track.GetPosition());
275 t->SetTouchableHandle(track.GetTouchableHandle());
276 t->SetCreatorModelID(fSecID);
277 fParticleChange.AddSecondary(t);
278
279 // Change kinematics of primary particle
280 kinEnergy -= deltaKinEnergy;
281 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
282 finalP = finalP.unit();
283
284 fParticleChange.SetProposedKineticEnergy(kinEnergy);
285 fParticleChange.SetProposedMomentumDirection(finalP);
286 return &fParticleChange;
287 }
288
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
290
291 G4double G4DynamicParticleIonisation::ComputeDEDX(G4double ekin)
292 {
293 G4double tau = ekin/fMass;
294 G4double gam = tau + 1.0;
295 G4double bg2 = tau * (tau + 2.0);
296 G4double beta2 = bg2/(gam*gam);
297 G4double xc = fCut/fTmax;
298
299 G4double exc = fMaterial->GetIonisation()->GetMeanExcitationEnergy();
300 G4double exc2 = exc*exc;
301
302 // general Bethe-Bloch formula
303 G4double dedx = G4Log(2.0*CLHEP::electron_mass_c2*bg2*fCut/exc2) - (1.0 + xc)*beta2;
304
305 // density correction
306 G4double x = G4Log(bg2)/twoln10;
307 dedx -= fMaterial->GetIonisation()->DensityCorrection(x);
308
309 // now compute the total ionization loss per volume
310 dedx *= CLHEP::twopi_mc2_rcl2*fCharge*fCharge*fMaterial->GetElectronDensity()/beta2;
311 dedx = std::max(dedx, 0.0);
312 return dedx;
313 }
314
315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
316
317 G4double G4DynamicParticleIonisation::ComputeCrossSection(G4double ekin)
318 {
319 G4double cross = 0.0;
320 if (fCut < fTmax) {
321
322 G4double totEnergy = ekin + fMass;
323 G4double energy2 = totEnergy*totEnergy;
324 G4double beta2 = ekin*(ekin + 2.0*fMass)/energy2;
325
326 cross = (fTmax - fCut)/(fCut*fTmax*beta2) - G4Log(fTmax/fCut)/fTmax;
327 cross *= CLHEP::twopi_mc2_rcl2*fCharge*fCharge*fMaterial->GetElectronDensity();
328 }
329 return cross;
330 }
331
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333
334 G4double G4DynamicParticleIonisation::GetMeanFreePath(const G4Track& /* track */, G4double,
335 G4ForceCondition* condition)
336 {
337 // Note: this method is not used at run-time, so its implementation is simplified.
338 // It might be eventually refined later.
339 *condition = NotForced;
340 return DBL_MAX;
341 }
342
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
344
345 G4double G4DynamicParticleIonisation::GetContinuousStepLimit(const G4Track&, G4double,
346 G4double, G4double&)
347 {
348 return DBL_MAX;
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352
353 void G4DynamicParticleIonisation::ProcessDescription(std::ostream& out) const
354 {
355 out << "G4DynamicParticleIonisation: dynamic ionisation" << G4endl;
356 }
357
358 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359