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


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