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 8.0)


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