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 9.6.p2)


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