Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4hCoulombScatteringModel.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/standard/src/G4hCoulombScatteringModel.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4hCoulombScatteringModel.cc (Version 10.1.p1)


  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: G4hCoulombScatteringModel.cc 81579 2014-06-03 10:15:54Z gcosmo $
                                                   >>  27 //
 26 // -------------------------------------------     28 // -------------------------------------------------------------------
 27 //                                                 29 //
 28 // GEANT4 Class file                               30 // GEANT4 Class file
 29 //                                                 31 //
 30 //                                                 32 //
 31 // File name:     G4hCoulombScatteringModel        33 // File name:     G4hCoulombScatteringModel
 32 //                                                 34 //
 33 // Author:        Vladimir Ivanchenko              35 // Author:        Vladimir Ivanchenko 
 34 //                                                 36 //
 35 // Creation date: 08.06.2012 from G4eCoulombSc     37 // Creation date: 08.06.2012 from G4eCoulombScatteringModel
 36 //                                                 38 //
 37 // Modifications:                                  39 // Modifications:
 38 //                                                 40 //
 39 //                                                 41 //
 40 // Class Description:                              42 // Class Description:
 41 //                                                 43 //
 42 // -------------------------------------------     44 // -------------------------------------------------------------------
 43 //                                                 45 //
 44 //....oooOO0OOooo........oooOO0OOooo........oo     46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45 //....oooOO0OOooo........oooOO0OOooo........oo     47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 46                                                    48 
 47 #include "G4hCoulombScatteringModel.hh"            49 #include "G4hCoulombScatteringModel.hh"
 48 #include "G4PhysicalConstants.hh"                  50 #include "G4PhysicalConstants.hh"
 49 #include "G4SystemOfUnits.hh"                      51 #include "G4SystemOfUnits.hh"
 50 #include "Randomize.hh"                            52 #include "Randomize.hh"
 51 #include "G4DataVector.hh"                         53 #include "G4DataVector.hh"
 52 #include "G4ElementTable.hh"                       54 #include "G4ElementTable.hh"
 53 #include "G4ParticleChangeForGamma.hh"             55 #include "G4ParticleChangeForGamma.hh"
 54 #include "G4Proton.hh"                             56 #include "G4Proton.hh"
 55 #include "G4ParticleTable.hh"                      57 #include "G4ParticleTable.hh"
 56 #include "G4IonTable.hh"                           58 #include "G4IonTable.hh"
 57 #include "G4ProductionCutsTable.hh"                59 #include "G4ProductionCutsTable.hh"
 58 #include "G4NucleiProperties.hh"                   60 #include "G4NucleiProperties.hh"
 59 #include "G4Pow.hh"                                61 #include "G4Pow.hh"
                                                   >>  62 #include "G4LossTableManager.hh"
                                                   >>  63 #include "G4LossTableBuilder.hh"
 60 #include "G4NistManager.hh"                        64 #include "G4NistManager.hh"
 61                                                    65 
 62 //....oooOO0OOooo........oooOO0OOooo........oo     66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 63                                                    67 
                                                   >>  68 using namespace std;
                                                   >>  69 
 64 G4hCoulombScatteringModel::G4hCoulombScatterin     70 G4hCoulombScatteringModel::G4hCoulombScatteringModel(G4bool combined)
 65   : G4VEmModel("hCoulombScattering"),          <<  71   : G4VEmModel("eCoulombScattering"),
 66     cosThetaMin(1.0),                              72     cosThetaMin(1.0),
 67     cosThetaMax(-1.0),                             73     cosThetaMax(-1.0),
 68     isCombined(combined)                       <<  74     isCombined(combined),
                                                   >>  75     isInitialised(false)
 69 {                                                  76 {
 70   fParticleChange = nullptr;                   <<  77   fParticleChange = 0;
 71   fNistManager = G4NistManager::Instance();        78   fNistManager = G4NistManager::Instance();
 72   theIonTable  = G4ParticleTable::GetParticleT <<  79   theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
 73   theProton    = G4Proton::Proton();           <<  80   theProton   = G4Proton::Proton();
 74   currentMaterial = nullptr;                   <<  81   currentMaterial = 0; 
 75   fixedCut = -1.0;                             <<  82 
 76                                                <<  83   pCuts = 0;
 77   pCuts = nullptr;                             <<  84 
 78                                                <<  85   lowEnergyThreshold = 1*keV;  // particle will be killed for lower energy
 79   recoilThreshold = 0.0; // by default does no <<  86   recoilThreshold = 0.*keV; // by default does not work
 80                                                <<  87 
 81   particle = nullptr;                          <<  88   particle = 0;
 82   currentCouple = nullptr;                     <<  89   currentCouple = 0;
 83   wokvi = new G4WentzelVIRelXSection();        <<  90   wokvi = new G4WentzelVIRelXSection(combined);
 84                                                    91 
 85   currentMaterialIndex = 0;                        92   currentMaterialIndex = 0;
 86   mass = CLHEP::proton_mass_c2;                <<  93 
                                                   >>  94   cosTetMinNuc = 1.0;
                                                   >>  95   cosTetMaxNuc = -1.0;
 87   elecRatio = 0.0;                                 96   elecRatio = 0.0;
                                                   >>  97   mass = proton_mass_c2;
 88 }                                                  98 }
 89                                                    99 
 90 //....oooOO0OOooo........oooOO0OOooo........oo    100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 91                                                   101 
 92 G4hCoulombScatteringModel::~G4hCoulombScatteri    102 G4hCoulombScatteringModel::~G4hCoulombScatteringModel()
 93 {                                                 103 {
 94   delete wokvi;                                   104   delete wokvi;
 95 }                                                 105 }
 96                                                   106 
 97 //....oooOO0OOooo........oooOO0OOooo........oo    107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 98                                                   108 
 99 void G4hCoulombScatteringModel::Initialise(con << 109 void G4hCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
100              const G4DataVector& cuts)            110              const G4DataVector& cuts)
101 {                                                 111 {
102   SetupParticle(part);                         << 112   SetupParticle(p);
103   currentCouple = nullptr;                     << 113   currentCouple = 0;
104                                                   114 
105   // defined theta limit between single and mu << 115   if(isCombined) {
106   isCombined = true;                           << 116     cosThetaMin = 1.0;
107   G4double tet = PolarAngleLimit();            << 117     G4double tet = PolarAngleLimit();
108                                                << 118     if(tet > pi)       { cosThetaMin = -1.0; }
109   if(tet <= 0.0) {                             << 119     else if(tet > 0.0) { cosThetaMin = cos(tet); }
110     cosThetaMin = 1.0;                         << 
111     isCombined = false;                        << 
112   } else if(tet >= CLHEP::pi) {                << 
113     cosThetaMin = -1.0;                        << 
114   } else {                                     << 
115     cosThetaMin = std::cos(tet);               << 
116   }                                               120   }
117                                                   121 
118   wokvi->Initialise(part, cosThetaMin);        << 122   wokvi->Initialise(p, cosThetaMin);
119   /*                                           << 123   /*  
120   G4cout << "G4hCoulombScatteringModel: " << p    124   G4cout << "G4hCoulombScatteringModel: " << particle->GetParticleName()
121          << "  1-cos(ThetaLimit)= " << 1 - cos    125          << "  1-cos(ThetaLimit)= " << 1 - cosThetaMin
122    << "  cos(thetaMax)= " <<  cosThetaMax         126    << "  cos(thetaMax)= " <<  cosThetaMax
123    << G4endl;                                     127    << G4endl;
124   */                                              128   */
125   pCuts = &cuts;                                  129   pCuts = &cuts;
126   //G4ProductionCutsTable::GetProductionCutsTa << 130   //  G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
127   /*                                           << 131   //G4cout << "!!! G4hCoulombScatteringModel::Initialise for " 
128   G4cout << "!!! G4hCoulombScatteringModel::In << 132   //   << p->GetParticleName() << "  cos(TetMin)= " << cosThetaMin 
129      << part->GetParticleName() << "  cos(TetM << 133   //   << "  cos(TetMax)= " << cosThetaMax <<G4endl;
130      << "  cos(TetMax)= " << cosThetaMax <<G4e << 134   // G4cout << "cut0= " << cuts[0] << "  cut1= " << cuts[1] << G4endl;
131   G4cout << "cut= " << (*pCuts)[0] << "  cut1= << 135   if(!isInitialised) {
132   */                                           << 136     isInitialised = true;
133   if(!fParticleChange) {                       << 
134     fParticleChange = GetParticleChangeForGamm    137     fParticleChange = GetParticleChangeForGamma();
135   }                                               138   }
136   if(IsMaster() && mass < CLHEP::GeV && part-> << 139   if(mass < GeV) {
137     InitialiseElementSelectors(part, cuts);    << 140     InitialiseElementSelectors(p,cuts);
138   }                                            << 
139 }                                              << 
140                                                << 
141 //....oooOO0OOooo........oooOO0OOooo........oo << 
142                                                << 
143 void G4hCoulombScatteringModel::InitialiseLoca << 
144             G4VEmModel* masterModel)           << 
145 {                                              << 
146   SetElementSelectors(masterModel->GetElementS << 
147 }                                              << 
148                                                << 
149 //....oooOO0OOooo........oooOO0OOooo........oo << 
150                                                << 
151 G4double                                       << 
152 G4hCoulombScatteringModel::MinPrimaryEnergy(co << 
153               const G4ParticleDefinition* part << 
154               G4double)                        << 
155 {                                              << 
156   SetupParticle(part);                         << 
157                                                << 
158   // define cut using cuts for proton          << 
159   G4double cut =                               << 
160     std::max(recoilThreshold, (*pCuts)[Current << 
161                                                << 
162   // find out lightest element                 << 
163   const G4ElementVector* theElementVector = ma << 
164   std::size_t nelm = material->GetNumberOfElem << 
165                                                << 
166   // select lightest element                   << 
167   G4int Z = 300;                               << 
168   for (std::size_t j=0; j<nelm; ++j) {         << 
169     Z = std::min(Z,(*theElementVector)[j]->Get << 
170   }                                               141   }
171   G4int A = G4lrint(fNistManager->GetAtomicMas << 
172   G4double targetMass = G4NucleiProperties::Ge << 
173   G4double t = std::max(cut, 0.5*(cut + std::s << 
174                                                << 
175   return t;                                    << 
176 }                                                 142 }
177                                                   143 
178 //....oooOO0OOooo........oooOO0OOooo........oo    144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179                                                   145 
180 G4double G4hCoulombScatteringModel::ComputeCro    146 G4double G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(
181                 const G4ParticleDefinition* p,    147                 const G4ParticleDefinition* p,
182     G4double kinEnergy,                           148     G4double kinEnergy,
183     G4double Z, G4double,                         149     G4double Z, G4double,
184     G4double cutEnergy, G4double)                 150     G4double cutEnergy, G4double)
185 {                                                 151 {
186   //G4cout << "### G4hCoulombScatteringModel::    152   //G4cout << "### G4hCoulombScatteringModel::ComputeCrossSectionPerAtom  for " 
187   //<< p->GetParticleName()<<" Z= "<<Z<<" e(Me << 153   //  << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl; 
188   G4double cross = 0.0;                           154   G4double cross = 0.0;
189   elecRatio = 0.0;                             << 
190   if(p != particle) { SetupParticle(p); }         155   if(p != particle) { SetupParticle(p); }
191                                                   156 
192   // cross section is set to zero to avoid pro    157   // cross section is set to zero to avoid problems in sample secondary
193   if(kinEnergy <= 0.0) { return cross; }          158   if(kinEnergy <= 0.0) { return cross; }
194   DefineMaterial(CurrentCouple());                159   DefineMaterial(CurrentCouple());
195                                                << 160   cosTetMinNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
196   G4int iz = G4lrint(Z);                       << 161   if(cosThetaMax < cosTetMinNuc) {
197   G4double tmass = (1 == iz) ? proton_mass_c2  << 162     G4int iz = G4int(Z);
198     fNistManager->GetAtomicMassAmu(iz)*amu_c2; << 163     cosTetMinNuc = wokvi->SetupTarget(iz, cutEnergy);
199   wokvi->SetTargetMass(tmass);                 << 164     cosTetMaxNuc = cosThetaMax; 
200                                                << 165     if(iz == 1 && cosTetMaxNuc < 0.0 && particle == theProton) { 
201   G4double costmin =                           << 166       cosTetMaxNuc = 0.0; 
202     wokvi->SetupKinematic(kinEnergy, currentMa << 
203                                                << 
204   if(cosThetaMax < costmin) {                  << 
205     G4double cut = (0.0 < fixedCut) ? fixedCut << 
206     costmin = wokvi->SetupTarget(iz, cut);     << 
207     G4double costmax =                         << 
208       (1 == iz && particle == theProton && cos << 
209       ? 0.0 : cosThetaMax;                     << 
210     if(costmin > costmax) {                    << 
211       cross = wokvi->ComputeNuclearCrossSectio << 
212   + wokvi->ComputeElectronCrossSection(costmin << 
213     }                                             167     }
214     /*                                         << 168     cross =  wokvi->ComputeNuclearCrossSection(cosTetMinNuc, cosTetMaxNuc);
215   if(p->GetParticleName() == "mu+")            << 169     elecRatio = wokvi->ComputeElectronCrossSection(cosTetMinNuc, cosThetaMax);
                                                   >> 170     cross += elecRatio;
                                                   >> 171     if(cross > 0.0) { elecRatio /= cross; }  
                                                   >> 172   }
                                                   >> 173   /*
                                                   >> 174   if(p->GetParticleName() == "e-") 
216   G4cout << "e(MeV)= " << kinEnergy/MeV << " c    175   G4cout << "e(MeV)= " << kinEnergy/MeV << " cross(b)= " << cross/barn  
217    << " 1-costmin= " << 1-costmin              << 176    << " 1-cosTetMinNuc= " << 1-cosTetMinNuc
218    << " 1-costmax= " << 1-costmax              << 177    << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
219    << " 1-cosThetaMax= " << 1-cosThetaMax      << 
220    << G4endl;                                     178    << G4endl;
221     */                                         << 179   */
222   }                                            << 
223   return cross;                                   180   return cross;  
224 }                                                 181 }
225                                                   182 
226 //....oooOO0OOooo........oooOO0OOooo........oo    183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227                                                   184 
228 void G4hCoulombScatteringModel::SampleSecondar    185 void G4hCoulombScatteringModel::SampleSecondaries(
229                 std::vector<G4DynamicParticle*    186                 std::vector<G4DynamicParticle*>* fvect,
230     const G4MaterialCutsCouple* couple,           187     const G4MaterialCutsCouple* couple,
231     const G4DynamicParticle* dp,                  188     const G4DynamicParticle* dp,
232     G4double cutEnergy,                           189     G4double cutEnergy,
233     G4double)                                     190     G4double)
234 {                                                 191 {
235   G4double kinEnergy = dp->GetKineticEnergy();    192   G4double kinEnergy = dp->GetKineticEnergy();
                                                   >> 193 
                                                   >> 194   // absorb particle below low-energy limit to avoid situation
                                                   >> 195   // when a particle has no energy loss
                                                   >> 196   if(kinEnergy < lowEnergyThreshold) { 
                                                   >> 197     fParticleChange->SetProposedKineticEnergy(0.0);
                                                   >> 198     fParticleChange->ProposeLocalEnergyDeposit(kinEnergy);
                                                   >> 199     fParticleChange->ProposeNonIonizingEnergyDeposit(kinEnergy);
                                                   >> 200     return; 
                                                   >> 201   }
236   SetupParticle(dp->GetDefinition());             202   SetupParticle(dp->GetDefinition());
237   DefineMaterial(couple);                         203   DefineMaterial(couple);
238                                                   204 
                                                   >> 205   //G4cout << "G4hCoulombScatteringModel::SampleSecondaries e(MeV)= " 
                                                   >> 206   //   << kinEnergy << "  " << particle->GetParticleName() 
                                                   >> 207   //   << " cut= " << cutEnergy<< G4endl;
                                                   >> 208  
239   // Choose nucleus                               209   // Choose nucleus
240   G4double cut = (0.0 < fixedCut) ? fixedCut : << 210   const G4Element* currentElement = 
                                                   >> 211     SelectRandomAtom(couple,particle,kinEnergy,cutEnergy,kinEnergy);
241                                                   212 
242   const G4Element* elm = SelectRandomAtom(coup << 213   G4double Z = currentElement->GetZ();
243             kinEnergy,cut,kinEnergy);          << 
244                                                   214 
245   G4int iz = elm->GetZasInt();                 << 215   if(ComputeCrossSectionPerAtom(particle,kinEnergy, Z,
246   G4int ia = SelectIsotopeNumber(elm);         << 216           kinEnergy, cutEnergy, kinEnergy) == 0.0) 
247   G4double mass2 = G4NucleiProperties::GetNucl << 217     { return; }
248                                                << 218 
249   wokvi->SetTargetMass(mass2);                 << 219   G4int iz = G4int(Z);
250   wokvi->SetupKinematic(kinEnergy, currentMate << 220   G4int ia = SelectIsotopeNumber(currentElement);
251   G4double costmin = wokvi->SetupTarget(iz, cu << 221   G4double targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
252   G4double costmax = (1 == iz && particle == t << 222   wokvi->SetTargetMass(targetMass);
253     ? 0.0 :  cosThetaMax;                      << 
254   if(costmin <= costmax) { return; }           << 
255                                                << 
256   G4double cross = wokvi->ComputeNuclearCrossS << 
257   G4double ecross = wokvi->ComputeElectronCros << 
258   G4double ratio = ecross/(cross + ecross);    << 
259                                                   223 
260   G4ThreeVector newDirection =                    224   G4ThreeVector newDirection = 
261     wokvi->SampleSingleScattering(costmin, cos << 225     wokvi->SampleSingleScattering(cosTetMinNuc, cosThetaMax, elecRatio);
                                                   >> 226   G4double cost = newDirection.z();
                                                   >> 227 
                                                   >> 228   G4ThreeVector direction = dp->GetMomentumDirection(); 
                                                   >> 229   newDirection.rotateUz(direction);   
262                                                   230 
263   // kinematics in the Lab system              << 
264   G4double ptot = std::sqrt(kinEnergy*(kinEner << 
265   G4double e1   = mass + kinEnergy;            << 
266                                                << 
267   // Lab. system kinematics along projectile d << 
268   G4LorentzVector v0 = G4LorentzVector(0, 0, p << 
269   G4LorentzVector v1 = G4LorentzVector(0, 0, p << 
270   G4ThreeVector bst = v0.boostVector();        << 
271   v1.boost(-bst);                              << 
272   // CM projectile                             << 
273   G4double momCM = v1.pz();                    << 
274                                                << 
275   // Momentum after scattering of incident par << 
276   v1.setX(momCM*newDirection.x());             << 
277   v1.setY(momCM*newDirection.y());             << 
278   v1.setZ(momCM*newDirection.z());             << 
279                                                << 
280   // CM--->Lab                                 << 
281   v1.boost(bst);                               << 
282                                                << 
283   G4ThreeVector dir = dp->GetMomentumDirection << 
284   newDirection = v1.vect().unit();             << 
285   newDirection.rotateUz(dir);                  << 
286                                                << 
287   fParticleChange->ProposeMomentumDirection(ne    231   fParticleChange->ProposeMomentumDirection(newDirection);   
288                                                << 232 
289   // recoil                                    << 233   // recoil sampling assuming a small recoil
290   v0 -= v1;                                    << 234   // and first order correction to primary 4-momentum
291   G4double trec = std::max(v0.e() - mass2, 0.0 << 235   G4double mom2 = wokvi->GetMomentumSquare();
                                                   >> 236   G4double trec = mom2*(1.0 - cost)
                                                   >> 237     /(targetMass + (mass + kinEnergy)*(1.0 - cost));
                                                   >> 238 
                                                   >> 239   // the check likely not needed
                                                   >> 240   if(trec > kinEnergy) { trec = kinEnergy; }
                                                   >> 241   G4double finalT = kinEnergy - trec; 
292   G4double edep = 0.0;                            242   G4double edep = 0.0;
                                                   >> 243   //G4cout<<"G4eCoulombScatteringModel: finalT= "<<finalT<<" Trec= "
                                                   >> 244   //  <<trec << " Z= " << iz << " A= " << ia<<G4endl;
293                                                   245 
294   G4double tcut = recoilThreshold;                246   G4double tcut = recoilThreshold;
295   if(pCuts) { tcut= std::max(tcut,(*pCuts)[cur    247   if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
296                                                << 248 
297   if(trec > tcut) {                               249   if(trec > tcut) {
298     G4ParticleDefinition* ion = theIonTable->G    250     G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
299     newDirection = v0.vect().unit();           << 251     G4ThreeVector dir = (direction*sqrt(mom2) - 
300     newDirection.rotateUz(dir);                << 252        newDirection*sqrt(finalT*(2*mass + finalT))).unit();
301     auto newdp = new G4DynamicParticle(ion, ne << 253     G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec);
302     fvect->push_back(newdp);                      254     fvect->push_back(newdp);
303   } else if(trec > 0.0) {                      << 255   } else {
304     edep = trec;                                  256     edep = trec;
305     fParticleChange->ProposeNonIonizingEnergyD    257     fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
306   }                                               258   }
307                                                   259 
308   // finelize primary energy and energy balanc    260   // finelize primary energy and energy balance
309   G4double finalT = v1.e() - mass;             << 261   // this threshold may be applied only because for low-enegry
310   if(finalT < 0.0) {                           << 262   // e+e- msc model is applied
311     edep += finalT;                            << 263   if(finalT <= lowEnergyThreshold) { 
                                                   >> 264     edep += finalT;  
312     finalT = 0.0;                                 265     finalT = 0.0;
313   }                                               266   } 
314   edep = std::max(edep, 0.0);                  << 
315   fParticleChange->SetProposedKineticEnergy(fi    267   fParticleChange->SetProposedKineticEnergy(finalT);
316   fParticleChange->ProposeLocalEnergyDeposit(e    268   fParticleChange->ProposeLocalEnergyDeposit(edep);
317 }                                                 269 }
318                                                   270 
319 //....oooOO0OOooo........oooOO0OOooo........oo    271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 272 
                                                   >> 273 
320                                                   274