ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ContinuousGainOfEnergy.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ContinuousGainOfEnergy.cc
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 
29 
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4Step.hh"
33 #include "G4ParticleDefinition.hh"
34 #include "G4VEmModel.hh"
35 #include "G4VEmFluctuationModel.hh"
36 #include "G4VParticleChange.hh"
37 #include "G4AdjointCSManager.hh"
38 #include "G4LossTableManager.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4PhysicalConstants.hh"
41 
43 //
45  G4ProcessType type): G4VContinuousProcess(name, type)
46 {
47 
48 
49  linLossLimit=0.05;
52  is_integral = false;
53 
54  //Will be properly set in SetDirectParticle()
55  IsIon=false;
56  massRatio =1.;
57  chargeSqRatio=1.;
59 
60  //Some initialization
61  currentCoupleIndex=9999999;
63  currentMaterialIndex=9999999;
64  currentTcut=0.;
66  preStepRange=0.;
68 
69  currentCouple=0;
70 }
71 
73 //
75 {
76 
77 }
79 //
80 
82  const G4ParticleDefinition& )
83 {//theDirectEnergyLossProcess->PreparePhysicsTable(part);
84 
85 ;
86 }
87 
89 //
90 
92 {//theDirectEnergyLossProcess->BuildPhysicsTable(part);
93 ;
94 }
95 
97 //
100  if (theDirectPartDef->GetParticleType()== "nucleus") {
101  IsIon=true;
104  chargeSqRatio=q*q;
105 
106 
107  }
108 
109 }
110 
112 //
113 //
115  const G4Step& step)
116 {
117 
118  //Caution in this method the step length should be the true step length
119  // 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
120  //Forward sim. Nothing we can really do against that at this time. This is inherent to the MS method
121  //
122 
123 
124 
126 
127  // Get the actual (true) Step length
128  //----------------------------------
129  G4double length = step.GetStepLength();
130  G4double degain = 0.0;
131 
132 
133 
134  // Compute this for weight change after continuous energy loss
135  //-------------------------------------------------------------
137 
138 
139 
140  // For the fluctuation we generate a new dynamic particle with energy =preEnergy+egain
141  // and then compute the fluctuation given in the direct case.
142  //-----------------------------------------------------------------------
143  G4DynamicParticle* dynParticle = new G4DynamicParticle();
144  *dynParticle = *(track.GetDynamicParticle());
145  dynParticle->SetDefinition(theDirectPartDef);
146  G4double Tkin = dynParticle->GetKineticEnergy();
147 
148 
149  size_t n=1;
150  if (is_integral ) n=10;
151  n=1;
152  G4double dlength= length/n;
153  for (size_t i=0;i<n;i++) {
154  if (Tkin != preStepKinEnergy && IsIon) {
157 
158  }
159 
161  if( dlength <= linLossLimit * r ) {
162  degain = DEDX_before*dlength;
163  }
164  else {
165  G4double x = r + dlength;
166  //degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
168  if (IsIon){
172 
173  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
174  G4int ii=0;
175  const G4int iimax = 100;
176  while (std::abs(x-x1)>0.01*x) {
181  ++ii;
182  if(ii >= iimax) { break; }
183  }
184  }
185 
186  degain=E-Tkin;
187 
188 
189 
190  }
191  //G4cout<<degain<<G4endl;
192  G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
193  tmax = std::min(tmax,currentTcut);
194 
195 
196  dynParticle->SetKineticEnergy(Tkin+degain);
197 
198  // Corrections, which cannot be tabulated for ions
199  //----------------------------------------
200  G4double esecdep=0;//not used in most models
201  currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength);
202 
203  // Sample fluctuations
204  //-------------------
205 
206 
207  G4double deltaE =0.;
208  if (lossFluctuationFlag ) {
210  SampleFluctuations(currentCouple,dynParticle,tmax,dlength,degain)-degain;
211  }
212 
213  G4double egain=degain+deltaE;
214  if (egain <=0) egain=degain;
215  Tkin+=egain;
216  dynParticle->SetKineticEnergy(Tkin);
217  }
218 
219 
220 
221 
222 
223  delete dynParticle;
224 
225  if (IsIon){
228 
229  }
230 
232 
233 
234  G4double weight_correction=DEDX_after/DEDX_before;
235 
236 
238 
239 
240  //Caution!!!
241  // It is important to select the weight of the post_step_point
242  // as the current weight and not the weight of the track, as t
243  // the weight of the track is changed after having applied all
244  // the along_step_do_it.
245 
246  // G4double new_weight=weight_correction*track.GetWeight(); //old
247  G4double new_weight=weight_correction*step.GetPostStepPoint()->GetWeight();
250 
251 
252  return &aParticleChange;
253 
254 }
256 //
258 {
259  if(val && !lossFluctuationArePossible) return;
260  lossFluctuationFlag = val;
261 }
263 //
264 
265 
266 
268  G4double , G4double , G4double& )
269 {
270  G4double x = DBL_MAX;
271  x=.1*mm;
272 
273 
275 
279  G4double emax_model=currentModel->HighEnergyLimit();
280  if (IsIon) {
284  }
285 
286 
288  /*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy;
289  else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy;
290  else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/
291 
292  if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE);
293 
294  maxE=std::min(emax_model*1.001,maxE);
295 
297 
298  if (IsIon) {
301  }
302 
304 
306 
307 
308 
309  x=r1-preStepRange;
310  x=std::max(r1-preStepRange,0.001*mm);
311 
312  return x;
313 
314 
315 }
316 #include "G4EmCorrections.hh"
318 //
319 
321 {
322 
325 }