Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/adjoint/src/G4ContinuousGainOfEnergy.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/adjoint/src/G4ContinuousGainOfEnergy.cc (Version 11.3.0) and /processes/electromagnetic/adjoint/src/G4ContinuousGainOfEnergy.cc (Version 11.2)


  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