Geant4 Cross Reference

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


  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: G4KleinNishinaModel.cc 82754 2014-07-08 14:06:13Z gcosmo $
 26 //                                                 27 //
 27 // -------------------------------------------     28 // -------------------------------------------------------------------
 28 //                                                 29 //
 29 // GEANT4 Class file                               30 // GEANT4 Class file
 30 //                                                 31 //
 31 //                                                 32 //
 32 // File name:     G4KleinNishinaModel              33 // File name:     G4KleinNishinaModel
 33 //                                                 34 //
 34 // Author:        Vladimir Ivanchenko on base      35 // Author:        Vladimir Ivanchenko on base of G4KleinNishinaCompton
 35 //                                                 36 //
 36 // Creation date: 13.06.2010                       37 // Creation date: 13.06.2010
 37 //                                                 38 //
 38 // Modifications:                                  39 // Modifications:
 39 //                                                 40 //
 40 // Class Description:                              41 // Class Description:
 41 //                                                 42 //
 42 // -------------------------------------------     43 // -------------------------------------------------------------------
 43 //                                                 44 //
 44 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45 //....oooOO0OOooo........oooOO0OOooo........oo     46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 46                                                    47 
 47 #include "G4KleinNishinaModel.hh"                  48 #include "G4KleinNishinaModel.hh"
 48 #include "G4PhysicalConstants.hh"                  49 #include "G4PhysicalConstants.hh"
 49 #include "G4SystemOfUnits.hh"                      50 #include "G4SystemOfUnits.hh"
 50 #include "G4Electron.hh"                           51 #include "G4Electron.hh"
 51 #include "G4Gamma.hh"                              52 #include "G4Gamma.hh"
 52 #include "Randomize.hh"                            53 #include "Randomize.hh"
 53 #include "G4RandomDirection.hh"                    54 #include "G4RandomDirection.hh"
 54 #include "G4DataVector.hh"                         55 #include "G4DataVector.hh"
 55 #include "G4ParticleChangeForGamma.hh"             56 #include "G4ParticleChangeForGamma.hh"
 56 #include "G4VAtomDeexcitation.hh"                  57 #include "G4VAtomDeexcitation.hh"
 57 #include "G4AtomicShells.hh"                       58 #include "G4AtomicShells.hh"
 58 #include "G4LossTableManager.hh"                   59 #include "G4LossTableManager.hh"
 59 #include "G4Log.hh"                                60 #include "G4Log.hh"
 60 #include "G4Exp.hh"                                61 #include "G4Exp.hh"
 61                                                    62 
 62 //....oooOO0OOooo........oooOO0OOooo........oo     63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 63                                                    64 
 64 using namespace std;                               65 using namespace std;
 65                                                    66 
                                                   >>  67 static const G4double
                                                   >>  68   d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLHEP::barn, 
                                                   >>  69   d3= 6.7527   *CLHEP::barn, d4=-1.9798e+1*CLHEP::barn,
                                                   >>  70   e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLHEP::barn, 
                                                   >>  71   e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLHEP::barn,
                                                   >>  72   f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLHEP::barn, 
                                                   >>  73   f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLHEP::barn;
                                                   >>  74 static const G4double dT0 = keV;
                                                   >>  75 static const G4int nlooplim = 1000;
                                                   >>  76 
 66 G4KleinNishinaModel::G4KleinNishinaModel(const     77 G4KleinNishinaModel::G4KleinNishinaModel(const G4String& nam)
 67   : G4VEmModel(nam),                           <<  78   : G4VEmModel(nam)
 68     lv1(0.,0.,0.,0.),                          << 
 69     lv2(0.,0.,0.,0.),                          << 
 70     bst(0.,0.,0.)                              << 
 71 {                                                  79 {
 72   theGamma = G4Gamma::Gamma();                     80   theGamma = G4Gamma::Gamma();
 73   theElectron = G4Electron::Electron();            81   theElectron = G4Electron::Electron();
 74   lowestSecondaryEnergy = 10*eV;                   82   lowestSecondaryEnergy = 10*eV;
 75   limitFactor       = 4;                           83   limitFactor       = 4;
 76   fProbabilities.resize(9,0.0);                    84   fProbabilities.resize(9,0.0);
 77   SetDeexcitationFlag(true);                       85   SetDeexcitationFlag(true);
 78   fParticleChange = nullptr;                   <<  86   fParticleChange = 0;
 79   fAtomDeexcitation = nullptr;                 <<  87   fAtomDeexcitation = 0;
 80 }                                                  88 }
 81                                                    89 
 82 //....oooOO0OOooo........oooOO0OOooo........oo     90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 83                                                    91 
 84 G4KleinNishinaModel::~G4KleinNishinaModel() =  <<  92 G4KleinNishinaModel::~G4KleinNishinaModel()
                                                   >>  93 {}
 85                                                    94 
 86 //....oooOO0OOooo........oooOO0OOooo........oo     95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 87                                                    96 
 88 void G4KleinNishinaModel::Initialise(const G4P     97 void G4KleinNishinaModel::Initialise(const G4ParticleDefinition* p,
 89                                      const G4D <<  98              const G4DataVector& cuts)
 90 {                                                  99 {
 91   fAtomDeexcitation = G4LossTableManager::Inst    100   fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
 92   if(IsMaster()) { InitialiseElementSelectors(    101   if(IsMaster()) { InitialiseElementSelectors(p, cuts); }
 93   if(nullptr == fParticleChange) {             << 102   if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
 94     fParticleChange = GetParticleChangeForGamm << 
 95   }                                            << 
 96 }                                                 103 }
 97                                                   104 
 98 //....oooOO0OOooo........oooOO0OOooo........oo    105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 99                                                   106 
100 void G4KleinNishinaModel::InitialiseLocal(cons    107 void G4KleinNishinaModel::InitialiseLocal(const G4ParticleDefinition*,
101                                           G4VE << 108             G4VEmModel* masterModel)
102 {                                                 109 {
103   SetElementSelectors(masterModel->GetElementS    110   SetElementSelectors(masterModel->GetElementSelectors());
104 }                                                 111 }
105                                                   112 
106 //....oooOO0OOooo........oooOO0OOooo........oo    113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107                                                   114 
108 G4double                                          115 G4double 
109 G4KleinNishinaModel::ComputeCrossSectionPerAto    116 G4KleinNishinaModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
110                                                << 117             G4double gammaEnergy,
111                                                << 118             G4double Z, G4double,
112                                                << 119             G4double, G4double)
113 {                                                 120 {
114   G4double xSection = 0.0 ;                       121   G4double xSection = 0.0 ;
115   if (gammaEnergy <= LowEnergyLimit()) { retur    122   if (gammaEnergy <= LowEnergyLimit()) { return xSection; }
116                                                   123 
117   static const G4double a = 20.0 , b = 230.0 ,    124   static const G4double a = 20.0 , b = 230.0 , c = 440.0;
118                                                << 
119 static const G4double                          << 
120   d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLH << 
121   d3= 6.7527   *CLHEP::barn, d4=-1.9798e+1*CLH << 
122   e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLH << 
123   e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLH << 
124   f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLH << 
125   f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLH << 
126                                                   125   
127   G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z =    126   G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
128            p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z =    127            p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
129                                                   128 
130   G4double T0  = 15.0*keV;                        129   G4double T0  = 15.0*keV; 
131   if (Z < 1.5) { T0 = 40.0*keV; }                 130   if (Z < 1.5) { T0 = 40.0*keV; } 
132                                                   131 
133   G4double X   = max(gammaEnergy, T0) / electr    132   G4double X   = max(gammaEnergy, T0) / electron_mass_c2;
134   xSection = p1Z*G4Log(1.+2.*X)/X                 133   xSection = p1Z*G4Log(1.+2.*X)/X
135                + (p2Z + p3Z*X + p4Z*X*X)/(1. +    134                + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
136                                                << 135     
137   //  modification for low energy. (special ca    136   //  modification for low energy. (special case for Hydrogen)
138   static const G4double dT0 = keV;             << 
139   if (gammaEnergy < T0) {                         137   if (gammaEnergy < T0) {
140     X = (T0+dT0) / electron_mass_c2 ;             138     X = (T0+dT0) / electron_mass_c2 ;
141     G4double sigma = p1Z*G4Log(1.+2*X)/X          139     G4double sigma = p1Z*G4Log(1.+2*X)/X
142                     + (p2Z + p3Z*X + p4Z*X*X)/    140                     + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
143     G4double   c1 = -T0*(sigma-xSection)/(xSec    141     G4double   c1 = -T0*(sigma-xSection)/(xSection*dT0);             
144     G4double   c2 = 0.150;                        142     G4double   c2 = 0.150; 
145     if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z);    143     if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); }
146     G4double    y = G4Log(gammaEnergy/T0);        144     G4double    y = G4Log(gammaEnergy/T0);
147     xSection *= G4Exp(-y*(c1+c2*y));              145     xSection *= G4Exp(-y*(c1+c2*y));          
148   }                                               146   }
149                                                   147 
150   if(xSection < 0.0) { xSection = 0.0; }          148   if(xSection < 0.0) { xSection = 0.0; }
151   //  G4cout << "e= " << GammaEnergy << " Z= "    149   //  G4cout << "e= " << GammaEnergy << " Z= " << Z 
152   //  << " cross= " << xSection << G4endl;        150   //  << " cross= " << xSection << G4endl;
153   return xSection;                                151   return xSection;
154 }                                                 152 }
155                                                   153 
156 //....oooOO0OOooo........oooOO0OOooo........oo    154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157                                                   155 
158 void G4KleinNishinaModel::SampleSecondaries(      156 void G4KleinNishinaModel::SampleSecondaries(
159                              std::vector<G4Dyn << 157            std::vector<G4DynamicParticle*>* fvect,
160                              const G4MaterialC << 158            const G4MaterialCutsCouple* couple,
161                              const G4DynamicPa << 159            const G4DynamicParticle* aDynamicGamma,
162                              G4double,         << 160            G4double,
163                              G4double)         << 161            G4double)
164 {                                                 162 {
165   // primary gamma                                163   // primary gamma
166   G4double energy = aDynamicGamma->GetKineticE    164   G4double energy = aDynamicGamma->GetKineticEnergy();
167                                                   165 
168   // do nothing below the threshold               166   // do nothing below the threshold
169   if(energy <= LowEnergyLimit()) { return; }      167   if(energy <= LowEnergyLimit()) { return; }
170                                                   168 
171   G4ThreeVector direction = aDynamicGamma->Get    169   G4ThreeVector direction = aDynamicGamma->GetMomentumDirection();
172                                                   170 
173   // select atom                                  171   // select atom
174   const G4Element* elm = SelectRandomAtom(coup    172   const G4Element* elm = SelectRandomAtom(couple, theGamma, energy);
175                                                   173 
176   // select shell first                           174   // select shell first
177   G4int nShells = elm->GetNbOfAtomicShells();     175   G4int nShells = elm->GetNbOfAtomicShells();
178   if(nShells > (G4int)fProbabilities.size()) {    176   if(nShells > (G4int)fProbabilities.size()) { fProbabilities.resize(nShells); }
179   G4double totprob = 0.0;                         177   G4double totprob = 0.0;
180   G4int i;                                        178   G4int i;
181   for(i=0; i<nShells; ++i) {                      179   for(i=0; i<nShells; ++i) {
182     //G4double bindingEnergy = elm->GetAtomicS    180     //G4double bindingEnergy = elm->GetAtomicShell(i);
183     totprob += elm->GetNbOfShellElectrons(i);     181     totprob += elm->GetNbOfShellElectrons(i);
184     //totprob += elm->GetNbOfShellElectrons(i)    182     //totprob += elm->GetNbOfShellElectrons(i)/(bindingEnergy*bindingEnergy);
185     fProbabilities[i] = totprob;                  183     fProbabilities[i] = totprob; 
186   }                                               184   }
187                                                   185 
188   // Loop on sampling                             186   // Loop on sampling
189   static const G4int nlooplim = 1000;          << 
190   G4int nloop = 0;                                187   G4int nloop = 0;
191                                                   188 
192   G4double bindingEnergy, ePotEnergy, eKinEner    189   G4double bindingEnergy, ePotEnergy, eKinEnergy;
193   G4double gamEnergy0, gamEnergy1;                190   G4double gamEnergy0, gamEnergy1;
194                                                   191 
195   CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 192   //static const G4double eminus2 =  1.0 - G4Exp(-2.0);
196   G4double rndm[4];                            << 
197                                                   193 
198   do {                                            194   do {
199     ++nloop;                                      195     ++nloop;
200                                                << 196     G4double xprob = totprob*rndmEngineMod->flat();
201     // 4 random numbers to select e-           << 
202     rndmEngineMod->flatArray(4, rndm);         << 
203     G4double xprob = totprob*rndm[0];          << 
204                                                   197 
205     // select shell                               198     // select shell
206     for(i=0; i<nShells; ++i) { if(xprob <= fPr    199     for(i=0; i<nShells; ++i) { if(xprob <= fProbabilities[i]) { break; } }
207                                                   200    
208     bindingEnergy = elm->GetAtomicShell(i);       201     bindingEnergy = elm->GetAtomicShell(i);
209     lv1.set(0.0,0.0,energy,energy);               202     lv1.set(0.0,0.0,energy,energy);
210     /*                                         << 203 
211     G4cout << "nShells= " << nShells << " i= " << 204     //G4cout << "nShells= " << nShells << " i= " << i 
212        << " Egamma= " << energy << " Ebind= "  << 205     //   << " Egamma= " << energy << " Ebind= " << bindingEnergy
213        << G4endl;                              << 206     //   << " Elim= " << limitEnergy 
214     */                                         << 207     //   << G4endl;
                                                   >> 208 
215     // for rest frame of the electron             209     // for rest frame of the electron
216     G4double x = -G4Log(rndm[1]);              << 210     G4double x = -G4Log(rndmEngineMod->flat());
217     eKinEnergy = bindingEnergy*x;                 211     eKinEnergy = bindingEnergy*x;
218     ePotEnergy = bindingEnergy*(1.0 + x);         212     ePotEnergy = bindingEnergy*(1.0 + x);
219                                                   213 
220     // for rest frame of the electron             214     // for rest frame of the electron
221     G4double eTotMomentum = sqrt(eKinEnergy*(e    215     G4double eTotMomentum = sqrt(eKinEnergy*(eKinEnergy + 2*electron_mass_c2));
222     G4double phi = rndm[2]*twopi;              << 216     G4double phi = rndmEngineMod->flat()*twopi;
223     G4double costet = 2*rndm[3] - 1;           << 217     G4double costet = 2*rndmEngineMod->flat() - 1;
224     G4double sintet = sqrt((1 - costet)*(1 + c    218     G4double sintet = sqrt((1 - costet)*(1 + costet));
225     lv2.set(eTotMomentum*sintet*cos(phi),eTotM    219     lv2.set(eTotMomentum*sintet*cos(phi),eTotMomentum*sintet*sin(phi),
226             eTotMomentum*costet,eKinEnergy + e << 220       eTotMomentum*costet,eKinEnergy + electron_mass_c2);
227     bst = lv2.boostVector();                      221     bst = lv2.boostVector();
228     lv1.boost(-bst);                              222     lv1.boost(-bst);
229                                                   223 
230     gamEnergy0 = lv1.e();                         224     gamEnergy0 = lv1.e();
231                                                   225    
232     // In the rest frame of the electron          226     // In the rest frame of the electron
233     // The scattered gamma energy is sampled a    227     // The scattered gamma energy is sampled according to Klein-Nishina formula
234     // The random number techniques of Butcher    228     // The random number techniques of Butcher & Messel are used 
235     // (Nuc Phys 20(1960),15).                    229     // (Nuc Phys 20(1960),15). 
236     G4double E0_m = gamEnergy0/electron_mass_c    230     G4double E0_m = gamEnergy0/electron_mass_c2;
237                                                   231 
238     //G4cout << "Nloop= "<< nloop << " Ecm(keV << 
239     //                                            232     //
240     // sample the energy rate of the scattered    233     // sample the energy rate of the scattered gamma 
241     //                                            234     //
242                                                   235 
243     G4double epsilon, epsilonsq, onecost, sint    236     G4double epsilon, epsilonsq, onecost, sint2, greject ;
244                                                   237 
245     G4double eps0       = 1./(1 + 2*E0_m);        238     G4double eps0       = 1./(1 + 2*E0_m);
246     G4double epsilon0sq = eps0*eps0;              239     G4double epsilon0sq = eps0*eps0;
247     G4double alpha1     = - G4Log(eps0);          240     G4double alpha1     = - G4Log(eps0);
248     G4double alpha2     = alpha1 + 0.5*(1 - ep << 241     G4double alpha2     = 0.5*(1 - epsilon0sq);
249                                                   242 
250     do {                                          243     do {
251       ++nloop;                                    244       ++nloop;
252       // false interaction if too many iterati    245       // false interaction if too many iterations
253       if(nloop > nlooplim) { return; }            246       if(nloop > nlooplim) { return; }
254                                                   247 
255       // 3 random numbers to sample scattering << 248       if ( alpha1/(alpha1+alpha2) > rndmEngineMod->flat() ) {
256       rndmEngineMod->flatArray(3, rndm);       << 249   epsilon   = G4Exp(-alpha1*rndmEngineMod->flat());   // epsilon0**r
257                                                << 250   epsilonsq = epsilon*epsilon; 
258       if ( alpha1 > alpha2*rndm[0] ) {         << 
259         epsilon   = G4Exp(-alpha1*rndm[1]);    << 
260         epsilonsq = epsilon*epsilon;           << 
261                                                   251 
262       } else {                                    252       } else {
263         epsilonsq = epsilon0sq + (1.- epsilon0 << 253   epsilonsq = epsilon0sq + (1.- epsilon0sq)*rndmEngineMod->flat();
264         epsilon   = sqrt(epsilonsq);           << 254   epsilon   = sqrt(epsilonsq);
265       }                                           255       }
266                                                   256 
267       onecost = (1.- epsilon)/(epsilon*E0_m);     257       onecost = (1.- epsilon)/(epsilon*E0_m);
268       sint2   = onecost*(2.-onecost);             258       sint2   = onecost*(2.-onecost);
269       greject = 1. - epsilon*sint2/(1.+ epsilo    259       greject = 1. - epsilon*sint2/(1.+ epsilonsq);
270                                                   260 
271       // Loop checking, 03-Aug-2015, Vladimir  << 261     } while (greject < rndmEngineMod->flat());
272     } while (greject < rndm[2]);               << 
273     gamEnergy1 = epsilon*gamEnergy0;              262     gamEnergy1 = epsilon*gamEnergy0;
274                                                   263  
275     // before scattering total 4-momentum in e    264     // before scattering total 4-momentum in e- system
276     lv2.set(0.0,0.0,0.0,electron_mass_c2);        265     lv2.set(0.0,0.0,0.0,electron_mass_c2);
277     lv2 += lv1;                                   266     lv2 += lv1;
278                                                   267  
279     //                                            268     //
280     // scattered gamma angles. ( Z - axis alon    269     // scattered gamma angles. ( Z - axis along the parent gamma)
281     //                                            270     //
282     if(sint2 < 0.0) { sint2 = 0.0; }              271     if(sint2 < 0.0) { sint2 = 0.0; }
283     costet = 1. - onecost;                        272     costet = 1. - onecost; 
284     sintet = sqrt(sint2);                         273     sintet = sqrt(sint2);
285     phi  = twopi * rndmEngineMod->flat();         274     phi  = twopi * rndmEngineMod->flat();
286                                                   275 
287     // e- recoil                                  276     // e- recoil
288     //                                            277     //
289     // in  rest frame of the electron             278     // in  rest frame of the electron
290     G4ThreeVector gamDir = lv1.vect().unit();     279     G4ThreeVector gamDir = lv1.vect().unit();
291     G4ThreeVector v = G4ThreeVector(sintet*cos    280     G4ThreeVector v = G4ThreeVector(sintet*cos(phi),sintet*sin(phi),costet);
292     v.rotateUz(gamDir);                           281     v.rotateUz(gamDir);
293     lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(),    282     lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(),gamEnergy1*v.z(),gamEnergy1);
294     lv2 -= lv1;                                   283     lv2 -= lv1;
295     //G4cout<<"Egam(keV)= " << lv1.e()/keV     << 284     //G4cout<<"Egam= "<<lv1.e()<<" Ee= "<< lv2.e()-electron_mass_c2 << G4endl;
296     //          <<" Ee(keV)= " << (lv2.e()-ele << 
297     lv2.boost(bst);                               285     lv2.boost(bst);
298     eKinEnergy = lv2.e() - electron_mass_c2 -     286     eKinEnergy = lv2.e() - electron_mass_c2 - ePotEnergy;   
299     //G4cout << "Nloop= " << nloop << " eKinEn << 287     //G4cout << "eKinEnergy= " << eKinEnergy << G4endl;
300                                                   288 
301     // Loop checking, 03-Aug-2015, Vladimir Iv << 
302   } while ( eKinEnergy < 0.0 );                   289   } while ( eKinEnergy < 0.0 );
303                                                   290 
304   //                                              291   //
305   // update G4VParticleChange for the scattere    292   // update G4VParticleChange for the scattered gamma
306   //                                              293   //
307                                                   294    
308   lv1.boost(bst);                                 295   lv1.boost(bst);
309   gamEnergy1 = lv1.e();                           296   gamEnergy1 = lv1.e();
310   if(gamEnergy1 > lowestSecondaryEnergy) {        297   if(gamEnergy1 > lowestSecondaryEnergy) {
311     G4ThreeVector gamDirection1 = lv1.vect().u    298     G4ThreeVector gamDirection1 = lv1.vect().unit();
312     gamDirection1.rotateUz(direction);            299     gamDirection1.rotateUz(direction);
313     fParticleChange->ProposeMomentumDirection(    300     fParticleChange->ProposeMomentumDirection(gamDirection1);
314   } else {                                        301   } else { 
315     fParticleChange->ProposeTrackStatus(fStopA    302     fParticleChange->ProposeTrackStatus(fStopAndKill);
316     gamEnergy1 = 0.0;                             303     gamEnergy1 = 0.0;
317   }                                               304   }
318   fParticleChange->SetProposedKineticEnergy(ga    305   fParticleChange->SetProposedKineticEnergy(gamEnergy1);
319                                                   306 
320   //                                              307   //
321   // kinematic of the scattered electron          308   // kinematic of the scattered electron
322   //                                              309   //
323                                                   310 
324   if(eKinEnergy > lowestSecondaryEnergy) {        311   if(eKinEnergy > lowestSecondaryEnergy) {
325     G4ThreeVector eDirection = lv2.vect().unit    312     G4ThreeVector eDirection = lv2.vect().unit();
326     eDirection.rotateUz(direction);               313     eDirection.rotateUz(direction);
327     auto dp = new G4DynamicParticle(theElectro << 314     G4DynamicParticle* dp = 
                                                   >> 315       new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
328     fvect->push_back(dp);                         316     fvect->push_back(dp);
329   } else { eKinEnergy = 0.0; }                    317   } else { eKinEnergy = 0.0; }
330                                                   318 
331   G4double edep = energy - gamEnergy1 - eKinEn    319   G4double edep = energy - gamEnergy1 - eKinEnergy;
332   G4double esec = 0.0;                            320   G4double esec = 0.0;
333                                                   321   
334   // sample deexcitation                          322   // sample deexcitation
335   //                                              323   //
336   if(nullptr != fAtomDeexcitation) {           << 324   if(fAtomDeexcitation) {
337     G4int index = couple->GetIndex();             325     G4int index = couple->GetIndex();
338     if(fAtomDeexcitation->CheckDeexcitationAct    326     if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
339       G4int Z = elm->GetZasInt();              << 327       G4int Z = G4lrint(elm->GetZ());
340       auto as = (G4AtomicShellEnumerator)(i);  << 328       G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i);
341       const G4AtomicShell* shell = fAtomDeexci    329       const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
342       G4int nbefore = (G4int)fvect->size();    << 330       size_t nbefore = fvect->size();
343       fAtomDeexcitation->GenerateParticles(fve    331       fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
344       G4int nafter = (G4int)fvect->size();     << 332       size_t nafter = fvect->size();
345       //G4cout << "N1= " << nbefore << "  N2=  << 333       if(nafter > nbefore) {
346       for (G4int j=nbefore; j<nafter; ++j) {   << 334   for (size_t j=nbefore; j<nafter; ++j) {
347         G4double e = ((*fvect)[j])->GetKinetic << 335     G4double e = ((*fvect)[j])->GetKineticEnergy();
348         if(esec + e > edep) {                  << 336     if(esec + e > edep) {
349           // correct energy in order to have e << 337       /*   
350           e = edep - esec;                     << 338       G4cout << "### G4KleinNishinaModel Edep(eV)= " << edep/eV 
351           ((*fvect)[j])->SetKineticEnergy(e);  << 339        << " Esec(eV)= " << esec/eV 
352           esec += e;                           << 340        << " E["<< j << "](eV)= " << e/eV
353           /*                                   << 341        << " N= " << nafter
354             G4cout << "### G4KleinNishinaModel << 342        << " Z= " << Z << " shell= " << i 
355                    << " Esec(eV)= " << esec/eV << 343        << "  Ebind(keV)= " << bindingEnergy/keV 
356                    << " E["<< j << "](eV)= " < << 344        << "  Eshell(keV)= " << shell->BindingEnergy()/keV 
357                    << " N= " << nafter         << 345        << G4endl;
358                    << " Z= " << Z << " shell=  << 346       */
359                    << "  Ebind(keV)= " << bind << 347       for (size_t jj=nafter-1; jj>=j; --jj) { 
360                    << "  Eshell(keV)= " << she << 348         delete (*fvect)[jj]; 
361                    << G4endl;                  << 349         fvect->pop_back(); 
362           */                                   << 350       }
363           // delete the rest of secondaries (s << 351       break;        
364           for (G4int jj=nafter-1; jj>j; --jj)  << 352     }
365             delete (*fvect)[jj];               << 353     esec += e;
366             fvect->pop_back();                 << 354   } 
367           }                                    << 
368           break;                               << 
369         }                                      << 
370         esec += e;                             << 
371       }                                           355       }
372       edep -= esec;                               356       edep -= esec;
373     }                                             357     }
374   }                                               358   }
375   if(std::abs(energy - gamEnergy1 - eKinEnergy << 359   if(fabs(energy - gamEnergy1 - eKinEnergy - esec - edep) > eV) {
376     G4cout << "### G4KleinNishinaModel dE(eV)=    360     G4cout << "### G4KleinNishinaModel dE(eV)= " 
377            << (energy - gamEnergy1 - eKinEnerg << 361      << (energy - gamEnergy1 - eKinEnergy - esec - edep)/eV 
378            << " shell= " << i                  << 362      << " shell= " << i 
379            << "  E(keV)= " << energy/keV       << 363      << "  E(keV)= " << energy/keV 
380            << "  Ebind(keV)= " << bindingEnerg << 364      << "  Ebind(keV)= " << bindingEnergy/keV 
381            << "  Eg(keV)= " << gamEnergy1/keV  << 365      << "  Eg(keV)= " << gamEnergy1/keV 
382            << "  Ee(keV)= " << eKinEnergy/keV  << 366      << "  Ee(keV)= " << eKinEnergy/keV 
383            << "  Esec(keV)= " << esec/keV      << 367      << "  Esec(keV)= " << esec/keV 
384            << "  Edep(keV)= " << edep/keV      << 368      << "  Edep(keV)= " << edep/keV 
385            << G4endl;                          << 369      << G4endl;
386   }                                               370   }
387   // energy balance                               371   // energy balance
388   if(edep > 0.0) {                                372   if(edep > 0.0) { 
389     fParticleChange->ProposeLocalEnergyDeposit    373     fParticleChange->ProposeLocalEnergyDeposit(edep);
390   }                                               374   }
391 }                                                 375 }
392                                                   376 
393 //....oooOO0OOooo........oooOO0OOooo........oo    377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
394                                                   378 
395                                                   379