| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 27
28 #include "G4ContinuousGainOfEnergy.hh" 28 #include "G4ContinuousGainOfEnergy.hh"
29 29
30 #include "G4EmCorrections.hh" 30 #include "G4EmCorrections.hh"
31 #include "G4LossTableManager.hh" 31 #include "G4LossTableManager.hh"
32 #include "G4Material.hh" 32 #include "G4Material.hh"
33 #include "G4MaterialCutsCouple.hh" 33 #include "G4MaterialCutsCouple.hh"
34 #include "G4ParticleChange.hh" 34 #include "G4ParticleChange.hh"
35 #include "G4ParticleDefinition.hh" 35 #include "G4ParticleDefinition.hh"
36 #include "G4PhysicalConstants.hh" 36 #include "G4PhysicalConstants.hh"
37 #include "G4Step.hh" 37 #include "G4Step.hh"
38 #include "G4SystemOfUnits.hh" 38 #include "G4SystemOfUnits.hh"
39 #include "G4VEmFluctuationModel.hh" 39 #include "G4VEmFluctuationModel.hh"
40 #include "G4VEmModel.hh" 40 #include "G4VEmModel.hh"
41 #include "G4VEnergyLossProcess.hh" 41 #include "G4VEnergyLossProcess.hh"
42 #include "G4VParticleChange.hh" 42 #include "G4VParticleChange.hh"
43 43
44 ////////////////////////////////////////////// 44 ///////////////////////////////////////////////////////
45 G4ContinuousGainOfEnergy::G4ContinuousGainOfEn 45 G4ContinuousGainOfEnergy::G4ContinuousGainOfEnergy(const G4String& name,
46 46 G4ProcessType type)
47 : G4VContinuousProcess(name, type) 47 : G4VContinuousProcess(name, type)
48 {} 48 {}
49 49
50 ////////////////////////////////////////////// 50 ///////////////////////////////////////////////////////
51 G4ContinuousGainOfEnergy::~G4ContinuousGainOfE 51 G4ContinuousGainOfEnergy::~G4ContinuousGainOfEnergy() {}
52 52
53 ////////////////////////////////////////////// 53 ///////////////////////////////////////////////////////
54 void G4ContinuousGainOfEnergy::ProcessDescript 54 void G4ContinuousGainOfEnergy::ProcessDescription(std::ostream& out) const
55 { 55 {
56 out << "Continuous process acting on adjoint 56 out << "Continuous process acting on adjoint particles to compute the "
57 "continuous gain of energy of charged 57 "continuous gain of energy of charged particles when they are "
58 "tracked back.\n"; 58 "tracked back.\n";
59 } 59 }
60 60
61 ////////////////////////////////////////////// 61 ///////////////////////////////////////////////////////
62 void G4ContinuousGainOfEnergy::SetDirectPartic 62 void G4ContinuousGainOfEnergy::SetDirectParticle(G4ParticleDefinition* p)
63 { 63 {
64 fDirectPartDef = p; 64 fDirectPartDef = p;
65 if(fDirectPartDef->GetParticleType() == "nuc 65 if(fDirectPartDef->GetParticleType() == "nucleus")
66 { 66 {
67 fIsIon = true; 67 fIsIon = true;
68 fMassRatio = proton_mass_c2 / fDirectPartD 68 fMassRatio = proton_mass_c2 / fDirectPartDef->GetPDGMass();
69 } 69 }
70 } 70 }
71 71
72 ////////////////////////////////////////////// 72 ///////////////////////////////////////////////////////
73 G4VParticleChange* G4ContinuousGainOfEnergy::A 73 G4VParticleChange* G4ContinuousGainOfEnergy::AlongStepDoIt(const G4Track& track,
74 74 const G4Step& step)
75 { 75 {
76 // Caution in this method the step length sh 76 // Caution in this method the step length should be the true step length
77 // A problem is that this is computed by the 77 // A problem is that this is computed by the multiple scattering that does
78 // not know the energy at the end of the adj 78 // not know the energy at the end of the adjoint step. This energy is used
79 // during the forward sim. Nothing we can re 79 // during the forward sim. Nothing we can really do against that at this
80 // time. This is inherent to the MS method 80 // time. This is inherent to the MS method
81 81
82 aParticleChange.Initialize(track); 82 aParticleChange.Initialize(track);
83 83
84 // Get the actual (true) Step length 84 // Get the actual (true) Step length
85 G4double length = step.GetStepLength(); 85 G4double length = step.GetStepLength();
86 G4double degain = 0.0; 86 G4double degain = 0.0;
87 87
88 // Compute this for weight change after cont 88 // Compute this for weight change after continuous energy loss
89 G4double DEDX_before = 89 G4double DEDX_before =
90 fDirectEnergyLossProcess->GetDEDX(fPreStep 90 fDirectEnergyLossProcess->GetDEDX(fPreStepKinEnergy, fCurrentCouple);
91 91
92 // For the fluctuation we generate a new dyn 92 // For the fluctuation we generate a new dynamic particle with energy
93 // = preEnergy+egain and then compute the fl 93 // = preEnergy+egain and then compute the fluctuation given in the direct
94 // case. 94 // case.
95 G4DynamicParticle* dynParticle = new G4Dynam 95 G4DynamicParticle* dynParticle = new G4DynamicParticle();
96 *dynParticle = *(track.Get 96 *dynParticle = *(track.GetDynamicParticle());
97 dynParticle->SetDefinition(fDirectPartDef); 97 dynParticle->SetDefinition(fDirectPartDef);
98 G4double Tkin = dynParticle->GetKineticEnerg 98 G4double Tkin = dynParticle->GetKineticEnergy();
99 99
100 G4double dlength = length; 100 G4double dlength = length;
101 if(Tkin != fPreStepKinEnergy && fIsIon) 101 if(Tkin != fPreStepKinEnergy && fIsIon)
102 { 102 {
103 G4double chargeSqRatio = fCurrentModel->Ge 103 G4double chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
104 fDirectPartDef, fCurrentMaterial, Tkin); 104 fDirectPartDef, fCurrentMaterial, Tkin);
105 fDirectEnergyLossProcess->SetDynamicMassCh 105 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio, chargeSqRatio);
106 } 106 }
107 107
108 G4double r = fDirectEnergyLossProcess->GetRa 108 G4double r = fDirectEnergyLossProcess->GetRange(Tkin, fCurrentCouple);
109 if(dlength <= fLinLossLimit * r) 109 if(dlength <= fLinLossLimit * r)
110 { 110 {
111 degain = DEDX_before * dlength; 111 degain = DEDX_before * dlength;
112 } 112 }
113 else 113 else
114 { 114 {
115 G4double x = r + dlength; 115 G4double x = r + dlength;
116 G4double E = fDirectEnergyLossProcess->Get 116 G4double E = fDirectEnergyLossProcess->GetKineticEnergy(x, fCurrentCouple);
117 if(fIsIon) 117 if(fIsIon)
118 { 118 {
119 G4double chargeSqRatio = fCurrentModel-> 119 G4double chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
120 fDirectPartDef, fCurrentMaterial, E); 120 fDirectPartDef, fCurrentMaterial, E);
121 fDirectEnergyLossProcess->SetDynamicMass 121 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio, chargeSqRatio);
122 G4double x1 = fDirectEnergyLossProcess-> 122 G4double x1 = fDirectEnergyLossProcess->GetRange(E, fCurrentCouple);
123 123
124 G4int ii = 0; 124 G4int ii = 0;
125 constexpr G4int iimax = 100; 125 constexpr G4int iimax = 100;
126 while(std::abs(x - x1) > 0.01 * x) 126 while(std::abs(x - x1) > 0.01 * x)
127 { 127 {
128 E = fDirectEnergyLossProcess->GetKinet 128 E = fDirectEnergyLossProcess->GetKineticEnergy(x, fCurrentCouple);
129 chargeSqRatio = fCurrentModel->GetChar 129 chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
130 fDirectPartDef, fCurrentMaterial, E) 130 fDirectPartDef, fCurrentMaterial, E);
131 fDirectEnergyLossProcess->SetDynamicMa 131 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio,
132 132 chargeSqRatio);
133 x1 = fDirectEnergyLossProcess->GetRang 133 x1 = fDirectEnergyLossProcess->GetRange(E, fCurrentCouple);
134 ++ii; 134 ++ii;
135 if(ii >= iimax) 135 if(ii >= iimax)
136 { 136 {
137 break; 137 break;
138 } 138 }
139 } 139 }
140 } 140 }
141 141
142 degain = E - Tkin; 142 degain = E - Tkin;
143 } 143 }
144 G4double tmax = fCurrentModel->MaxSecondaryK 144 G4double tmax = fCurrentModel->MaxSecondaryKinEnergy(dynParticle);
145 fCurrentTcut = std::min(fCurrentTcut, tmax); 145 fCurrentTcut = std::min(fCurrentTcut, tmax);
146 146
147 dynParticle->SetKineticEnergy(Tkin + degain) 147 dynParticle->SetKineticEnergy(Tkin + degain);
148 148
149 // Corrections, which cannot be tabulated fo 149 // Corrections, which cannot be tabulated for ions
150 fCurrentModel->CorrectionsAlongStep(fCurrent 150 fCurrentModel->CorrectionsAlongStep(fCurrentCouple, dynParticle, dlength, degain);
151 151
152 // Sample fluctuations 152 // Sample fluctuations
153 G4double deltaE = 0.; 153 G4double deltaE = 0.;
154 if(fLossFluctuationFlag) 154 if(fLossFluctuationFlag)
155 { 155 {
156 deltaE = fCurrentModel->GetModelOfFluctuat 156 deltaE = fCurrentModel->GetModelOfFluctuations()->SampleFluctuations(
157 fCurrentCouple, dynParticle, fCurrentTcu 157 fCurrentCouple, dynParticle, fCurrentTcut, tmax, dlength, degain)
158 - degain; 158 - degain;
159 } 159 }
160 160
161 G4double egain = degain + deltaE; 161 G4double egain = degain + deltaE;
162 if(egain <= 0.) 162 if(egain <= 0.)
163 egain = degain; 163 egain = degain;
164 Tkin += egain; 164 Tkin += egain;
165 dynParticle->SetKineticEnergy(Tkin); 165 dynParticle->SetKineticEnergy(Tkin);
166 166
167 delete dynParticle; 167 delete dynParticle;
168 168
169 if(fIsIon) 169 if(fIsIon)
170 { 170 {
171 G4double chargeSqRatio = fCurrentModel->Ge 171 G4double chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
172 fDirectPartDef, fCurrentMaterial, Tkin); 172 fDirectPartDef, fCurrentMaterial, Tkin);
173 fDirectEnergyLossProcess->SetDynamicMassCh 173 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio, chargeSqRatio);
174 } 174 }
175 175
176 G4double DEDX_after = fDirectEnergyLossProce 176 G4double DEDX_after = fDirectEnergyLossProcess->GetDEDX(Tkin, fCurrentCouple);
177 G4double weight_correction = DEDX_after / DE 177 G4double weight_correction = DEDX_after / DEDX_before;
178 178
179 aParticleChange.ProposeEnergy(Tkin); 179 aParticleChange.ProposeEnergy(Tkin);
180 180
181 // Caution!!! It is important to select the 181 // Caution!!! It is important to select the weight of the post_step_point
182 // as the current weight and not the weight 182 // as the current weight and not the weight of the track, as the weight of
183 // the track is changed after having applied 183 // the track is changed after having applied all the along_step_do_it.
184 184
185 G4double new_weight = 185 G4double new_weight =
186 weight_correction * step.GetPostStepPoint( 186 weight_correction * step.GetPostStepPoint()->GetWeight();
187 aParticleChange.SetParentWeightByProcess(fal 187 aParticleChange.SetParentWeightByProcess(false);
188 aParticleChange.ProposeParentWeight(new_weig 188 aParticleChange.ProposeParentWeight(new_weight);
189 189
190 return &aParticleChange; 190 return &aParticleChange;
191 } 191 }
192 192
193 ////////////////////////////////////////////// 193 ///////////////////////////////////////////////////////
194 void G4ContinuousGainOfEnergy::SetLossFluctuat 194 void G4ContinuousGainOfEnergy::SetLossFluctuations(G4bool val)
195 { 195 {
196 if(val && !fLossFluctuationArePossible) 196 if(val && !fLossFluctuationArePossible)
197 return; 197 return;
198 fLossFluctuationFlag = val; 198 fLossFluctuationFlag = val;
199 } 199 }
200 200
201 ////////////////////////////////////////////// 201 ///////////////////////////////////////////////////////
202 G4double G4ContinuousGainOfEnergy::GetContinuo 202 G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
203 203 G4double, G4double,
204 204 G4double&)
205 { 205 {
206 DefineMaterial(track.GetMaterialCutsCouple() 206 DefineMaterial(track.GetMaterialCutsCouple());
207 207
208 fPreStepKinEnergy = track.GetKineticEnergy() 208 fPreStepKinEnergy = track.GetKineticEnergy();
209 fCurrentModel = fDirectEnergyLossProcess 209 fCurrentModel = fDirectEnergyLossProcess->SelectModelForMaterial(
210 track.GetKineticEnergy() * fMassRatio, fCu 210 track.GetKineticEnergy() * fMassRatio, fCurrentCoupleIndex);
211 G4double emax_model = fCurrentMode 211 G4double emax_model = fCurrentModel->HighEnergyLimit();
212 G4double preStepChargeSqRatio = 0.; 212 G4double preStepChargeSqRatio = 0.;
213 if(fIsIon) 213 if(fIsIon)
214 { 214 {
215 G4double chargeSqRatio = fCurrentModel->Ge 215 G4double chargeSqRatio = fCurrentModel->GetChargeSquareRatio(
216 fDirectPartDef, fCurrentMaterial, fPreSt 216 fDirectPartDef, fCurrentMaterial, fPreStepKinEnergy);
217 preStepChargeSqRatio = chargeSqRatio; 217 preStepChargeSqRatio = chargeSqRatio;
218 fDirectEnergyLossProcess->SetDynamicMassCh 218 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio,
219 219 preStepChargeSqRatio);
220 } 220 }
221 221
222 G4double maxE = 1.1 * fPreStepKinEnergy; 222 G4double maxE = 1.1 * fPreStepKinEnergy;
223 223
224 if(fPreStepKinEnergy < fCurrentTcut) 224 if(fPreStepKinEnergy < fCurrentTcut)
225 maxE = std::min(fCurrentTcut, maxE); 225 maxE = std::min(fCurrentTcut, maxE);
226 226
227 maxE = std::min(emax_model * 1.001, maxE); 227 maxE = std::min(emax_model * 1.001, maxE);
228 228
229 G4double preStepRange = 229 G4double preStepRange =
230 fDirectEnergyLossProcess->GetRange(fPreSte 230 fDirectEnergyLossProcess->GetRange(fPreStepKinEnergy, fCurrentCouple);
231 231
232 if(fIsIon) 232 if(fIsIon)
233 { 233 {
234 G4double chargeSqRatioAtEmax = fCurrentMod 234 G4double chargeSqRatioAtEmax = fCurrentModel->GetChargeSquareRatio(
235 fDirectPartDef, fCurrentMaterial, maxE); 235 fDirectPartDef, fCurrentMaterial, maxE);
236 fDirectEnergyLossProcess->SetDynamicMassCh 236 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio,
237 237 chargeSqRatioAtEmax);
238 } 238 }
239 239
240 G4double r1 = fDirectEnergyLossProcess->GetR 240 G4double r1 = fDirectEnergyLossProcess->GetRange(maxE, fCurrentCouple);
241 241
242 if(fIsIon) 242 if(fIsIon)
243 fDirectEnergyLossProcess->SetDynamicMassCh 243 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio,
244 244 preStepChargeSqRatio);
245 245
246 return std::max(r1 - preStepRange, 0.001 * m 246 return std::max(r1 - preStepRange, 0.001 * mm);
247 } 247 }
248 248
249 ////////////////////////////////////////////// 249 ///////////////////////////////////////////////////////
250 void G4ContinuousGainOfEnergy::SetDynamicMassC 250 void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track&,
251 251 G4double energy)
252 { 252 {
253 G4double ChargeSqRatio = 253 G4double ChargeSqRatio =
254 G4LossTableManager::Instance()->EmCorrecti 254 G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(
255 fDirectPartDef, fCurrentMaterial, energy 255 fDirectPartDef, fCurrentMaterial, energy);
256 if(fDirectEnergyLossProcess) 256 if(fDirectEnergyLossProcess)
257 fDirectEnergyLossProcess->SetDynamicMassCh 257 fDirectEnergyLossProcess->SetDynamicMassCharge(fMassRatio, ChargeSqRatio);
258 } 258 }
259 259