| 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 // $Id: G4ContinuousGainOfEnergy.cc 91870 2015-08-07 15:21:40Z gcosmo $
26 // 27 //
27 28
28 #include "G4ContinuousGainOfEnergy.hh" 29 #include "G4ContinuousGainOfEnergy.hh"
29 30
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" 31 #include "G4PhysicalConstants.hh"
37 #include "G4Step.hh" <<
38 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
39 #include "G4VEmFluctuationModel.hh" << 33 #include "G4Step.hh"
>> 34 #include "G4ParticleDefinition.hh"
40 #include "G4VEmModel.hh" 35 #include "G4VEmModel.hh"
41 #include "G4VEnergyLossProcess.hh" << 36 #include "G4VEmFluctuationModel.hh"
42 #include "G4VParticleChange.hh" 37 #include "G4VParticleChange.hh"
>> 38 #include "G4AdjointCSManager.hh"
>> 39 #include "G4LossTableManager.hh"
>> 40 #include "G4SystemOfUnits.hh"
>> 41 #include "G4PhysicalConstants.hh"
43 42
44 ////////////////////////////////////////////// 43 ///////////////////////////////////////////////////////
45 G4ContinuousGainOfEnergy::G4ContinuousGainOfEn << 44 //
46 << 45 G4ContinuousGainOfEnergy::G4ContinuousGainOfEnergy(const G4String& name,
47 : G4VContinuousProcess(name, type) << 46 G4ProcessType type): G4VContinuousProcess(name, type)
48 {} << 47 {
49 48
50 ////////////////////////////////////////////// << 49
51 G4ContinuousGainOfEnergy::~G4ContinuousGainOfE << 50 linLossLimit=0.05;
>> 51 lossFluctuationArePossible =true;
>> 52 lossFluctuationFlag=true;
>> 53 is_integral = false;
>> 54
>> 55 //Will be properly set in SetDirectParticle()
>> 56 IsIon=false;
>> 57 massRatio =1.;
>> 58 chargeSqRatio=1.;
>> 59 preStepChargeSqRatio=1.;
>> 60
>> 61 //Some initialization
>> 62 currentCoupleIndex=9999999;
>> 63 currentCutInRange=0.;
>> 64 currentMaterialIndex=9999999;
>> 65 currentTcut=0.;
>> 66 preStepKinEnergy=0.;
>> 67 preStepRange=0.;
>> 68 preStepScaledKinEnergy=0.;
>> 69
>> 70 currentCouple=0;
>> 71 }
52 72
53 ////////////////////////////////////////////// 73 ///////////////////////////////////////////////////////
54 void G4ContinuousGainOfEnergy::ProcessDescript << 74 //
>> 75 G4ContinuousGainOfEnergy::~G4ContinuousGainOfEnergy()
55 { 76 {
56 out << "Continuous process acting on adjoint << 77
57 "continuous gain of energy of charged << 78 }
58 "tracked back.\n"; << 79 ///////////////////////////////////////////////////////
>> 80 //
>> 81
>> 82 void G4ContinuousGainOfEnergy::PreparePhysicsTable(
>> 83 const G4ParticleDefinition& )
>> 84 {//theDirectEnergyLossProcess->PreparePhysicsTable(part);
>> 85
>> 86 ;
59 } 87 }
60 88
61 ////////////////////////////////////////////// 89 ///////////////////////////////////////////////////////
62 void G4ContinuousGainOfEnergy::SetDirectPartic << 90 //
63 { << 91
64 fDirectPartDef = p; << 92 void G4ContinuousGainOfEnergy::BuildPhysicsTable(const G4ParticleDefinition&)
65 if(fDirectPartDef->GetParticleType() == "nuc << 93 {//theDirectEnergyLossProcess->BuildPhysicsTable(part);
66 { << 94 ;
67 fIsIon = true; << 95 }
68 fMassRatio = proton_mass_c2 / fDirectPartD << 96
69 } << 97 ///////////////////////////////////////////////////////
>> 98 //
>> 99 void G4ContinuousGainOfEnergy::SetDirectParticle(G4ParticleDefinition* p)
>> 100 {theDirectPartDef=p;
>> 101 if (theDirectPartDef->GetParticleType()== "nucleus") {
>> 102 IsIon=true;
>> 103 massRatio = proton_mass_c2/theDirectPartDef->GetPDGMass();
>> 104 G4double q=theDirectPartDef->GetPDGCharge();
>> 105 chargeSqRatio=q*q;
>> 106
>> 107
>> 108 }
>> 109
70 } 110 }
71 111
72 ////////////////////////////////////////////// 112 ///////////////////////////////////////////////////////
>> 113 //
>> 114 //
73 G4VParticleChange* G4ContinuousGainOfEnergy::A 115 G4VParticleChange* G4ContinuousGainOfEnergy::AlongStepDoIt(const G4Track& track,
74 << 116 const G4Step& step)
75 { 117 {
76 // Caution in this method the step length sh << 118
77 // A problem is that this is computed by the << 119 //Caution in this method the step length should be the true step length
78 // not know the energy at the end of the adj << 120 // A problem is that this is compute by the multiple scattering that does not know the energy at the end of the adjoint step. This energy is used during the
79 // during the forward sim. Nothing we can re << 121 //Forward sim. Nothing we can really do against that at this time. This is inherent to the MS method
80 // time. This is inherent to the MS method << 122 //
81 << 123
>> 124
>> 125
82 aParticleChange.Initialize(track); 126 aParticleChange.Initialize(track);
83 << 127
84 // Get the actual (true) Step length 128 // Get the actual (true) Step length
>> 129 //----------------------------------
85 G4double length = step.GetStepLength(); 130 G4double length = step.GetStepLength();
86 G4double degain = 0.0; << 131 G4double degain = 0.0;
87 << 132
>> 133
>> 134
88 // Compute this for weight change after cont 135 // Compute this for weight change after continuous energy loss
89 G4double DEDX_before = << 136 //-------------------------------------------------------------
90 fDirectEnergyLossProcess->GetDEDX(fPreStep << 137 G4double DEDX_before = theDirectEnergyLossProcess->GetDEDX(preStepKinEnergy, currentCouple);
91 << 138
92 // For the fluctuation we generate a new dyn << 139
93 // = preEnergy+egain and then compute the fl << 140
94 // case. << 141 // For the fluctuation we generate a new dynamic particle with energy =preEnergy+egain
>> 142 // and then compute the fluctuation given in the direct case.
>> 143 //-----------------------------------------------------------------------
95 G4DynamicParticle* dynParticle = new G4Dynam 144 G4DynamicParticle* dynParticle = new G4DynamicParticle();
96 *dynParticle = *(track.Get << 145 *dynParticle = *(track.GetDynamicParticle());
97 dynParticle->SetDefinition(fDirectPartDef); << 146 dynParticle->SetDefinition(theDirectPartDef);
98 G4double Tkin = dynParticle->GetKineticEnerg << 147 G4double Tkin = dynParticle->GetKineticEnergy();
99 << 148
100 G4double dlength = length; << 149
101 if(Tkin != fPreStepKinEnergy && fIsIon) << 150 size_t n=1;
102 { << 151 if (is_integral ) n=10;
103 G4double chargeSqRatio = fCurrentModel->Ge << 152 n=1;
104 fDirectPartDef, fCurrentMaterial, Tkin); << 153 G4double dlength= length/n;
105 fDirectEnergyLossProcess->SetDynamicMassCh << 154 for (size_t i=0;i<n;i++) {
106 } << 155 if (Tkin != preStepKinEnergy && IsIon) {
107 << 156 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
108 G4double r = fDirectEnergyLossProcess->GetRa << 157 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
109 if(dlength <= fLinLossLimit * r) << 158
110 { << 159 }
111 degain = DEDX_before * dlength; << 160
112 } << 161 G4double r = theDirectEnergyLossProcess->GetRange(Tkin, currentCouple);
113 else << 162 if( dlength <= linLossLimit * r ) {
114 { << 163 degain = DEDX_before*dlength;
115 G4double x = r + dlength; << 164 }
116 G4double E = fDirectEnergyLossProcess->Get << 165 else {
117 if(fIsIon) << 166 G4double x = r + dlength;
118 { << 167 //degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
119 G4double chargeSqRatio = fCurrentModel-> << 168 G4double E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
120 fDirectPartDef, fCurrentMaterial, E); << 169 if (IsIon){
121 fDirectEnergyLossProcess->SetDynamicMass << 170 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
122 G4double x1 = fDirectEnergyLossProcess-> << 171 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
123 << 172 G4double x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
124 G4int ii = 0; << 173
125 constexpr G4int iimax = 100; << 174 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
126 while(std::abs(x - x1) > 0.01 * x) << 175 G4int ii=0;
127 { << 176 const G4int iimax = 100;
128 E = fDirectEnergyLossProcess->GetKinet << 177 while (std::abs(x-x1)>0.01*x) {
129 chargeSqRatio = fCurrentModel->GetChar << 178 E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
130 fDirectPartDef, fCurrentMaterial, E) << 179 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
131 fDirectEnergyLossProcess->SetDynamicMa << 180 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
132 << 181 x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
133 x1 = fDirectEnergyLossProcess->GetRang << 182 ++ii;
134 ++ii; << 183 if(ii >= iimax) { break; }
135 if(ii >= iimax) << 184 }
136 { << 185 }
137 break; << 186
138 } << 187 degain=E-Tkin;
139 } << 188
140 } << 189
141 << 190
142 degain = E - Tkin; << 191 }
143 } << 192 //G4cout<<degain<<G4endl;
144 G4double tmax = fCurrentModel->MaxSecondaryK << 193 G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
145 fCurrentTcut = std::min(fCurrentTcut, tmax); << 194 tmax = std::min(tmax,currentTcut);
146 << 195
147 dynParticle->SetKineticEnergy(Tkin + degain) << 196
148 << 197 dynParticle->SetKineticEnergy(Tkin+degain);
149 // Corrections, which cannot be tabulated fo << 198
150 fCurrentModel->CorrectionsAlongStep(fCurrent << 199 // Corrections, which cannot be tabulated for ions
151 << 200 //----------------------------------------
152 // Sample fluctuations << 201 G4double esecdep=0;//not used in most models
153 G4double deltaE = 0.; << 202 currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength);
154 if(fLossFluctuationFlag) << 203
155 { << 204 // Sample fluctuations
156 deltaE = fCurrentModel->GetModelOfFluctuat << 205 //-------------------
157 fCurrentCouple, dynParticle, fCurrentTcu << 206
158 - degain; << 207
159 } << 208 G4double deltaE =0.;
160 << 209 if (lossFluctuationFlag ) {
161 G4double egain = degain + deltaE; << 210 deltaE = currentModel->GetModelOfFluctuations()->
162 if(egain <= 0.) << 211 SampleFluctuations(currentCouple,dynParticle,tmax,dlength,degain)-degain;
163 egain = degain; << 212 }
164 Tkin += egain; << 213
165 dynParticle->SetKineticEnergy(Tkin); << 214 G4double egain=degain+deltaE;
>> 215 if (egain <=0) egain=degain;
>> 216 Tkin+=egain;
>> 217 dynParticle->SetKineticEnergy(Tkin);
>> 218 }
>> 219
>> 220
>> 221
166 222
>> 223
167 delete dynParticle; 224 delete dynParticle;
168 << 225
169 if(fIsIon) << 226 if (IsIon){
170 { << 227 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
171 G4double chargeSqRatio = fCurrentModel->Ge << 228 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
172 fDirectPartDef, fCurrentMaterial, Tkin); << 229
173 fDirectEnergyLossProcess->SetDynamicMassCh <<
174 } 230 }
175 << 231
176 G4double DEDX_after = fDirectEnergyLossProce << 232 G4double DEDX_after = theDirectEnergyLossProcess->GetDEDX(Tkin, currentCouple);
177 G4double weight_correction = DEDX_after / DE << 233
178 << 234
>> 235 G4double weight_correction=DEDX_after/DEDX_before;
>> 236
>> 237
179 aParticleChange.ProposeEnergy(Tkin); 238 aParticleChange.ProposeEnergy(Tkin);
180 239
181 // Caution!!! It is important to select the <<
182 // as the current weight and not the weight <<
183 // the track is changed after having applied <<
184 240
185 G4double new_weight = << 241 //Caution!!!
186 weight_correction * step.GetPostStepPoint( << 242 // It is important to select the weight of the post_step_point
>> 243 // as the current weight and not the weight of the track, as t
>> 244 // the weight of the track is changed after having applied all
>> 245 // the along_step_do_it.
>> 246
>> 247 // G4double new_weight=weight_correction*track.GetWeight(); //old
>> 248 G4double new_weight=weight_correction*step.GetPostStepPoint()->GetWeight();
187 aParticleChange.SetParentWeightByProcess(fal 249 aParticleChange.SetParentWeightByProcess(false);
188 aParticleChange.ProposeParentWeight(new_weig 250 aParticleChange.ProposeParentWeight(new_weight);
189 251
>> 252
190 return &aParticleChange; 253 return &aParticleChange;
191 } <<
192 254
>> 255 }
193 ////////////////////////////////////////////// 256 ///////////////////////////////////////////////////////
>> 257 //
194 void G4ContinuousGainOfEnergy::SetLossFluctuat 258 void G4ContinuousGainOfEnergy::SetLossFluctuations(G4bool val)
195 { 259 {
196 if(val && !fLossFluctuationArePossible) << 260 if(val && !lossFluctuationArePossible) return;
197 return; << 261 lossFluctuationFlag = val;
198 fLossFluctuationFlag = val; <<
199 } 262 }
200 <<
201 ////////////////////////////////////////////// 263 ///////////////////////////////////////////////////////
202 G4double G4ContinuousGainOfEnergy::GetContinuo << 264 //
203 <<
204 <<
205 { <<
206 DefineMaterial(track.GetMaterialCutsCouple() <<
207 <<
208 fPreStepKinEnergy = track.GetKineticEnergy() <<
209 fCurrentModel = fDirectEnergyLossProcess <<
210 track.GetKineticEnergy() * fMassRatio, fCu <<
211 G4double emax_model = fCurrentMode <<
212 G4double preStepChargeSqRatio = 0.; <<
213 if(fIsIon) <<
214 { <<
215 G4double chargeSqRatio = fCurrentModel->Ge <<
216 fDirectPartDef, fCurrentMaterial, fPreSt <<
217 preStepChargeSqRatio = chargeSqRatio; <<
218 fDirectEnergyLossProcess->SetDynamicMassCh <<
219 <<
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(fPreSte <<
231 <<
232 if(fIsIon) <<
233 { <<
234 G4double chargeSqRatioAtEmax = fCurrentMod <<
235 fDirectPartDef, fCurrentMaterial, maxE); <<
236 fDirectEnergyLossProcess->SetDynamicMassCh <<
237 <<
238 } <<
239 265
240 G4double r1 = fDirectEnergyLossProcess->GetR <<
241 266
242 if(fIsIon) <<
243 fDirectEnergyLossProcess->SetDynamicMassCh <<
244 <<
245 267
246 return std::max(r1 - preStepRange, 0.001 * m << 268 G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
>> 269 G4double , G4double , G4double& )
>> 270 {
>> 271 G4double x = DBL_MAX;
>> 272 x=.1*mm;
>> 273
>> 274
>> 275 DefineMaterial(track.GetMaterialCutsCouple());
>> 276
>> 277 preStepKinEnergy = track.GetKineticEnergy();
>> 278 preStepScaledKinEnergy = track.GetKineticEnergy()*massRatio;
>> 279 currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(preStepScaledKinEnergy,currentCoupleIndex);
>> 280 G4double emax_model=currentModel->HighEnergyLimit();
>> 281 if (IsIon) {
>> 282 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,preStepKinEnergy);
>> 283 preStepChargeSqRatio = chargeSqRatio;
>> 284 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
>> 285 }
>> 286
>> 287
>> 288 G4double maxE =1.1*preStepKinEnergy;
>> 289 /*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy;
>> 290 else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy;
>> 291 else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/
>> 292
>> 293 if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE);
>> 294
>> 295 maxE=std::min(emax_model*1.001,maxE);
>> 296
>> 297 preStepRange = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple);
>> 298
>> 299 if (IsIon) {
>> 300 G4double chargeSqRatioAtEmax = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,maxE);
>> 301 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatioAtEmax);
>> 302 }
>> 303
>> 304 G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple);
>> 305
>> 306 if (IsIon) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
>> 307
>> 308
>> 309
>> 310 x=r1-preStepRange;
>> 311 x=std::max(r1-preStepRange,0.001*mm);
>> 312
>> 313 return x;
>> 314
>> 315
247 } 316 }
248 << 317 #include "G4EmCorrections.hh"
249 ////////////////////////////////////////////// 318 ///////////////////////////////////////////////////////
250 void G4ContinuousGainOfEnergy::SetDynamicMassC << 319 //
251 << 320
252 { << 321 void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track& ,G4double energy)
253 G4double ChargeSqRatio = << 322 {
254 G4LossTableManager::Instance()->EmCorrecti << 323
255 fDirectPartDef, fCurrentMaterial, energy << 324 G4double ChargeSqRatio= G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(theDirectPartDef,currentMaterial,energy);
256 if(fDirectEnergyLossProcess) << 325 if (theDirectEnergyLossProcess) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,ChargeSqRatio);
257 fDirectEnergyLossProcess->SetDynamicMassCh <<
258 } 326 }
259 327