Geant4 Cross Reference

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


  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: G4MollerBhabhaModel.cc,v 1.30 2007/05/22 17:34:36 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-00 $
 26 //                                                 28 //
 27 // -------------------------------------------     29 // -------------------------------------------------------------------
 28 //                                                 30 //
 29 // GEANT4 Class file                               31 // GEANT4 Class file
 30 //                                                 32 //
 31 //                                                 33 //
 32 // File name:   G4MollerBhabhaModel                34 // File name:   G4MollerBhabhaModel
 33 //                                                 35 //
 34 // Author:        Vladimir Ivanchenko on base      36 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 35 //                                                 37 //
 36 // Creation date: 03.01.2002                       38 // Creation date: 03.01.2002
 37 //                                                 39 //
 38 // Modifications:                                  40 // Modifications:
 39 //                                                 41 //
 40 // 13-11-02 Minor fix - use normalised directi     42 // 13-11-02 Minor fix - use normalised direction (V.Ivanchenko)
 41 // 04-12-02 Change G4DynamicParticle construct     43 // 04-12-02 Change G4DynamicParticle constructor in PostStepDoIt (V.Ivanchenko)
 42 // 23-12-02 Change interface in order to move      44 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
 43 // 27-01-03 Make models region aware (V.Ivanch     45 // 27-01-03 Make models region aware (V.Ivanchenko)
 44 // 13-02-03 Add name (V.Ivanchenko)                46 // 13-02-03 Add name (V.Ivanchenko)
 45 // 08-04-05 Major optimisation of internal int     47 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
 46 // 25-07-05 Add protection in calculation of r     48 // 25-07-05 Add protection in calculation of recoil direction for the case 
 47 //          of complete energy transfer from e     49 //          of complete energy transfer from e+ to e- (V.Ivanchenko)
 48 // 06-02-06 ComputeCrossSectionPerElectron, Co     50 // 06-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
 49 // 15-05-06 Fix MinEnergyCut (V.Ivanchenko)        51 // 15-05-06 Fix MinEnergyCut (V.Ivanchenko)
 50 //                                                 52 //
 51 //                                                 53 //
 52 // Class Description:                              54 // Class Description:
 53 //                                                 55 //
 54 // Implementation of energy loss and delta-ele     56 // Implementation of energy loss and delta-electron production by e+/e-
 55 //                                                 57 //
 56 // -------------------------------------------     58 // -------------------------------------------------------------------
 57 //                                                 59 //
 58 //....oooOO0OOooo........oooOO0OOooo........oo     60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 59 //....oooOO0OOooo........oooOO0OOooo........oo     61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 60                                                    62 
 61 #include "G4MollerBhabhaModel.hh"                  63 #include "G4MollerBhabhaModel.hh"
 62 #include "G4PhysicalConstants.hh"              << 
 63 #include "G4SystemOfUnits.hh"                  << 
 64 #include "G4Electron.hh"                           64 #include "G4Electron.hh"
 65 #include "G4Positron.hh"                           65 #include "G4Positron.hh"
 66 #include "Randomize.hh"                            66 #include "Randomize.hh"
 67 #include "G4ParticleChangeForLoss.hh"              67 #include "G4ParticleChangeForLoss.hh"
 68 #include "G4Log.hh"                            << 
 69 #include "G4DeltaAngle.hh"                     << 
 70                                                    68 
 71 //....oooOO0OOooo........oooOO0OOooo........oo     69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 72                                                    70 
 73 using namespace std;                               71 using namespace std;
 74                                                    72 
 75 G4MollerBhabhaModel::G4MollerBhabhaModel(const     73 G4MollerBhabhaModel::G4MollerBhabhaModel(const G4ParticleDefinition* p,
 76                                          const     74                                          const G4String& nam)
 77   : G4VEmModel(nam),                               75   : G4VEmModel(nam),
 78     particle(nullptr),                         <<  76   particle(0),
 79     isElectron(true),                          <<  77   isElectron(true),
 80     twoln10(2.0*G4Log(10.0)),                  <<  78   twoln10(2.0*log(10.0)),
 81     lowLimit(0.02*keV),                        <<  79   lowLimit(0.2*keV)
 82     isInitialised(false)                       << 
 83 {                                                  80 {
 84   theElectron = G4Electron::Electron();            81   theElectron = G4Electron::Electron();
 85   if(nullptr != p) { SetParticle(p); }         <<  82   if(p) SetParticle(p);
 86   fParticleChange = nullptr;                   << 
 87 }                                                  83 }
 88                                                    84 
 89 //....oooOO0OOooo........oooOO0OOooo........oo     85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 90                                                    86 
 91 G4MollerBhabhaModel::~G4MollerBhabhaModel() =  <<  87 G4MollerBhabhaModel::~G4MollerBhabhaModel()
                                                   >>  88 {}
 92                                                    89 
 93 //....oooOO0OOooo........oooOO0OOooo........oo     90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 94                                                    91 
 95 G4double G4MollerBhabhaModel::MaxSecondaryEner <<  92 void G4MollerBhabhaModel::SetParticle(const G4ParticleDefinition* p)
 96                                                << 
 97 {                                                  93 {
 98   G4double tmax = kinEnergy;                   <<  94   particle = p;
 99   if(isElectron) { tmax *= 0.5; }              <<  95   if(p != theElectron) isElectron = false;
100   return tmax;                                 << 
101 }                                                  96 }
102                                                    97 
103 //....oooOO0OOooo........oooOO0OOooo........oo     98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104                                                    99 
105 void G4MollerBhabhaModel::Initialise(const G4P << 100 G4double G4MollerBhabhaModel::MinEnergyCut(const G4ParticleDefinition*,
106                                      const G4D << 101                                            const G4MaterialCutsCouple* couple)
107 {                                                 102 {
108   if(p != particle) { SetParticle(p); }        << 103   G4double electronDensity = couple->GetMaterial()->GetElectronDensity();
                                                   >> 104   G4double Zeff  = electronDensity/couple->GetMaterial()->GetTotNbOfAtomsPerVolume();
                                                   >> 105   return 0.25*sqrt(Zeff)*keV;
                                                   >> 106 }
109                                                   107 
110   if(isInitialised) { return; }                << 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111                                                   109 
112   isInitialised = true;                        << 110 void G4MollerBhabhaModel::Initialise(const G4ParticleDefinition* p,
113   fParticleChange = GetParticleChangeForLoss() << 111                                      const G4DataVector&)
114   if(UseAngularGeneratorFlag() && !GetAngularD << 112 {
115     SetAngularDistribution(new G4DeltaAngle()) << 113   if(!particle) SetParticle(p);
116   }                                            << 114   if(pParticleChange)
                                                   >> 115     fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>
                                                   >> 116                                                      (pParticleChange);
                                                   >> 117   else
                                                   >> 118     fParticleChange = new G4ParticleChangeForLoss();
117 }                                                 119 }
118                                                   120 
119 //....oooOO0OOooo........oooOO0OOooo........oo    121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120                                                   122 
121 G4double G4MollerBhabhaModel::ComputeCrossSect    123 G4double G4MollerBhabhaModel::ComputeCrossSectionPerElectron(
122          const G4ParticleDefinition* p, G4doub << 124                                            const G4ParticleDefinition* p,
123    G4double cutEnergy, G4double maxEnergy)     << 125                                                  G4double kineticEnergy,
                                                   >> 126                                                  G4double cutEnergy,
                                                   >> 127                                                  G4double maxEnergy)
124 {                                                 128 {
125   if(p != particle) { SetParticle(p); }        << 129   if(!particle) SetParticle(p);
126                                                   130 
127   G4double cross = 0.0;                           131   G4double cross = 0.0;
128   G4double tmax = MaxSecondaryEnergy(p, kineti    132   G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
129   tmax = std::min(maxEnergy, tmax);            << 133   tmax = min(maxEnergy, tmax);
130   //G4cout << "E= " << kineticEnergy << " cut= << 134 
131   //         << " Emax= " << tmax << G4endl;   << 
132   if(cutEnergy < tmax) {                          135   if(cutEnergy < tmax) {
133                                                   136 
134     G4double xmin  = cutEnergy/kineticEnergy;     137     G4double xmin  = cutEnergy/kineticEnergy;
135     G4double xmax  = tmax/kineticEnergy;          138     G4double xmax  = tmax/kineticEnergy;
136     G4double tau   = kineticEnergy/electron_ma << 139     G4double gam   = kineticEnergy/electron_mass_c2 + 1.0;
137     G4double gam   = tau + 1.0;                << 
138     G4double gamma2= gam*gam;                     140     G4double gamma2= gam*gam;
139     G4double beta2 = tau*(tau + 2)/gamma2;     << 141     G4double beta2 = 1.0 - 1.0/gamma2;
140                                                   142 
141     //Moller (e-e-) scattering                    143     //Moller (e-e-) scattering
142     if (isElectron) {                             144     if (isElectron) {
143                                                   145 
144       G4double gg = (2.0*gam - 1.0)/gamma2;    << 146       G4double g     = (2.0*gam - 1.0)/gamma2;
145       cross = ((xmax - xmin)*(1.0 - gg + 1.0/( << 147       cross = ((xmax - xmin)*(1.0 - g + 1.0/(xmin*xmax)
146                               + 1.0/((1.0-xmin << 148             + 1.0/((1.0-xmin)*(1.0 - xmax)))
147             - gg*G4Log( xmax*(1.0 - xmin)/(xmi << 149             - g*log( xmax*(1.0 - xmin)/(xmin*(1.0 - xmax)) ) ) / beta2;
148                                                   150 
149     //Bhabha (e+e-) scattering                    151     //Bhabha (e+e-) scattering
150     } else {                                      152     } else {
151                                                   153 
152       G4double y   = 1.0/(1.0 + gam);             154       G4double y   = 1.0/(1.0 + gam);
153       G4double y2  = y*y;                         155       G4double y2  = y*y;
154       G4double y12 = 1.0 - 2.0*y;                 156       G4double y12 = 1.0 - 2.0*y;
155       G4double b1  = 2.0 - y2;                    157       G4double b1  = 2.0 - y2;
156       G4double b2  = y12*(3.0 + y2);              158       G4double b2  = y12*(3.0 + y2);
157       G4double y122= y12*y12;                     159       G4double y122= y12*y12;
158       G4double b4  = y122*y12;                    160       G4double b4  = y122*y12;
159       G4double b3  = b4 + y122;                   161       G4double b3  = b4 + y122;
160                                                   162 
161       cross = (xmax - xmin)*(1.0/(beta2*xmin*x    163       cross = (xmax - xmin)*(1.0/(beta2*xmin*xmax) + b2
162             - 0.5*b3*(xmin + xmax)                164             - 0.5*b3*(xmin + xmax)
163             + b4*(xmin*xmin + xmin*xmax + xmax << 165       + b4*(xmin*xmin + xmin*xmax + xmax*xmax)/3.0)
164             - b1*G4Log(xmax/xmin);             << 166             - b1*log(xmax/xmin);
165     }                                             167     }
166                                                   168 
167     cross *= twopi_mc2_rcl2/kineticEnergy;        169     cross *= twopi_mc2_rcl2/kineticEnergy;
168   }                                               170   }
169   return cross;                                   171   return cross;
170 }                                                 172 }
171                                                   173 
172 //....oooOO0OOooo........oooOO0OOooo........oo    174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173                                                   175 
174 G4double G4MollerBhabhaModel::ComputeCrossSect    176 G4double G4MollerBhabhaModel::ComputeCrossSectionPerAtom(
175                                            con    177                                            const G4ParticleDefinition* p,
176                                                   178                                                  G4double kineticEnergy,
177                                                << 179              G4double Z, G4double,
178                                                   180                                                  G4double cutEnergy,
179                                                   181                                                  G4double maxEnergy)
180 {                                                 182 {
181   return Z*ComputeCrossSectionPerElectron(p,ki << 183   G4double cross = Z*ComputeCrossSectionPerElectron
                                                   >> 184                                          (p,kineticEnergy,cutEnergy,maxEnergy);
                                                   >> 185   return cross;
182 }                                                 186 }
183                                                   187 
184 //....oooOO0OOooo........oooOO0OOooo........oo    188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185                                                   189 
186 G4double G4MollerBhabhaModel::CrossSectionPerV    190 G4double G4MollerBhabhaModel::CrossSectionPerVolume(
187                                            con << 191              const G4Material* material,
188                                            con    192                                            const G4ParticleDefinition* p,
189                                                << 193                                                  G4double kineticEnergy,
190                                                   194                                                  G4double cutEnergy,
191                                                   195                                                  G4double maxEnergy)
192 {                                                 196 {
193   G4double eDensity = material->GetElectronDen    197   G4double eDensity = material->GetElectronDensity();
194   return eDensity*ComputeCrossSectionPerElectr << 198   G4double cross = eDensity*ComputeCrossSectionPerElectron
                                                   >> 199                                          (p,kineticEnergy,cutEnergy,maxEnergy);
                                                   >> 200   return cross;
195 }                                                 201 }
196                                                   202 
197 //....oooOO0OOooo........oooOO0OOooo........oo    203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
198                                                   204 
199 G4double G4MollerBhabhaModel::ComputeDEDXPerVo    205 G4double G4MollerBhabhaModel::ComputeDEDXPerVolume(
200                                           cons << 206             const G4Material* material,
201                                           cons    207                                           const G4ParticleDefinition* p,
202                                                   208                                                 G4double kineticEnergy,
203                                                << 209                                                 G4double cutEnergy)
204 {                                                 210 {
205   if(p != particle) { SetParticle(p); }        << 211   if(!particle) SetParticle(p);
206   // calculate the dE/dx due to the ionization    212   // calculate the dE/dx due to the ionization by Seltzer-Berger formula
207   // checl low-energy limit                    << 
208   G4double electronDensity = material->GetElec << 
209                                                   213   
210   G4double Zeff  = material->GetIonisation()-> << 214   G4double electronDensity = material->GetElectronDensity();
                                                   >> 215   G4double Zeff  = electronDensity/material->GetTotNbOfAtomsPerVolume();
211   G4double th    = 0.25*sqrt(Zeff)*keV;           216   G4double th    = 0.25*sqrt(Zeff)*keV;
212   G4double tkin = std::max(kineticEnergy, th); << 217   G4double tkin  = kineticEnergy;
213                                                << 218   G4bool   lowEnergy = false;
                                                   >> 219   if (kineticEnergy < th) {
                                                   >> 220     tkin = th;
                                                   >> 221     lowEnergy = true;
                                                   >> 222   }
214   G4double tau   = tkin/electron_mass_c2;         223   G4double tau   = tkin/electron_mass_c2;
215   G4double gam   = tau + 1.0;                     224   G4double gam   = tau + 1.0;
216   G4double gamma2= gam*gam;                       225   G4double gamma2= gam*gam;
217   G4double bg2   = tau*(tau + 2);              << 226   G4double beta2 = 1. - 1./gamma2;
218   G4double beta2 = bg2/gamma2;                 << 227   G4double bg2   = beta2*gamma2;
219                                                   228 
220   G4double eexc  = material->GetIonisation()->    229   G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
221   eexc          /= electron_mass_c2;              230   eexc          /= electron_mass_c2;
222   G4double eexc2 = eexc*eexc;                     231   G4double eexc2 = eexc*eexc; 
223                                                   232 
224   G4double d = std::min(cut, MaxSecondaryEnerg << 233   G4double d = min(cutEnergy, MaxSecondaryEnergy(p, tkin))/electron_mass_c2;
225   G4double dedx;                                  234   G4double dedx;
226                                                   235 
227   // electron                                     236   // electron
228   if (isElectron) {                               237   if (isElectron) {
229                                                   238 
230     dedx = G4Log(2.0*(tau + 2.0)/eexc2) - 1.0  << 239     dedx = log(2.0*(tau + 2.0)/eexc2) - 1.0 - beta2
231          + G4Log((tau-d)*d) + tau/(tau-d)      << 240          + log((tau-d)*d) + tau/(tau-d)
232          + (0.5*d*d + (2.0*tau + 1.)*G4Log(1.  << 241          + (0.5*d*d + (2.0*tau + 1.)*log(1. - d/tau))/gamma2;
233                                                   242    
234   //positron                                      243   //positron
235   } else {                                        244   } else {
236                                                   245 
237     G4double d2 = d*d*0.5;                        246     G4double d2 = d*d*0.5;
238     G4double d3 = d2*d/1.5;                       247     G4double d3 = d2*d/1.5;
239     G4double d4 = d3*d*0.75;                   << 248     G4double d4 = d3*d*3.75;
240     G4double y  = 1.0/(1.0 + gam);                249     G4double y  = 1.0/(1.0 + gam);
241     dedx = G4Log(2.0*(tau + 2.0)/eexc2) + G4Lo << 250     dedx = log(2.0*(tau + 2.0)/eexc2) + log(tau*d)
242          - beta2*(tau + 2.0*d - y*(3.0*d2         251          - beta2*(tau + 2.0*d - y*(3.0*d2 
243          + y*(d - d3 + y*(d2 - tau*d3 + d4))))    252          + y*(d - d3 + y*(d2 - tau*d3 + d4))))/tau;
244   }                                               253   } 
245                                                   254 
246   //density correction                            255   //density correction 
247   G4double x = G4Log(bg2)/twoln10;             << 256   G4double cden  = material->GetIonisation()->GetCdensity();
248   dedx -= material->GetIonisation()->DensityCo << 257   G4double mden  = material->GetIonisation()->GetMdensity();
                                                   >> 258   G4double aden  = material->GetIonisation()->GetAdensity();
                                                   >> 259   G4double x0den = material->GetIonisation()->GetX0density();
                                                   >> 260   G4double x1den = material->GetIonisation()->GetX1density();
                                                   >> 261   G4double x     = log(bg2)/twoln10;
                                                   >> 262 
                                                   >> 263   if (x >= x0den) {
                                                   >> 264     dedx -= twoln10*x - cden;
                                                   >> 265     if (x < x1den) dedx -= aden*pow(x1den-x, mden);
                                                   >> 266   } 
249                                                   267 
250   // now you can compute the total ionization     268   // now you can compute the total ionization loss
251   dedx *= twopi_mc2_rcl2*electronDensity/beta2    269   dedx *= twopi_mc2_rcl2*electronDensity/beta2;
252   if (dedx < 0.0) { dedx = 0.0; }              << 270   if (dedx < 0.0) dedx = 0.0;
253                                                << 271 
254   // lowenergy extrapolation                      272   // lowenergy extrapolation
255                                                << 273 
256   if (kineticEnergy < th) {                    << 274   if (lowEnergy) {
257     x = kineticEnergy/th;                      << 275 
258     if(x > 0.25) { dedx /= sqrt(x); }          << 276     if (kineticEnergy >= lowLimit) dedx *= sqrt(tkin/kineticEnergy);
259     else { dedx *= 1.4*sqrt(x)/(0.1 + x); }    << 277     else                           dedx *= sqrt(tkin*kineticEnergy)/lowLimit;
                                                   >> 278 
260   }                                               279   }
261   return dedx;                                    280   return dedx;
262 }                                                 281 }
263                                                   282 
264 //....oooOO0OOooo........oooOO0OOooo........oo    283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
265                                                   284 
266 void                                           << 285 void G4MollerBhabhaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
267 G4MollerBhabhaModel::SampleSecondaries(std::ve << 286               const G4MaterialCutsCouple*,
268                                        const G << 287               const G4DynamicParticle* dp,
269                                        const G << 288               G4double tmin,
270                                        G4doubl << 289               G4double maxEnergy)
271                                        G4doubl << 
272 {                                                 290 {
273   G4double kineticEnergy = dp->GetKineticEnerg << 291   G4double tmax = std::min(maxEnergy, MaxSecondaryKinEnergy(dp));
274   //G4cout << "G4MollerBhabhaModel::SampleSeco << 292   if(tmin >= tmax) return;
275   //       << " in " << couple->GetMaterial()- << 
276   G4double tmax;                               << 
277   G4double tmin = cutEnergy;                   << 
278   if(isElectron) {                             << 
279     tmax = 0.5*kineticEnergy;                  << 
280   } else {                                     << 
281     tmax = kineticEnergy;                      << 
282   }                                            << 
283   if(maxEnergy < tmax) { tmax = maxEnergy; }   << 
284   if(tmin >= tmax) { return; }                 << 
285                                                   293 
                                                   >> 294   G4double kineticEnergy = dp->GetKineticEnergy();
286   G4double energy = kineticEnergy + electron_m    295   G4double energy = kineticEnergy + electron_mass_c2;
                                                   >> 296   G4double totalMomentum = sqrt(kineticEnergy*(energy + electron_mass_c2));
287   G4double xmin   = tmin/kineticEnergy;           297   G4double xmin   = tmin/kineticEnergy;
288   G4double xmax   = tmax/kineticEnergy;           298   G4double xmax   = tmax/kineticEnergy;
289   G4double gam    = energy/electron_mass_c2;      299   G4double gam    = energy/electron_mass_c2;
290   G4double gamma2 = gam*gam;                      300   G4double gamma2 = gam*gam;
291   G4double beta2  = 1.0 - 1.0/gamma2;             301   G4double beta2  = 1.0 - 1.0/gamma2;
292   G4double x, z, grej;                         << 302   G4double x, z, q, grej;
293   CLHEP::HepRandomEngine* rndmEngine = G4Rando << 303 
294   G4double rndm[2];                            << 304   G4ThreeVector direction = dp->GetMomentumDirection();
295                                                   305 
296   //Moller (e-e-) scattering                      306   //Moller (e-e-) scattering
297   if (isElectron) {                               307   if (isElectron) {
298                                                   308 
299     G4double gg = (2.0*gam - 1.0)/gamma2;      << 309     G4double g = (2.0*gam - 1.0)/gamma2;
300     G4double y = 1.0 - xmax;                      310     G4double y = 1.0 - xmax;
301     grej = 1.0 - gg*xmax + xmax*xmax*(1.0 - gg << 311     grej = 1.0 - g*xmax + xmax*xmax*(1.0 - g + (1.0 - g*y)/(y*y));
302                                                   312 
303     do {                                          313     do {
304       rndmEngine->flatArray(2, rndm);          << 314       q = G4UniformRand();
305       x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xm << 315       x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
306       y = 1.0 - x;                                316       y = 1.0 - x;
307       z = 1.0 - gg*x + x*x*(1.0 - gg + (1.0 -  << 317       z = 1.0 - g*x + x*x*(1.0 - g + (1.0 - g*y)/(y*y));
308       /*                                          318       /*
309       if(z > grej) {                              319       if(z > grej) {
310         G4cout << "G4MollerBhabhaModel::Sample    320         G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
311                << "Majorant " << grej << " < "    321                << "Majorant " << grej << " < "
312                << z << " for x= " << x            322                << z << " for x= " << x
313                << " e-e- scattering"              323                << " e-e- scattering"
314                << G4endl;                         324                << G4endl;
315       }                                           325       }
316       */                                          326       */
317       // Loop checking, 03-Aug-2015, Vladimir  << 327     } while(grej * G4UniformRand() > z);
318     } while(grej * rndm[1] > z);               << 
319                                                   328 
320   //Bhabha (e+e-) scattering                      329   //Bhabha (e+e-) scattering
321   } else {                                        330   } else {
322                                                   331 
323     G4double y   = 1.0/(1.0 + gam);               332     G4double y   = 1.0/(1.0 + gam);
324     G4double y2  = y*y;                           333     G4double y2  = y*y;
325     G4double y12 = 1.0 - 2.0*y;                   334     G4double y12 = 1.0 - 2.0*y;
326     G4double b1  = 2.0 - y2;                      335     G4double b1  = 2.0 - y2;
327     G4double b2  = y12*(3.0 + y2);                336     G4double b2  = y12*(3.0 + y2);
328     G4double y122= y12*y12;                       337     G4double y122= y12*y12;
329     G4double b4  = y122*y12;                      338     G4double b4  = y122*y12;
330     G4double b3  = b4 + y122;                     339     G4double b3  = b4 + y122;
331                                                   340 
332     y    = xmax*xmax;                          << 341     y     = xmax*xmax;
333     grej = 1.0 + (y*y*b4 - xmin*xmin*xmin*b3 + << 342     grej  = -xmin*b1;
                                                   >> 343     grej += y*b2;
                                                   >> 344     grej -= xmin*xmin*xmin*b3;
                                                   >> 345     grej += y*y*b4;
                                                   >> 346     grej *= beta2;
                                                   >> 347     grej += 1.0;
334     do {                                          348     do {
335       rndmEngine->flatArray(2, rndm);          << 349       q  = G4UniformRand();
336       x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xm << 350       x  = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
337       y = x*x;                                 << 351       z  = -x*b1;
338       z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1 << 352       y  = x*x;
                                                   >> 353       z += y*b2;
                                                   >> 354       y *= x;
                                                   >> 355       z -= y*b3;
                                                   >> 356       y *= x;
                                                   >> 357       z += y*b4;
                                                   >> 358       z *= beta2;
                                                   >> 359       z += 1.0;
339       /*                                          360       /*
340       if(z > grej) {                              361       if(z > grej) {
341         G4cout << "G4MollerBhabhaModel::Sample    362         G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
342                << "Majorant " << grej << " < "    363                << "Majorant " << grej << " < "
343                << z << " for x= " << x            364                << z << " for x= " << x
344                << " e+e- scattering"              365                << " e+e- scattering"
345                << G4endl;                         366                << G4endl;
346       }                                           367       }
347       */                                          368       */
348       // Loop checking, 03-Aug-2015, Vladimir  << 369     } while(grej * G4UniformRand() > z);
349     } while(grej * rndm[1] > z);               << 
350   }                                               370   }
351                                                   371 
352   G4double deltaKinEnergy = x * kineticEnergy;    372   G4double deltaKinEnergy = x * kineticEnergy;
353                                                   373 
354   G4ThreeVector deltaDirection;                << 374   G4double deltaMomentum =
                                                   >> 375            sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
                                                   >> 376   G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
                                                   >> 377                                    (deltaMomentum * totalMomentum);
                                                   >> 378   G4double sint = 1.0 - cost*cost;
                                                   >> 379   if(sint > 0.0) sint = sqrt(sint);
355                                                   380 
356   if(UseAngularGeneratorFlag()) {              << 381   G4double phi = twopi * G4UniformRand() ;
357     const G4Material* mat =  couple->GetMateri << 
358     G4int Z = SelectRandomAtomNumber(mat);     << 
359                                                   382 
360     deltaDirection =                           << 383   G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
361       GetAngularDistribution()->SampleDirectio << 384   deltaDirection.rotateUz(direction);
362                                                << 
363   } else {                                     << 
364                                                << 
365     G4double deltaMomentum =                   << 
366       sqrt(deltaKinEnergy * (deltaKinEnergy +  << 
367     G4double cost = deltaKinEnergy * (energy + << 
368       (deltaMomentum * dp->GetTotalMomentum()) << 
369     if(cost > 1.0) { cost = 1.0; }             << 
370     G4double sint = sqrt((1.0 - cost)*(1.0 + c << 
371                                                << 
372     G4double phi = twopi * rndmEngine->flat()  << 
373                                                << 
374     deltaDirection.set(sint*cos(phi),sint*sin( << 
375     deltaDirection.rotateUz(dp->GetMomentumDir << 
376   }                                            << 
377                                                << 
378   // create G4DynamicParticle object for delta << 
379   auto delta = new G4DynamicParticle(theElectr << 
380   vdp->push_back(delta);                       << 
381                                                   385 
382   // primary change                               386   // primary change
383   kineticEnergy -= deltaKinEnergy;                387   kineticEnergy -= deltaKinEnergy;
384   G4ThreeVector finalP = dp->GetMomentum() - d << 
385   finalP               = finalP.unit();        << 
386                                                << 
387   fParticleChange->SetProposedKineticEnergy(ki    388   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
388   fParticleChange->SetProposedMomentumDirectio << 389 
                                                   >> 390   if(kineticEnergy > DBL_MIN) {
                                                   >> 391     G4ThreeVector dir = totalMomentum*direction - deltaMomentum*deltaDirection;
                                                   >> 392     direction = dir.unit();
                                                   >> 393     fParticleChange->SetProposedMomentumDirection(direction);
                                                   >> 394   }
                                                   >> 395 
                                                   >> 396   // create G4DynamicParticle object for delta ray
                                                   >> 397   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
                                                   >> 398              deltaDirection,deltaKinEnergy);
                                                   >> 399   vdp->push_back(delta);
389 }                                                 400 }
390                                                   401 
391 //....oooOO0OOooo........oooOO0OOooo........oo    402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
392                                                   403