| 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 75591 2013-11-04 12:33:11Z 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 while (std::abs(x-x1)>0.01*x) {
125 constexpr G4int iimax = 100; << 174 E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
126 while(std::abs(x - x1) > 0.01 * x) << 175 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
127 { << 176 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
128 E = fDirectEnergyLossProcess->GetKinet << 177 x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
129 chargeSqRatio = fCurrentModel->GetChar << 178
130 fDirectPartDef, fCurrentMaterial, E) << 179 }
131 fDirectEnergyLossProcess->SetDynamicMa << 180 }
132 << 181
133 x1 = fDirectEnergyLossProcess->GetRang << 182 degain=E-Tkin;
134 ++ii; << 183
135 if(ii >= iimax) << 184
136 { << 185
137 break; << 186 }
138 } << 187 //G4cout<<degain<<G4endl;
139 } << 188 G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
140 } << 189 tmax = std::min(tmax,currentTcut);
141 << 190
142 degain = E - Tkin; << 191
143 } << 192 dynParticle->SetKineticEnergy(Tkin+degain);
144 G4double tmax = fCurrentModel->MaxSecondaryK << 193
145 fCurrentTcut = std::min(fCurrentTcut, tmax); << 194 // Corrections, which cannot be tabulated for ions
146 << 195 //----------------------------------------
147 dynParticle->SetKineticEnergy(Tkin + degain) << 196 G4double esecdep=0;//not used in most models
148 << 197 currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength);
149 // Corrections, which cannot be tabulated fo << 198
150 fCurrentModel->CorrectionsAlongStep(fCurrent << 199 // Sample fluctuations
151 << 200 //-------------------
152 // Sample fluctuations << 201
153 G4double deltaE = 0.; << 202
154 if(fLossFluctuationFlag) << 203 G4double deltaE =0.;
155 { << 204 if (lossFluctuationFlag ) {
156 deltaE = fCurrentModel->GetModelOfFluctuat << 205 deltaE = currentModel->GetModelOfFluctuations()->
157 fCurrentCouple, dynParticle, fCurrentTcu << 206 SampleFluctuations(currentCouple,dynParticle,tmax,dlength,degain)-degain;
158 - degain; << 207 }
159 } << 208
160 << 209 G4double egain=degain+deltaE;
161 G4double egain = degain + deltaE; << 210 if (egain <=0) egain=degain;
162 if(egain <= 0.) << 211 Tkin+=egain;
163 egain = degain; << 212 dynParticle->SetKineticEnergy(Tkin);
164 Tkin += egain; << 213 }
165 dynParticle->SetKineticEnergy(Tkin); << 214
>> 215
>> 216
166 217
>> 218
167 delete dynParticle; 219 delete dynParticle;
168 << 220
169 if(fIsIon) << 221 if (IsIon){
170 { << 222 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
171 G4double chargeSqRatio = fCurrentModel->Ge << 223 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
172 fDirectPartDef, fCurrentMaterial, Tkin); << 224
173 fDirectEnergyLossProcess->SetDynamicMassCh <<
174 } 225 }
175 << 226
176 G4double DEDX_after = fDirectEnergyLossProce << 227 G4double DEDX_after = theDirectEnergyLossProcess->GetDEDX(Tkin, currentCouple);
177 G4double weight_correction = DEDX_after / DE << 228
178 << 229
>> 230 G4double weight_correction=DEDX_after/DEDX_before;
>> 231
>> 232
179 aParticleChange.ProposeEnergy(Tkin); 233 aParticleChange.ProposeEnergy(Tkin);
180 234
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 235
185 G4double new_weight = << 236 //Caution!!!
186 weight_correction * step.GetPostStepPoint( << 237 // It is important to select the weight of the post_step_point
>> 238 // as the current weight and not the weight of the track, as t
>> 239 // the weight of the track is changed after having applied all
>> 240 // the along_step_do_it.
>> 241
>> 242 // G4double new_weight=weight_correction*track.GetWeight(); //old
>> 243 G4double new_weight=weight_correction*step.GetPostStepPoint()->GetWeight();
187 aParticleChange.SetParentWeightByProcess(fal 244 aParticleChange.SetParentWeightByProcess(false);
188 aParticleChange.ProposeParentWeight(new_weig 245 aParticleChange.ProposeParentWeight(new_weight);
189 246
>> 247
190 return &aParticleChange; 248 return &aParticleChange;
191 } <<
192 249
>> 250 }
193 ////////////////////////////////////////////// 251 ///////////////////////////////////////////////////////
>> 252 //
194 void G4ContinuousGainOfEnergy::SetLossFluctuat 253 void G4ContinuousGainOfEnergy::SetLossFluctuations(G4bool val)
195 { 254 {
196 if(val && !fLossFluctuationArePossible) << 255 if(val && !lossFluctuationArePossible) return;
197 return; << 256 lossFluctuationFlag = val;
198 fLossFluctuationFlag = val; <<
199 } 257 }
200 <<
201 ////////////////////////////////////////////// 258 ///////////////////////////////////////////////////////
202 G4double G4ContinuousGainOfEnergy::GetContinuo << 259 //
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 260
240 G4double r1 = fDirectEnergyLossProcess->GetR <<
241 261
242 if(fIsIon) <<
243 fDirectEnergyLossProcess->SetDynamicMassCh <<
244 <<
245 262
246 return std::max(r1 - preStepRange, 0.001 * m << 263 G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
>> 264 G4double , G4double , G4double& )
>> 265 {
>> 266 G4double x = DBL_MAX;
>> 267 x=.1*mm;
>> 268
>> 269
>> 270 DefineMaterial(track.GetMaterialCutsCouple());
>> 271
>> 272 preStepKinEnergy = track.GetKineticEnergy();
>> 273 preStepScaledKinEnergy = track.GetKineticEnergy()*massRatio;
>> 274 currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(preStepScaledKinEnergy,currentCoupleIndex);
>> 275 G4double emax_model=currentModel->HighEnergyLimit();
>> 276 if (IsIon) {
>> 277 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,preStepKinEnergy);
>> 278 preStepChargeSqRatio = chargeSqRatio;
>> 279 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
>> 280 }
>> 281
>> 282
>> 283 G4double maxE =1.1*preStepKinEnergy;
>> 284 /*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy;
>> 285 else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy;
>> 286 else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/
>> 287
>> 288 if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE);
>> 289
>> 290 maxE=std::min(emax_model*1.001,maxE);
>> 291
>> 292 preStepRange = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple);
>> 293
>> 294 if (IsIon) {
>> 295 G4double chargeSqRatioAtEmax = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,maxE);
>> 296 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatioAtEmax);
>> 297 }
>> 298
>> 299 G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple);
>> 300
>> 301 if (IsIon) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
>> 302
>> 303
>> 304
>> 305 x=r1-preStepRange;
>> 306 x=std::max(r1-preStepRange,0.001*mm);
>> 307
>> 308 return x;
>> 309
>> 310
247 } 311 }
248 << 312 #include "G4EmCorrections.hh"
249 ////////////////////////////////////////////// 313 ///////////////////////////////////////////////////////
250 void G4ContinuousGainOfEnergy::SetDynamicMassC << 314 //
251 << 315
252 { << 316 void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track& ,G4double energy)
253 G4double ChargeSqRatio = << 317 {
254 G4LossTableManager::Instance()->EmCorrecti << 318
255 fDirectPartDef, fCurrentMaterial, energy << 319 G4double ChargeSqRatio= G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(theDirectPartDef,currentMaterial,energy);
256 if(fDirectEnergyLossProcess) << 320 if (theDirectEnergyLossProcess) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,ChargeSqRatio);
257 fDirectEnergyLossProcess->SetDynamicMassCh <<
258 } 321 }
259 322