Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4MicroElecElasticModel.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/lowenergy/src/G4MicroElecElasticModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4MicroElecElasticModel.cc (Version 10.7)


  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 //                                                 26 //
 27 // G4MicroElecElasticModel.cc, 2011/08/29 A.Va     27 // G4MicroElecElasticModel.cc, 2011/08/29 A.Valentin, M. Raine
 28 //                                                 28 //
 29 // Based on the following publications             29 // Based on the following publications
 30 //      - Geant4 physics processes for microdo     30 //      - Geant4 physics processes for microdosimetry simulation:
 31 //      very low energy electromagnetic models     31 //      very low energy electromagnetic models for electrons in Si,
 32 //      NIM B, vol. 288, pp. 66 - 73, 2012.        32 //      NIM B, vol. 288, pp. 66 - 73, 2012.
 33 //                                                 33 //
 34 //                                                 34 //
 35 //....oooOO0OOooo........oooOO0OOooo........oo     35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 36                                                    36 
 37                                                    37 
 38 #include "G4MicroElecElasticModel.hh"              38 #include "G4MicroElecElasticModel.hh"
 39 #include "G4PhysicalConstants.hh"                  39 #include "G4PhysicalConstants.hh"
 40 #include "G4SystemOfUnits.hh"                      40 #include "G4SystemOfUnits.hh"
 41 #include "G4Exp.hh"                                41 #include "G4Exp.hh"
 42                                                    42 
 43 //....oooOO0OOooo........oooOO0OOooo........oo     43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 44                                                    44 
 45 using namespace std;                               45 using namespace std;
 46                                                    46 
 47 //....oooOO0OOooo........oooOO0OOooo........oo     47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 48                                                    48 
 49 G4MicroElecElasticModel::G4MicroElecElasticMod     49 G4MicroElecElasticModel::G4MicroElecElasticModel(const G4ParticleDefinition*,
 50              const G4String& nam)              <<  50                                              const G4String& nam)
 51  :G4VEmModel(nam),isInitialised(false)         <<  51 :G4VEmModel(nam),isInitialised(false)
 52 {                                                  52 {
 53   nistSi = G4NistManager::Instance()->FindOrBu     53   nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
 54                                                    54 
 55   killBelowEnergy = 16.7 * eV; // Minimum e- e     55   killBelowEnergy = 16.7 * eV; // Minimum e- energy for energy loss by excitation
 56   lowEnergyLimit = 0 * eV;                         56   lowEnergyLimit = 0 * eV;
 57   lowEnergyLimitOfModel = 5 * eV; // The model     57   lowEnergyLimitOfModel = 5 * eV; // The model lower energy is 5 eV
 58   highEnergyLimit = 100. * MeV;                    58   highEnergyLimit = 100. * MeV;
 59   SetLowEnergyLimit(lowEnergyLimit);               59   SetLowEnergyLimit(lowEnergyLimit);
 60   SetHighEnergyLimit(highEnergyLimit);             60   SetHighEnergyLimit(highEnergyLimit);
 61                                                    61 
 62   verboseLevel= 0;                                 62   verboseLevel= 0;
 63   // Verbosity scale:                              63   // Verbosity scale:
 64   // 0 = nothing                                   64   // 0 = nothing
 65   // 1 = warning for energy non-conservation       65   // 1 = warning for energy non-conservation
 66   // 2 = details of energy budget                  66   // 2 = details of energy budget
 67   // 3 = calculation of cross sections, file o     67   // 3 = calculation of cross sections, file openings, sampling of atoms
 68   // 4 = entering in methods                       68   // 4 = entering in methods
 69                                                    69 
 70   if( verboseLevel>0 )                             70   if( verboseLevel>0 )
 71     {                                          <<  71   {
 72       G4cout << "MicroElec Elastic model is co <<  72     G4cout << "MicroElec Elastic model is constructed " << G4endl
 73        << "Energy range: "                     <<  73            << "Energy range: "
 74        << lowEnergyLimit / eV << " eV - "      <<  74            << lowEnergyLimit / eV << " eV - "
 75        << highEnergyLimit / MeV << " MeV"      <<  75            << highEnergyLimit / MeV << " MeV"
 76        << G4endl;                              <<  76            << G4endl;
 77     }                                          <<  77   }
 78   fParticleChangeForGamma = 0;                     78   fParticleChangeForGamma = 0;
 79 }                                                  79 }
 80                                                    80 
 81 //....oooOO0OOooo........oooOO0OOooo........oo     81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 82                                                    82 
 83 G4MicroElecElasticModel::~G4MicroElecElasticMo     83 G4MicroElecElasticModel::~G4MicroElecElasticModel()
 84 {                                                  84 {
 85   // For total cross section                       85   // For total cross section
 86   for (auto & pos : tableData)                 << 
 87     {                                          << 
 88       G4MicroElecCrossSectionDataSet* table =  << 
 89       delete table;                            << 
 90     }                                          << 
 91                                                    86 
 92   // For final state                           <<  87   std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
 93   eVecm.clear();                               <<  88   for (pos = tableData.begin(); pos != tableData.end(); ++pos)
                                                   >>  89   {
                                                   >>  90     G4MicroElecCrossSectionDataSet* table = pos->second;
                                                   >>  91     delete table;
                                                   >>  92   }
                                                   >>  93 
                                                   >>  94    // For final state
                                                   >>  95 
                                                   >>  96    eVecm.clear();
                                                   >>  97 
 94 }                                                  98 }
 95                                                    99 
 96 //....oooOO0OOooo........oooOO0OOooo........oo    100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 97                                                   101 
 98 void G4MicroElecElasticModel::Initialise(const    102 void G4MicroElecElasticModel::Initialise(const G4ParticleDefinition* /*particle*/,
 99            const G4DataVector& /*cuts*/)       << 103                                        const G4DataVector& /*cuts*/)
100 {                                                 104 {
                                                   >> 105 
101   if (verboseLevel > 3)                           106   if (verboseLevel > 3)
102     G4cout << "Calling G4MicroElecElasticModel    107     G4cout << "Calling G4MicroElecElasticModel::Initialise()" << G4endl;
103                                                   108 
104   // Energy limits                                109   // Energy limits
                                                   >> 110 
105   if (LowEnergyLimit() < lowEnergyLimit)          111   if (LowEnergyLimit() < lowEnergyLimit)
106     {                                          << 112   {
107       G4cout << "G4MicroElecElasticModel: low  << 113     G4cout << "G4MicroElecElasticModel: low energy limit increased from " <<
108   LowEnergyLimit()/eV << " eV to " << lowEnerg    114   LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
109       SetLowEnergyLimit(lowEnergyLimit);       << 115     SetLowEnergyLimit(lowEnergyLimit);
110     }                                             116     }
111                                                   117 
112   if (HighEnergyLimit() > highEnergyLimit)        118   if (HighEnergyLimit() > highEnergyLimit)
113     {                                          << 119   {
114       G4cout << "G4MicroElecElasticModel: high << 120     G4cout << "G4MicroElecElasticModel: high energy limit decreased from " <<
115         HighEnergyLimit()/MeV << " MeV to " <<    121         HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
116       SetHighEnergyLimit(highEnergyLimit);     << 122     SetHighEnergyLimit(highEnergyLimit);
117     }                                          << 123   }
118                                                   124 
119   // Reading of data files                        125   // Reading of data files
120                                                   126 
121   G4double scaleFactor = 1e-18 * cm * cm;         127   G4double scaleFactor = 1e-18 * cm * cm;
                                                   >> 128 
122   G4String fileElectron("microelec/sigma_elast    129   G4String fileElectron("microelec/sigma_elastic_e_Si");
123                                                   130 
124   G4ParticleDefinition* electronDef = G4Electr    131   G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
125   G4String electron;                              132   G4String electron;
126                                                   133 
127   // For total cross section                   << 134     // For total cross section
128   electron = electronDef->GetParticleName();   << 
129   tableFile[electron] = fileElectron;          << 
130                                                   135 
131   G4MicroElecCrossSectionDataSet* tableE = new << 136     electron = electronDef->GetParticleName();
132   tableE->LoadData(fileElectron);              << 
133   tableData[electron] = tableE;                << 
134                                                   137 
135   // For final state                           << 138     tableFile[electron] = fileElectron;
136   const char* path = G4FindDataDir("G4LEDATA") << 
137                                                   139 
138   if (!path)                                   << 140     G4MicroElecCrossSectionDataSet* tableE = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
                                                   >> 141     tableE->LoadData(fileElectron);
                                                   >> 142     tableData[electron] = tableE;
                                                   >> 143 
                                                   >> 144     // For final state
                                                   >> 145 
                                                   >> 146     char *path = std::getenv("G4LEDATA");
                                                   >> 147 
                                                   >> 148     if (!path)
139     {                                             149     {
140       G4Exception("G4MicroElecElasticModel::In    150       G4Exception("G4MicroElecElasticModel::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
141       return;                                     151       return;
142     }                                             152     }
143                                                   153 
144   std::ostringstream eFullFileName;            << 154     std::ostringstream eFullFileName;
145   eFullFileName << path << "/microelec/sigmadi << 155     eFullFileName << path << "/microelec/sigmadiff_cumulated_elastic_e_Si.dat";
146   std::ifstream eDiffCrossSection(eFullFileNam << 156     std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
147                                                << 157 
148   if (!eDiffCrossSection)                      << 158     if (!eDiffCrossSection)
149     G4Exception("G4MicroElecElasticModel::Init << 159   G4Exception("G4MicroElecElasticModel::Initialise","em0003",FatalException,"Missing data file: /microelec/sigmadiff_cumulated_elastic_e_Si.dat");
150     "Missing data file: /microelec/sigmadiff_c << 160 
151                                                << 161 
152   // Added clear for MT                        << 162     // October 21th, 2014 - Melanie Raine
153   eTdummyVec.clear();                          << 163     // Added clear for MT
154   eVecm.clear();                               << 164 
155   eDiffCrossSectionData.clear();               << 165     eTdummyVec.clear();
                                                   >> 166     eVecm.clear();
                                                   >> 167     eDiffCrossSectionData.clear();
156                                                   168 
157   //                                           << 169     //
158   eTdummyVec.push_back(0.);                    << 
159                                                   170 
160   while(!eDiffCrossSection.eof())              << 171 
                                                   >> 172     eTdummyVec.push_back(0.);
                                                   >> 173 
                                                   >> 174     while(!eDiffCrossSection.eof())
161     {                                             175     {
162       double tDummy;                           << 176   double tDummy;
163       double eDummy;                           << 177   double eDummy;
164       eDiffCrossSection>>tDummy>>eDummy;       << 178   eDiffCrossSection>>tDummy>>eDummy;
165                                                   179 
166       if (tDummy != eTdummyVec.back())         << 180   // SI : mandatory eVecm initialization
                                                   >> 181 
                                                   >> 182         if (tDummy != eTdummyVec.back())
167         {                                         183         {
168           eTdummyVec.push_back(tDummy);           184           eTdummyVec.push_back(tDummy);
169           eVecm[tDummy].push_back(0.);            185           eVecm[tDummy].push_back(0.);
170         }                                         186         }
171                                                   187 
172       eDiffCrossSection>>eDiffCrossSectionData << 188         eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
                                                   >> 189 
                                                   >> 190         if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
173                                                   191 
174       if (eDummy != eVecm[tDummy].back()) eVec << 
175     }                                             192     }
176   // End final state                           << 193 
                                                   >> 194     // End final state
177                                                   195 
178   if (verboseLevel > 2)                           196   if (verboseLevel > 2)
179     G4cout << "Loaded cross section files for     197     G4cout << "Loaded cross section files for MicroElec Elastic model" << G4endl;
180                                                   198 
181   if( verboseLevel>0 )                            199   if( verboseLevel>0 )
182     {                                          << 200   {
183       G4cout << "MicroElec Elastic model is in << 201     G4cout << "MicroElec Elastic model is initialized " << G4endl
184        << "Energy range: "                     << 202            << "Energy range: "
185        << LowEnergyLimit() / eV << " eV - "    << 203            << LowEnergyLimit() / eV << " eV - "
186        << HighEnergyLimit() / MeV << " MeV"    << 204            << HighEnergyLimit() / MeV << " MeV"
187        << G4endl;                              << 205            << G4endl;
188     }                                          << 206   }
189                                                   207 
190   if (isInitialised) { return; }                  208   if (isInitialised) { return; }
191   fParticleChangeForGamma = GetParticleChangeF    209   fParticleChangeForGamma = GetParticleChangeForGamma();
192   isInitialised = true;                           210   isInitialised = true;
                                                   >> 211 
193 }                                                 212 }
194                                                   213 
195 //....oooOO0OOooo........oooOO0OOooo........oo    214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196                                                   215 
197 G4double G4MicroElecElasticModel::CrossSection    216 G4double G4MicroElecElasticModel::CrossSectionPerVolume(const G4Material* material,
198               const G4ParticleDefinition* p,   << 217              const G4ParticleDefinition* p,
199               G4double ekin,                   << 218              G4double ekin,
200               G4double,                        << 219              G4double,
201               G4double)                        << 220              G4double)
202 {                                                 221 {
203   if (verboseLevel > 3)                           222   if (verboseLevel > 3)
204     G4cout << "Calling CrossSectionPerVolume()    223     G4cout << "Calling CrossSectionPerVolume() of G4MicroElecElasticModel" << G4endl;
205                                                   224 
206   // Calculate total cross section for model   << 225  // Calculate total cross section for model
207   G4double sigma=0;                            << 
208   G4double density = material->GetTotNbOfAtoms << 
209                                                   226 
210   if (material == nistSi || material->GetBaseM << 227  G4double sigma=0;
211     {                                          << 228 
212       const G4String& particleName = p->GetPar << 229  G4double density = material->GetTotNbOfAtomsPerVolume();
                                                   >> 230 
                                                   >> 231  if (material == nistSi || material->GetBaseMaterial() == nistSi)
                                                   >> 232  {
                                                   >> 233   const G4String& particleName = p->GetParticleName();
213                                                   234 
214       if (ekin < highEnergyLimit)              << 235   if (ekin < highEnergyLimit)
                                                   >> 236   {
                                                   >> 237       //SI : XS must not be zero otherwise sampling of secondaries method ignored
                                                   >> 238       if (ekin < killBelowEnergy) return DBL_MAX;
                                                   >> 239       //
                                                   >> 240 
                                                   >> 241   std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
                                                   >> 242   pos = tableData.find(particleName);
                                                   >> 243 
                                                   >> 244   if (pos != tableData.end())
215   {                                               245   {
216     //SI : XS must not be zero otherwise sampl << 246     G4MicroElecCrossSectionDataSet* table = pos->second;
217     if (ekin < killBelowEnergy) return DBL_MAX << 247     if (table != 0)
218     //                                         << 248     {
219                                                << 249       sigma = table->FindValue(ekin);
220     auto pos = tableData.find(particleName);   << 250     }
221     if (pos != tableData.end())                << 
222       {                                        << 
223         G4MicroElecCrossSectionDataSet* table  << 
224         if (table != nullptr)                  << 
225     {                                          << 
226       sigma = table->FindValue(ekin);          << 
227     }                                          << 
228       }                                        << 
229     else                                       << 
230       {                                        << 
231         G4Exception("G4MicroElecElasticModel:: << 
232         FatalException,"Model not applicable t << 
233       }                                        << 
234   }                                               251   }
235                                                << 252   else
236       if (verboseLevel > 3)                    << 
237   {                                               253   {
238     G4cout << "---> Kinetic energy(eV)=" << ek << 254       G4Exception("G4MicroElecElasticModel::ComputeCrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
239     G4cout << " - Cross section per Si atom (c << 
240     G4cout << " - Cross section per Si atom (c << 
241   }                                               255   }
242     }                                          << 256   }
243   return sigma*density;                        << 257 
                                                   >> 258   if (verboseLevel > 3)
                                                   >> 259   {
                                                   >> 260     G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
                                                   >> 261     G4cout << " - Cross section per Si atom (cm^2)=" << sigma/cm/cm << G4endl;
                                                   >> 262     G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density/(1./cm) << G4endl;
                                                   >> 263   }
                                                   >> 264 
                                                   >> 265  }
                                                   >> 266 
                                                   >> 267  return sigma*density;
244 }                                                 268 }
245                                                   269 
246 //....oooOO0OOooo........oooOO0OOooo........oo    270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
247                                                   271 
248 void G4MicroElecElasticModel::SampleSecondarie    272 void G4MicroElecElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
249             const G4MaterialCutsCouple* /*coup << 273                 const G4MaterialCutsCouple* /*couple*/,
250             const G4DynamicParticle* aDynamicE << 274                 const G4DynamicParticle* aDynamicElectron,
251             G4double,                          << 275                 G4double,
252             G4double)                          << 276                 G4double)
253 {                                                 277 {
254                                                   278 
255   if (verboseLevel > 3)                           279   if (verboseLevel > 3)
256     G4cout << "Calling SampleSecondaries() of     280     G4cout << "Calling SampleSecondaries() of G4MicroElecElasticModel" << G4endl;
257                                                   281 
258   G4double electronEnergy0 = aDynamicElectron-    282   G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
259                                                   283 
260   if (electronEnergy0 < killBelowEnergy)          284   if (electronEnergy0 < killBelowEnergy)
261     {                                          << 285   {
262       fParticleChangeForGamma->SetProposedKine << 286     fParticleChangeForGamma->SetProposedKineticEnergy(0.);
263       fParticleChangeForGamma->ProposeTrackSta << 287     fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
264       fParticleChangeForGamma->ProposeLocalEne << 288     fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
265       return ;                                 << 289     return ;
266     }                                          << 290   }
267                                                   291 
268   if (electronEnergy0>= killBelowEnergy && ele    292   if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
269     {                                          << 293   {
270       G4double cosTheta = RandomizeCosTheta(el << 294     G4double cosTheta = RandomizeCosTheta(electronEnergy0);
271       G4double phi = 2. * pi * G4UniformRand() << 
272       G4ThreeVector zVers = aDynamicElectron-> << 
273       G4ThreeVector xVers = zVers.orthogonal() << 
274       G4ThreeVector yVers = zVers.cross(xVers) << 
275                                                << 
276       G4double xDir = std::sqrt(1. - cosTheta* << 
277       G4double yDir = xDir;                    << 
278       xDir *= std::cos(phi);                   << 
279       yDir *= std::sin(phi);                   << 
280                                                   295 
281       G4ThreeVector zPrimeVers((xDir*xVers + y << 296     G4double phi = 2. * pi * G4UniformRand();
                                                   >> 297 
                                                   >> 298     G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
                                                   >> 299     G4ThreeVector xVers = zVers.orthogonal();
                                                   >> 300     G4ThreeVector yVers = zVers.cross(xVers);
                                                   >> 301 
                                                   >> 302     G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
                                                   >> 303     G4double yDir = xDir;
                                                   >> 304     xDir *= std::cos(phi);
                                                   >> 305     yDir *= std::sin(phi);
                                                   >> 306 
                                                   >> 307     G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
                                                   >> 308 
                                                   >> 309     fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
                                                   >> 310 
                                                   >> 311     fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
                                                   >> 312   }
282                                                   313 
283       fParticleChangeForGamma->ProposeMomentum << 
284       fParticleChangeForGamma->SetProposedKine << 
285     }                                          << 
286 }                                                 314 }
287                                                   315 
288 //....oooOO0OOooo........oooOO0OOooo........oo    316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289                                                   317 
290 G4double G4MicroElecElasticModel::Theta           318 G4double G4MicroElecElasticModel::Theta
291 (G4ParticleDefinition * particleDefinition, G4 << 319   (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
292 {                                                 320 {
                                                   >> 321 
293   G4double theta = 0.;                            322   G4double theta = 0.;
294   G4double valueT1 = 0;                           323   G4double valueT1 = 0;
295   G4double valueT2 = 0;                           324   G4double valueT2 = 0;
296   G4double valueE21 = 0;                          325   G4double valueE21 = 0;
297   G4double valueE22 = 0;                          326   G4double valueE22 = 0;
298   G4double valueE12 = 0;                          327   G4double valueE12 = 0;
299   G4double valueE11 = 0;                          328   G4double valueE11 = 0;
300   G4double xs11 = 0;                              329   G4double xs11 = 0;
301   G4double xs12 = 0;                              330   G4double xs12 = 0;
302   G4double xs21 = 0;                              331   G4double xs21 = 0;
303   G4double xs22 = 0;                              332   G4double xs22 = 0;
304                                                   333 
                                                   >> 334 
305   if (particleDefinition == G4Electron::Electr    335   if (particleDefinition == G4Electron::ElectronDefinition())
306     {                                          << 336   {
307       auto t2 = std::upper_bound(eTdummyVec.be << 337     std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
308       auto t1 = t2-1;                          << 338     std::vector<double>::iterator t1 = t2-1;
309       auto e12 = std::upper_bound(eVecm[(*t1)] << 339 
310       auto e11 = e12-1;                        << 340     std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff);
311       auto e22 = std::upper_bound(eVecm[(*t2)] << 341     std::vector<double>::iterator e11 = e12-1;
312       auto e21 = e22-1;                        << 342 
313                                                << 343     std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff);
314       valueT1  =*t1;                           << 344     std::vector<double>::iterator e21 = e22-1;
315       valueT2  =*t2;                           << 345 
316       valueE21 =*e21;                          << 346     valueT1  =*t1;
317       valueE22 =*e22;                          << 347     valueT2  =*t2;
318       valueE12 =*e12;                          << 348     valueE21 =*e21;
319       valueE11 =*e11;                          << 349     valueE22 =*e22;
320                                                << 350     valueE12 =*e12;
321       xs11 = eDiffCrossSectionData[valueT1][va << 351     valueE11 =*e11;
322       xs12 = eDiffCrossSectionData[valueT1][va << 352 
323       xs21 = eDiffCrossSectionData[valueT2][va << 353     xs11 = eDiffCrossSectionData[valueT1][valueE11];
324       xs22 = eDiffCrossSectionData[valueT2][va << 354     xs12 = eDiffCrossSectionData[valueT1][valueE12];
325     }                                          << 355     xs21 = eDiffCrossSectionData[valueT2][valueE21];
                                                   >> 356     xs22 = eDiffCrossSectionData[valueT2][valueE22];
                                                   >> 357 
                                                   >> 358 }
                                                   >> 359 
326   if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)     360   if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
327                                                   361 
328   theta = QuadInterpolator(  valueE11, valueE1    362   theta = QuadInterpolator(  valueE11, valueE12,
329                valueE21, valueE22,                363                valueE21, valueE22,
330            xs11, xs12,                            364            xs11, xs12,
331            xs21, xs22,                            365            xs21, xs22,
332            valueT1, valueT2,                      366            valueT1, valueT2,
333            k, integrDiff );                       367            k, integrDiff );
334                                                   368 
335   return theta;                                   369   return theta;
336 }                                                 370 }
337                                                   371 
338 //....oooOO0OOooo........oooOO0OOooo........oo    372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339                                                   373 
340 G4double G4MicroElecElasticModel::LinLogInterp    374 G4double G4MicroElecElasticModel::LinLogInterpolate(G4double e1,
341                 G4double e2,                   << 375                     G4double e2,
342                 G4double e,                    << 376                     G4double e,
343                 G4double xs1,                  << 377                     G4double xs1,
344                 G4double xs2)                  << 378                     G4double xs2)
345 {                                                 379 {
346   G4double d1 = std::log(xs1);                    380   G4double d1 = std::log(xs1);
347   G4double d2 = std::log(xs2);                    381   G4double d2 = std::log(xs2);
348   G4double value = G4Exp(d1 + (d2 - d1)*(e - e    382   G4double value = G4Exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
349   return value;                                   383   return value;
350 }                                                 384 }
351                                                   385 
352 //....oooOO0OOooo........oooOO0OOooo........oo    386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
353                                                   387 
354 G4double G4MicroElecElasticModel::LinLinInterp    388 G4double G4MicroElecElasticModel::LinLinInterpolate(G4double e1,
355                 G4double e2,                   << 389                     G4double e2,
356                 G4double e,                    << 390                     G4double e,
357                 G4double xs1,                  << 391                     G4double xs1,
358                 G4double xs2)                  << 392                     G4double xs2)
359 {                                                 393 {
360   G4double d1 = xs1;                              394   G4double d1 = xs1;
361   G4double d2 = xs2;                              395   G4double d2 = xs2;
362   G4double value = (d1 + (d2 - d1)*(e - e1)/ (    396   G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
363   return value;                                   397   return value;
364 }                                                 398 }
365                                                   399 
366 //....oooOO0OOooo........oooOO0OOooo........oo    400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
367                                                   401 
368 G4double G4MicroElecElasticModel::LogLogInterp    402 G4double G4MicroElecElasticModel::LogLogInterpolate(G4double e1,
369                 G4double e2,                   << 403                     G4double e2,
370                 G4double e,                    << 404                     G4double e,
371                 G4double xs1,                  << 405                     G4double xs1,
372                 G4double xs2)                  << 406                     G4double xs2)
373 {                                                 407 {
374   G4double a = (std::log10(xs2)-std::log10(xs1    408   G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
375   G4double b = std::log10(xs2) - a*std::log10(    409   G4double b = std::log10(xs2) - a*std::log10(e2);
376   G4double sigma = a*std::log10(e) + b;           410   G4double sigma = a*std::log10(e) + b;
377   G4double value = (std::pow(10.,sigma));         411   G4double value = (std::pow(10.,sigma));
378   return value;                                   412   return value;
379 }                                                 413 }
380                                                   414 
381 //....oooOO0OOooo........oooOO0OOooo........oo    415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
382                                                   416 
383 G4double G4MicroElecElasticModel::QuadInterpol    417 G4double G4MicroElecElasticModel::QuadInterpolator(G4double e11, G4double e12,
384                G4double e21, G4double e22,     << 418                    G4double e21, G4double e22,
385                G4double xs11, G4double xs12,   << 419                    G4double xs11, G4double xs12,
386                G4double xs21, G4double xs22,   << 420                    G4double xs21, G4double xs22,
387                G4double t1, G4double t2,       << 421                    G4double t1, G4double t2,
388                G4double t, G4double e)         << 422                    G4double t, G4double e)
389 {                                              << 423 {
390   // Log-Log                                   << 424  // Log-Log
391   /*                                           << 425 /*
392     G4double interpolatedvalue1 = LogLogInterp << 426   G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
393     G4double interpolatedvalue2 = LogLogInterp << 427   G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
394     G4double value = LogLogInterpolate(t1, t2, << 428   G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
395                                                << 429 
396                                                << 430 
397     // Lin-Log                                 << 431   // Lin-Log
398     G4double interpolatedvalue1 = LinLogInterp << 432   G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
399     G4double interpolatedvalue2 = LinLogInterp << 433   G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
400     G4double value = LinLogInterpolate(t1, t2, << 434   G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
401   */                                           << 435 */
402                                                   436 
403   // Lin-Lin                                      437   // Lin-Lin
404   G4double interpolatedvalue1 = LinLinInterpol    438   G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
405   G4double interpolatedvalue2 = LinLinInterpol    439   G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
406   G4double value = LinLinInterpolate(t1, t2, t    440   G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
407                                                   441 
408   return value;                                   442   return value;
409 }                                                 443 }
410                                                   444 
411 //....oooOO0OOooo........oooOO0OOooo........oo    445 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
412                                                   446 
413 G4double G4MicroElecElasticModel::RandomizeCos    447 G4double G4MicroElecElasticModel::RandomizeCosTheta(G4double k)
414 {                                                 448 {
415   G4double integrdiff=0;                          449   G4double integrdiff=0;
416   G4double uniformRand=G4UniformRand();        << 450  G4double uniformRand=G4UniformRand();
417   integrdiff = uniformRand;                    << 451  integrdiff = uniformRand;
418                                                   452 
419   G4double theta=0.;                           << 453  G4double theta=0.;
420   G4double cosTheta=0.;                        << 454  G4double cosTheta=0.;
421   theta = Theta(G4Electron::ElectronDefinition << 455  theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
422                                                   456 
423   cosTheta= std::cos(theta*pi/180);            << 457  cosTheta= std::cos(theta*pi/180);
424                                                   458 
425   return cosTheta;                             << 459  return cosTheta;
426 }                                                 460 }
427                                                   461