| 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 #include "G4ContinuousGainOfEnergy.hh"
29
30 #include "G4EmCorrections.hh"
31 #include "G4LossTableManager.hh"
32 #include "G4Material.hh"
33 #include "G4MaterialCutsCouple.hh"
34 #include "G4ParticleChange.hh"
35 #include "G4ParticleDefinition.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4Step.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4VEmFluctuationModel.hh"
40 #include "G4VEmModel.hh"
41 #include "G4VEnergyLossProcess.hh"
42 #include "G4VParticleChange.hh"
43
44 ///////////////////////////////////////////////////////
45 G4ContinuousGainOfEnergy::G4ContinuousGainOfEnergy(const G4String& name,
46 G4ProcessType type)
47 : G4VContinuousProcess(name, type)
48 {}
49
50 ///////////////////////////////////////////////////////
51 G4ContinuousGainOfEnergy::~G4ContinuousGainOfEnergy() {}
52
53 ///////////////////////////////////////////////////////
54 void G4ContinuousGainOfEnergy::ProcessDescription(std::ostream& out) const
55 {
56 out << "Continuous process acting on adjoint particles to compute the "
57 "continuous gain of energy of charged particles when they are "
58 "tracked back.\n";
59 }
60
61 ///////////////////////////////////////////////////////
62 void G4ContinuousGainOfEnergy::SetDirectParticle(G4ParticleDefinition* p)
63 {
64 fDirectPartDef = p;
65 if(fDirectPartDef->GetParticleType() == "nucleus")
66 {
67 fIsIon = true;
68 fMassRatio = proton_mass_c2 / fDirectPartDef->GetPDGMass();
69 }
70 }
71
72 ///////////////////////////////////////////////////////
73 G4VParticleChange* G4ContinuousGainOfEnergy::AlongStepDoIt(const G4Track& track,
74 const G4Step& step)
75 {
76 // Caution in this method the step length should be the true step length
77 // A problem is that this is computed by the multiple scattering that does
78 // not know the energy at the end of the adjoint step. This energy is used
79 // during the forward sim. Nothing we can really do against that at this
80 // time. This is inherent to the MS method
81
82 aParticleChange.Initialize(track);
83
84 // Get the actual (true) Step length
85 G4double length = step.GetStepLength();
86 G4double degain = 0.0;
87
88 // Compute this for weight change after continuous energy loss
89 G4double DEDX_before =
90 fDirectEnergyLossProcess->GetDEDX(fPreStepKinEnergy, fCurrentCouple);
91
92 // For the fluctuation we generate a new dynamic particle with energy
93 // = preEnergy+egain and then compute the fluctuation given in the direct
94 // case.
95 G4DynamicParticle* dynParticle = new G4DynamicParticle();
96 *dynParticle = *(track.GetDynamicParticle());
97 dynParticle->SetDefinition(fDirectPartDef);
98 G4double Tkin = dynParticle->GetKineticEnergy();
99
100 G4double dlength = length;
101 if(Tkin != fPreStepKinEnergy && fIsIon)
102 {
103 G4double chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
104 fDirectPartDef, fCurrentMaterial, Tkin);
105 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio, chargeSqRatio);
106 }
107
108 G4double r = fDirectEnergyLossProcess->GetRange(Tkin, fCurrentCouple);
109 if(dlength <= fLinLossLimit * r)
110 {
111 degain = DEDX_before * dlength;
112 }
113 else
114 {
115 G4double x = r + dlength;
116 G4double E = fDirectEnergyLossProcess->GetKineticEnergy(x, fCurrentCouple);
117 if(fIsIon)
118 {
119 G4double chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
120 fDirectPartDef, fCurrentMaterial, E);
121 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio, chargeSqRatio);
122 G4double x1 = fDirectEnergyLossProcess->GetRange(E, fCurrentCouple);
123
124 G4int ii = 0;
125 constexpr G4int iimax = 100;
126 while(std::abs(x - x1) > 0.01 * x)
127 {
128 E = fDirectEnergyLossProcess->GetKineticEnergy(x, fCurrentCouple);
129 chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
130 fDirectPartDef, fCurrentMaterial, E);
131 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio,
132 chargeSqRatio);
133 x1 = fDirectEnergyLossProcess->GetRange(E, fCurrentCouple);
134 ++ii;
135 if(ii >= iimax)
136 {
137 break;
138 }
139 }
140 }
141
142 degain = E - Tkin;
143 }
144 G4double tmax = fCurrentModel->MaxSecondaryKinEnergy(dynParticle);
145 fCurrentTcut = std::min(fCurrentTcut, tmax);
146
147 dynParticle->SetKineticEnergy(Tkin + degain);
148
149 // Corrections, which cannot be tabulated for ions
150 fCurrentModel->CorrectionsAlongStep(fCurrentCouple, dynParticle, dlength, degain);
151
152 // Sample fluctuations
153 G4double deltaE = 0.;
154 if(fLossFluctuationFlag)
155 {
156 deltaE = fCurrentModel->GetModelOfFluctuations()->SampleFluctuations(
157 fCurrentCouple, dynParticle, fCurrentTcut, tmax, dlength, degain)
158 - degain;
159 }
160
161 G4double egain = degain + deltaE;
162 if(egain <= 0.)
163 egain = degain;
164 Tkin += egain;
165 dynParticle->SetKineticEnergy(Tkin);
166
167 delete dynParticle;
168
169 if(fIsIon)
170 {
171 G4double chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
172 fDirectPartDef, fCurrentMaterial, Tkin);
173 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio, chargeSqRatio);
174 }
175
176 G4double DEDX_after = fDirectEnergyLossProcess->GetDEDX(Tkin, fCurrentCouple);
177 G4double weight_correction = DEDX_after / DEDX_before;
178
179 aParticleChange.ProposeEnergy(Tkin);
180
181 // Caution!!! It is important to select the weight of the post_step_point
182 // as the current weight and not the weight of the track, as the weight of
183 // the track is changed after having applied all the along_step_do_it.
184
185 G4double new_weight =
186 weight_correction * step.GetPostStepPoint()->GetWeight();
187 aParticleChange.SetParentWeightByProcess(false);
188 aParticleChange.ProposeParentWeight(new_weight);
189
190 return &aParticleChange;
191 }
192
193 ///////////////////////////////////////////////////////
194 void G4ContinuousGainOfEnergy::SetLossFluctuations(G4bool val)
195 {
196 if(val && !fLossFluctuationArePossible)
197 return;
198 fLossFluctuationFlag = val;
199 }
200
201 ///////////////////////////////////////////////////////
202 G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
203 G4double, G4double,
204 G4double&)
205 {
206 DefineMaterial(track.GetMaterialCutsCouple());
207
208 fPreStepKinEnergy = track.GetKineticEnergy();
209 fCurrentModel = fDirectEnergyLossProcess->SelectModelForMaterial(
210 track.GetKineticEnergy() * fMassRatio, fCurrentCoupleIndex);
211 G4double emax_model = fCurrentModel->HighEnergyLimit();
212 G4double preStepChargeSqRatio = 0.;
213 if(fIsIon)
214 {
215 G4double chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
216 fDirectPartDef, fCurrentMaterial, fPreStepKinEnergy);
217 preStepChargeSqRatio = chargeSqRatio;
218 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio,
219 preStepChargeSqRatio);
220 }
221
222 G4double maxE = 1.1 * fPreStepKinEnergy;
223
224 if(fPreStepKinEnergy < fCurrentTcut)
225 maxE = std::min(fCurrentTcut, maxE);
226
227 maxE = std::min(emax_model * 1.001, maxE);
228
229 G4double preStepRange =
230 fDirectEnergyLossProcess->GetRange(fPreStepKinEnergy, fCurrentCouple);
231
232 if(fIsIon)
233 {
234 G4double chargeSqRatioAtEmax = fCurrentModel->GetChargeSquareRatio(
235 fDirectPartDef, fCurrentMaterial, maxE);
236 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio,
237 chargeSqRatioAtEmax);
238 }
239
240 G4double r1 = fDirectEnergyLossProcess->GetRange(maxE, fCurrentCouple);
241
242 if(fIsIon)
243 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio,
244 preStepChargeSqRatio);
245
246 return std::max(r1 - preStepRange, 0.001 * mm);
247 }
248
249 ///////////////////////////////////////////////////////
250 void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track&,
251 G4double energy)
252 {
253 G4double ChargeSqRatio =
254 G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(
255 fDirectPartDef, fCurrentMaterial, energy);
256 if(fDirectEnergyLossProcess)
257 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio, ChargeSqRatio);
258 }
259