Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4MicroElecSurface.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/G4MicroElecSurface.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4MicroElecSurface.cc (Version 11.1.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 //                                                 26 //
 27 // G4MicroElecSurface.cc,                          27 // G4MicroElecSurface.cc, 
 28 //               2020/05/20 P. Caron, C. Ingui <<  28 //                          2020/05/20 P. Caron, C. Inguimbert are with ONERA [b] 
 29 //                          Q. Gibaru is with  <<  29 //                   Q. Gibaru is with CEA [a], ONERA [b] and CNES [c]
 30 //                          M. Raine and D. La <<  30 //                   M. Raine and D. Lambert are with CEA [a]
 31 //                                                 31 //
 32 // A part of this work has been funded by the      32 // A part of this work has been funded by the French space agency(CNES[c])
 33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France       33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France
 34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 T     34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 Toulouse, France
 35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CED     35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CEDEX, France
 36 //                                                 36 //
 37 // Based on the following publications             37 // Based on the following publications
 38 //                                                 38 //
 39 // - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine, <<  39 //  - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine, D.Lambert, J.Puech, 
 40 //   Geant4 physics processes for microdosimet <<  40 //        Geant4 physics processes for microdosimetry and secondary electron emission simulation : 
 41 //   Extension of MicroElec to very low energi <<  41 //        Extension of MicroElec to very low energies and new materials
 42 //   NIM B, 2020, in review.                   <<  42 //        NIM B, 2020, in review.
 43 //                                                 43 //
 44 // - Modele de transport d'electrons a basse e <<  44 //
 45 //   applications spatiales (OSMOSEE, GEANT4), <<  45 //      - Modèle de transport d'électrons à basse énergie (10 eV- 2 keV) pour
                                                   >>  46 //        applications spatiales (OSMOSEE, GEANT4), PhD dissertation, 2017. 
                                                   >>  47 //      
 46 //                                                 48 //
 47 //////////////////////////////////////////////     49 ////////////////////////////////////////////////////////////////////////
 48                                                    50 
 49 #include "G4MicroElecSurface.hh"                   51 #include "G4MicroElecSurface.hh"
                                                   >>  52 
 50 #include "G4ios.hh"                                53 #include "G4ios.hh"
 51 #include "G4PhysicalConstants.hh"                  54 #include "G4PhysicalConstants.hh"
 52 #include "G4EmProcessSubType.hh"                   55 #include "G4EmProcessSubType.hh"
 53 #include "G4GeometryTolerance.hh"                  56 #include "G4GeometryTolerance.hh"
 54 #include "G4SystemOfUnits.hh"                  <<  57 #include "G4EmProcessSubType.hh"
 55                                                    58 
                                                   >>  59 #include "G4SystemOfUnits.hh"
 56                                                    60 
 57 //....oooOO0OOooo........oooOO0OOooo........oo     61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 58                                                    62 
 59 G4MicroElecSurface::G4MicroElecSurface(const G     63 G4MicroElecSurface::G4MicroElecSurface(const G4String& processName,G4ProcessType type)
 60   : G4VDiscreteProcess(processName, type),         64   : G4VDiscreteProcess(processName, type),
 61     oldMomentum(0.,0.,0.), previousMomentum(0.     65     oldMomentum(0.,0.,0.), previousMomentum(0.,0.,0.),
 62     theGlobalNormal(0.,0.,0.), theFacetNormal(     66     theGlobalNormal(0.,0.,0.), theFacetNormal(0.,0.,0.)
 63 {                                                  67 {
 64   if ( verboseLevel > 0)                           68   if ( verboseLevel > 0) 
 65   {                                            <<  69     {
 66     G4cout << GetProcessName() << " is created <<  70       G4cout << GetProcessName() << " is created " << G4endl;
 67   }                                            <<  71     }
 68                                                    72   
 69   isInitialised=false;                             73   isInitialised=false;
 70   SetProcessSubType(fSurfaceReflection);           74   SetProcessSubType(fSurfaceReflection);
 71                                                    75   
 72   theStatus = UndefinedSurf;                       76   theStatus = UndefinedSurf;
 73   material1 = nullptr;                             77   material1 = nullptr;
 74   material2 = nullptr;                             78   material2 = nullptr;
 75                                                    79   
 76   kCarTolerance = G4GeometryTolerance::GetInst     80   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 
 77   theParticleMomentum = 0.;                        81   theParticleMomentum = 0.;
 78                                                    82   
 79   flag_franchissement_surface = false;             83   flag_franchissement_surface = false;
 80   flag_normal = false;                             84   flag_normal = false;
 81   flag_reflexion = false;                          85   flag_reflexion = false;
 82   teleportToDo = teleportDone = false;             86   teleportToDo = teleportDone = false;
 83                                                    87 
 84   ekint = thetat = thetaft = energyThreshold =     88   ekint = thetat = thetaft = energyThreshold = crossingProbability = 0.0;
 85 }                                                  89 }
 86                                                    90 
 87 //....oooOO0OOooo........oooOO0OOooo........oo     91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 88                                                    92 
 89 G4MicroElecSurface::~G4MicroElecSurface()          93 G4MicroElecSurface::~G4MicroElecSurface()
 90 {}                                                 94 {}
 91                                                    95 
 92 //....oooOO0OOooo........oooOO0OOooo........oo     96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 93                                                    97 
 94 G4bool                                             98 G4bool 
 95 G4MicroElecSurface::IsApplicable(const G4Parti     99 G4MicroElecSurface::IsApplicable(const G4ParticleDefinition& aParticleType) 
 96 {                                                 100 { 
 97   return ( aParticleType.GetPDGEncoding() == 1    101   return ( aParticleType.GetPDGEncoding() == 11 ); 
 98 }                                                 102 } 
 99                                                   103 
100 //....oooOO0OOooo........oooOO0OOooo........oo    104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101                                                   105 
102 void G4MicroElecSurface::Initialise()          << 
103 {                                              << 
104   if (isInitialised) { return; }               << 
105                                                << 
106   G4ProductionCutsTable* theCoupleTable =      << 
107     G4ProductionCutsTable::GetProductionCutsTa << 
108   G4int numOfCouples = (G4int)theCoupleTable-> << 
109   G4cout << numOfCouples << G4endl;            << 
110                                                << 
111   for (G4int i = 0; i < numOfCouples; ++i)     << 
112   {                                            << 
113     const G4Material* material =               << 
114       theCoupleTable->GetMaterialCutsCouple(i) << 
115                                                << 
116     G4cout << this->GetProcessName() << ", Mat << 
117            << " / " << numOfCouples << " : " < << 
118     if (material->GetName() == "Vacuum")       << 
119     {                                          << 
120       tableWF[material->GetName()] = 0; contin << 
121     }                                          << 
122     G4String mat = material->GetName();        << 
123     G4MicroElecMaterialStructure str = G4Micro << 
124     tableWF[mat] = str.GetWorkFunction();      << 
125   }                                            << 
126                                                << 
127   isInitialised = true;                        << 
128 }                                              << 
129                                                << 
130 //....oooOO0OOooo........oooOO0OOooo........oo << 
131                                                << 
132 void G4MicroElecSurface::BuildPhysicsTable(con    106 void G4MicroElecSurface::BuildPhysicsTable(const G4ParticleDefinition&)
133 {                                                 107 {
134   if (isInitialised) { return; }                  108   if (isInitialised) { return; }
135                                                   109   
136   G4ProductionCutsTable* theCoupleTable =         110   G4ProductionCutsTable* theCoupleTable =
137     G4ProductionCutsTable::GetProductionCutsTa    111     G4ProductionCutsTable::GetProductionCutsTable();
138   G4int numOfCouples = (G4int)theCoupleTable->    112   G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
139   G4cout << "G4MicroElecSurface::Initialise: N    113   G4cout << "G4MicroElecSurface::Initialise: Ncouples= " 
140          << numOfCouples << G4endl;               114          << numOfCouples << G4endl;
141                                                   115   
142   for (G4int i = 0; i < numOfCouples; ++i)     << 116   for (G4int i = 0; i < numOfCouples; ++i) {
143   {                                            << 
144     const G4Material* material =                  117     const G4Material* material =
145       theCoupleTable->GetMaterialCutsCouple(i)    118       theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
146                                                   119     
147     G4cout << "G4Surface, Material " << i + 1  << 120     G4cout << "G4Surface, Material " << i + 1 << " / " << numOfCouples << " : " << material->GetName() << G4endl;
148            << numOfCouples << " : " << materia << 121     if (material->GetName() == "Vacuum") { tableWF[material->GetName()] = 0; continue; }
149     if (material->GetName() == "Vacuum")       << 
150     {                                          << 
151       tableWF[material->GetName()] = 0;        << 
152       continue;                                << 
153     }                                          << 
154     G4String mat = material->GetName();           122     G4String mat = material->GetName();
155     G4MicroElecMaterialStructure str = G4Micro    123     G4MicroElecMaterialStructure str = G4MicroElecMaterialStructure(mat);
156     tableWF[mat] = str.GetWorkFunction();         124     tableWF[mat] = str.GetWorkFunction();
157   }                                               125   }
158   isInitialised = true;                           126   isInitialised = true;
159 }                                                 127 }
160                                                   128 
161 //....oooOO0OOooo........oooOO0OOooo........oo    129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
162                                                   130 
163 G4VParticleChange* G4MicroElecSurface::PostSte    131 G4VParticleChange* G4MicroElecSurface::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
164 {                                                 132 {
165                                                << 
166   if (!isInitialised) { Initialise(); }        << 
167   theStatus = UndefinedSurf;                      133   theStatus = UndefinedSurf;
168                                                << 134   
169   // Definition of the parameters for the part << 135   //Definition of the parameters for the particle
170   aParticleChange.Initialize(aTrack);             136   aParticleChange.Initialize(aTrack);
171   aParticleChange.ProposeVelocity(aTrack.GetVe    137   aParticleChange.ProposeVelocity(aTrack.GetVelocity());
172                                                << 138   
173   G4StepPoint* pPreStepPoint  = aStep.GetPreSt    139   G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
174   G4StepPoint* pPostStepPoint = aStep.GetPostS    140   G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
175                                                << 141   
176   material1 = pPreStepPoint  -> GetMaterial();    142   material1 = pPreStepPoint  -> GetMaterial();
177   material2 = pPostStepPoint -> GetMaterial();    143   material2 = pPostStepPoint -> GetMaterial();
178                                                << 144   
179   theStatus = UndefinedSurf;                   << 
180                                                << 
181   const G4DynamicParticle* aParticle = aTrack.    145   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
182                                                << 146   
183   theParticleMomentum = aParticle->GetTotalMom    147   theParticleMomentum = aParticle->GetTotalMomentum();
184   previousMomentum = oldMomentum;                 148   previousMomentum = oldMomentum;
185   oldMomentum = aParticle->GetMomentumDirectio    149   oldMomentum = aParticle->GetMomentumDirection(); 
186                                                << 150     
187                                                << 151   //First case: not a boundary
188   // Fisrt case: not a boundary                << 152   if (pPostStepPoint->GetStepStatus() != fGeomBoundary || 
189   if (pPostStepPoint->GetStepStatus() != fGeom << 153       pPostStepPoint->GetPhysicalVolume() == pPreStepPoint->GetPhysicalVolume())
190       || pPostStepPoint->GetPhysicalVolume()-> << 154     {     
191   {                                            << 155       theStatus = NotAtBoundarySurf;
192     theStatus = NotAtBoundarySurf;             << 156       flag_franchissement_surface = false;
193     flag_franchissement_surface = false;       << 157       flag_reflexion = false;
194     flag_reflexion = false;                    << 158       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
195     return G4VDiscreteProcess::PostStepDoIt(aT << 159       
196   }                                            << 160     }
197   theStatus = UndefinedSurf;                      161   theStatus = UndefinedSurf;
198                                                << 162   
199   // Third case: same material                 << 163   //Third case: same material  
200   if (material1 == material2)                     164   if (material1 == material2)
201   {                                            << 165     {      
202     theStatus = SameMaterialSurf;              << 166       theStatus = SameMaterialSurf;
203     return G4VDiscreteProcess::PostStepDoIt(aT << 167       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
204   }                                            << 168     
205   if (verboseLevel > 3)                        << 169     }
206   {                                            << 170   if (verboseLevel > 0)
207     G4cout << G4endl << " Electron at Boundary << 171     {
208     G4VPhysicalVolume* thePrePV = pPreStepPoin << 172       G4cout << G4endl << " Electron at Boundary! " << G4endl;
209     G4VPhysicalVolume* thePostPV = pPostStepPo << 173       G4VPhysicalVolume* thePrePV = pPreStepPoint->GetPhysicalVolume();
210     if (thePrePV)  G4cout << " thePrePV:  " << << 174       G4VPhysicalVolume* thePostPV = pPostStepPoint->GetPhysicalVolume();
211     if (thePostPV) G4cout << " thePostPV: " << << 175       if (thePrePV)  G4cout << " thePrePV:  " << thePrePV->GetName() << G4endl;
212     G4cout << " Old Momentum Direction: " << o << 176       if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
213   }                                            << 177       G4cout << " Old Momentum Direction: " << oldMomentum << G4endl;
214                                                << 178     }
215   // Definition of the parameters for the surf << 179   
                                                   >> 180   //Definition of the parameters for the surface
216   G4ThreeVector theGlobalPoint = pPostStepPoin    181   G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
217                                                << 182   
218   G4Navigator* theNavigator = G4Transportation << 183   G4Navigator* theNavigator =
219     GetTransportationManager()->GetNavigatorFo << 184     G4TransportationManager::GetTransportationManager()->
220                                                << 185     GetNavigatorForTracking();
                                                   >> 186   
221   G4bool valid;                                   187   G4bool valid;
                                                   >> 188   
222   theGlobalNormal = theNavigator->GetGlobalExi    189   theGlobalNormal = theNavigator->GetGlobalExitNormal(theGlobalPoint, &valid);
223                                                << 190   //  G4cout << "Global exit normal = " << theGlobalNormal << " valid = " << valid << G4endl;
                                                   >> 191   
224   if (valid)                                      192   if (valid)
225   {                                            << 
226     theGlobalNormal = -theGlobalNormal;        << 
227   }                                            << 
228   else                                         << 
229   {                                            << 
230     G4ExceptionDescription ed;                 << 
231     ed << " G4MicroElecSurface/PostStepDoIt(): << 
232        << " The Navigator reports that it retu << 
233        << G4endl;                              << 
234     G4Exception("G4MuElecSurf::PostStepDoIt",  << 
235                 EventMustBeAborted, ed,        << 
236                 "Invalid Surface Normal - Geom << 
237   }                                            << 
238                                                << 
239   // Exception: the particle is not in the rig << 
240   if (oldMomentum * theGlobalNormal > 0.0)     << 
241   {                                            << 
242     theGlobalNormal = -theGlobalNormal;        << 
243   }                                            << 
244                                                << 
245   if (aTrack.GetStepLength()<=kCarTolerance/2  << 
246   {                                            << 
247     if (flag_reflexion == true)                << 
248     {                                          << 
249       flag_reflexion = false;                  << 
250       return G4VDiscreteProcess::PostStepDoIt( << 
251     }                                          << 
252                                                << 
253     theStatus = StepTooSmallSurf;              << 
254                                                << 
255     G4double energyThreshold_surface = 0.0*eV; << 
256                                                << 
257     WorkFunctionTable::iterator postStepWF;    << 
258     postStepWF = tableWF.find(pPostStepPoint-> << 
259     WorkFunctionTable::iterator preStepWF;     << 
260     preStepWF = tableWF.find(pPreStepPoint->Ge << 
261                                                << 
262     if (postStepWF == tableWF.end())           << 
263     {                                             193     {
264       G4String str = "Material ";              << 194       theGlobalNormal = -theGlobalNormal;
265       str += pPostStepPoint->GetMaterial()->Ge << 
266       G4Exception("G4Surface::G4Surface", "em0 << 
267       return nullptr;                          << 
268     }                                             195     }
269     else if (preStepWF == tableWF.end())       << 196   else
270     {                                             197     {
271       G4String str = "Material ";              << 198       G4ExceptionDescription ed;
272       str += pPreStepPoint->GetMaterial()->Get << 199       ed << " G4MicroElecSurface/PostStepDoIt(): "
273       G4Exception("G4Surface::G4Surface", "em0 << 200    << " The Navigator reports that it returned an invalid normal.\n"
274       return nullptr;                          << 201    << "PV: " << pPreStepPoint->GetPhysicalVolume()->GetName()
                                                   >> 202    << " TrackID= " << aTrack.GetTrackID()
                                                   >> 203    << " Ekin(MeV)= " << aTrack.GetKineticEnergy()
                                                   >> 204    << " position: " << theGlobalPoint
                                                   >> 205    << " direction: " << oldMomentum
                                                   >> 206    << G4endl;
                                                   >> 207       G4Exception("G4MuElecSurf::PostStepDoIt", "OpBoun01",
                                                   >> 208       FatalException, ed,
                                                   >> 209       "Invalid Surface Normal - Geometry must return valid surface normal");
                                                   >> 210       return 0;
275     }                                             211     }
276     else                                       << 212   
                                                   >> 213   //Exception: the particle is not in the right direction
                                                   >> 214   if (oldMomentum * theGlobalNormal > 0.0)
277     {                                             215     {
278       G4double thresholdNew_surface = postStep << 216       theGlobalNormal = -theGlobalNormal;
279       G4double thresholdOld_surface = preStepW << 
280       energyThreshold_surface = thresholdNew_s << 
281     }                                             217     }
                                                   >> 218   
                                                   >> 219   //Second case: step too small
                                                   >> 220   //Corrections bug rotation + réflexion
                                                   >> 221   if (aTrack.GetStepLength()<=kCarTolerance)
                                                   >> 222     { 
                                                   >> 223       theStatus = StepTooSmallSurf;       
                                                   >> 224   
                                                   >> 225       WorkFunctionTable::iterator postStepWF;
                                                   >> 226       postStepWF = tableWF.find(pPostStepPoint->GetMaterial()->GetName());
                                                   >> 227       WorkFunctionTable::iterator preStepWF;
                                                   >> 228       preStepWF = tableWF.find(pPreStepPoint->GetMaterial()->GetName());
                                                   >> 229   
                                                   >> 230       if (postStepWF == tableWF.end()) {
                                                   >> 231   G4String str = "Material ";
                                                   >> 232   str += pPostStepPoint->GetMaterial()->GetName() + " not found!";
                                                   >> 233   G4Exception("G4Surface::G4Surface", "em0002", FatalException, str);
                                                   >> 234   return 0;
                                                   >> 235       }
                                                   >> 236   
                                                   >> 237       else if (preStepWF == tableWF.end()) {
                                                   >> 238   G4String str = "Material ";
                                                   >> 239   str += pPreStepPoint->GetMaterial()->GetName() + " not found!";
                                                   >> 240   G4Exception("G4Surface::G4Surface", "em0002", FatalException, str);
                                                   >> 241   return 0;
                                                   >> 242       }
                                                   >> 243       
                                                   >> 244       if (pPreStepPoint->GetMaterial() != pPostStepPoint->GetMaterial()) {
                                                   >> 245   
                                                   >> 246   flag_franchissement_surface = false;
282                                                   247 
283     if (flag_franchissement_surface == true)   << 248   if (flag_reflexion == true && flag_normal == true) {
284     {                                          << 249     aParticleChange.ProposeMomentumDirection(-Reflexion(aStep.GetPostStepPoint()));
285       aParticleChange.ProposeEnergy(aStep.GetP << 250     flag_reflexion = false;
286       flag_franchissement_surface = false;     << 251     flag_normal = false;
287     }                                          << 252   }
288     if (flag_reflexion == true && flag_normal  << 253       }
289     {                                          << 254       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
290       aParticleChange.ProposeMomentumDirection << 
291       flag_reflexion = false;                  << 
292       flag_normal = false;                     << 
293     }                                             255     }
294     return G4VDiscreteProcess::PostStepDoIt(aT << 
295   }                                            << 
296                                                << 
297   flag_normal = (theGlobalNormal == G4ThreeVec << 
298               || theGlobalNormal == G4ThreeVec << 
299                                                   256 
                                                   >> 257   flag_normal = (theGlobalNormal.x() == 0.0 && theGlobalNormal.y() == 0.0);
                                                   >> 258   
300   G4LogicalSurface* Surface = nullptr;            259   G4LogicalSurface* Surface = nullptr;
301                                                << 260   
302   Surface = G4LogicalBorderSurface::GetSurface    261   Surface = G4LogicalBorderSurface::GetSurface
303                   (pPreStepPoint ->GetPhysical << 262     (pPreStepPoint ->GetPhysicalVolume(),
304                    pPostStepPoint->GetPhysical << 263      pPostStepPoint->GetPhysicalVolume());
305                                                << 264   
306   if (Surface == nullptr)                         265   if (Surface == nullptr)
307   {                                            << 
308     G4bool enteredDaughter=(pPostStepPoint->Ge << 
309                          == pPreStepPoint->Get << 
310     if(enteredDaughter)                        << 
311     {                                          << 
312       Surface = G4LogicalSkinSurface::GetSurfa << 
313                 (pPostStepPoint->GetPhysicalVo << 
314                                                << 
315       if(Surface == nullptr)                   << 
316       Surface = G4LogicalSkinSurface::GetSurfa << 
317                 (pPreStepPoint->GetPhysicalVol << 
318     }                                          << 
319     else                                       << 
320     {                                             266     {
321       Surface = G4LogicalSkinSurface::GetSurfa << 267       G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
322                 (pPreStepPoint->GetPhysicalVol << 268             ->GetMotherLogical() ==
323                                                << 269             pPreStepPoint->GetPhysicalVolume()
324       if(Surface == nullptr)                   << 270             ->GetLogicalVolume());
325       Surface = G4LogicalSkinSurface::GetSurfa << 271       if(enteredDaughter)
326                 (pPostStepPoint->GetPhysicalVo << 272   {
                                                   >> 273     Surface = G4LogicalSkinSurface::GetSurface
                                                   >> 274       (pPostStepPoint->GetPhysicalVolume()->
                                                   >> 275        GetLogicalVolume());
                                                   >> 276     
                                                   >> 277     if(Surface == nullptr)
                                                   >> 278       Surface = G4LogicalSkinSurface::GetSurface
                                                   >> 279         (pPreStepPoint->GetPhysicalVolume()->
                                                   >> 280          GetLogicalVolume());
                                                   >> 281   }
                                                   >> 282       else 
                                                   >> 283   {
                                                   >> 284     Surface = G4LogicalSkinSurface::GetSurface
                                                   >> 285       (pPreStepPoint->GetPhysicalVolume()->
                                                   >> 286        GetLogicalVolume());
                                                   >> 287     
                                                   >> 288     if(Surface == nullptr)
                                                   >> 289       Surface = G4LogicalSkinSurface::GetSurface
                                                   >> 290         (pPostStepPoint->GetPhysicalVolume()->
                                                   >> 291          GetLogicalVolume());
                                                   >> 292   }
327     }                                             293     }
328   }                                            << 294     
329                                                << 
330   // Definition of the parameters for the surf << 
331   G4VPhysicalVolume* thePrePV = pPreStepPoint-    295   G4VPhysicalVolume* thePrePV = pPreStepPoint->GetPhysicalVolume();
332   G4VPhysicalVolume* thePostPV = pPostStepPoin    296   G4VPhysicalVolume* thePostPV = pPostStepPoint->GetPhysicalVolume();
333                                                << 297   
334   energyThreshold = 0.0*eV;                    << 298   if (thePostPV)
335   G4double energyDelta = 0;                    << 
336                                                << 
337   if ((thePrePV)&&(thePostPV))                 << 
338   {                                            << 
339     WorkFunctionTable::iterator postStepWF;    << 
340     postStepWF = tableWF.find(thePostPV->GetLo << 
341     WorkFunctionTable::iterator preStepWF;     << 
342     preStepWF = tableWF.find(thePrePV->GetLogi << 
343                                                << 
344     if (postStepWF == tableWF.end())           << 
345     {                                          << 
346       G4String str = "Material ";              << 
347       str += thePostPV->GetLogicalVolume()->Ge << 
348       G4Exception("G4Surface::G4Surface", "em0 << 
349       return nullptr;                          << 
350     }                                          << 
351     else if (preStepWF == tableWF.end())       << 
352     {                                          << 
353       G4String str = "Material ";              << 
354       str += thePrePV->GetLogicalVolume()->Get << 
355       G4Exception("G4Surface::G4Surface", "em0 << 
356       return nullptr;                          << 
357     }                                          << 
358     else                                       << 
359     {                                             299     {
360       G4double thresholdNew = postStepWF->seco << 300       WorkFunctionTable::iterator postStepWF;
361       G4double thresholdOld = preStepWF->secon << 301       postStepWF = tableWF.find(thePostPV->GetLogicalVolume()->GetMaterial()->GetName());
362                                                << 302       WorkFunctionTable::iterator preStepWF;
363       energyThreshold = thresholdNew - thresho << 303       preStepWF = tableWF.find(thePrePV->GetLogicalVolume()->GetMaterial()->GetName());
364       energyDelta = thresholdOld- thresholdNew << 304       
                                                   >> 305       if (postStepWF == tableWF.end()) {
                                                   >> 306   G4String str = "Material ";
                                                   >> 307   str += thePostPV->GetLogicalVolume()->GetMaterial()->GetName() + " not found!";
                                                   >> 308   G4Exception("G4Surface::G4Surface", "em0002", FatalException, str);
                                                   >> 309   return 0;
                                                   >> 310       }
                                                   >> 311       else if (preStepWF == tableWF.end()) {
                                                   >> 312   G4String str = "Material ";
                                                   >> 313   str += thePrePV->GetLogicalVolume()->GetMaterial()->GetName() + " not found!";
                                                   >> 314   G4Exception("G4Surface::G4Surface", "em0002", FatalException, str);
                                                   >> 315   return 0;
                                                   >> 316       }
                                                   >> 317       else
                                                   >> 318   {
                                                   >> 319     G4double thresholdNew = postStepWF->second;
                                                   >> 320     G4double thresholdOld = preStepWF->second;
                                                   >> 321     energyThreshold = thresholdNew - thresholdOld;
                                                   >> 322   }
365     }                                             323     }
366   }                                            << 324   
367                                                << 325   ekint = pPreStepPoint->GetKineticEnergy();
368   ekint=aStep.GetPostStepPoint()->GetKineticEn << 
369   thetat= GetIncidentAngle(); //angle d'incide    326   thetat= GetIncidentAngle(); //angle d'incidence
370   G4double ekinNormalt=ekint*std::cos(thetat)*    327   G4double ekinNormalt=ekint*std::cos(thetat)*std::cos(thetat); 
371                                                << 328   
372   thetaft=std::asin(std::sqrt(ekint/(ekint+ene << 329   G4double atet = std::sqrt(ekint/(ekint+energyThreshold))*std::sin(thetat);
373   if(std::sqrt(ekint/(ekint+energyThreshold))* << 330   thetaft = (atet > 1.0) ? pi*0.5 : std::asin(atet);//Angle de réfraction
374   {                                            << 331   
375     thetaft=std::asin(1.0);                    << 
376   }                                            << 
377                                                << 
378   G4double aleat=G4UniformRand();                 332   G4double aleat=G4UniformRand();   
379                                                << 333   
380   G4double waveVectort=std::sqrt(2*9.1093826E- << 334   const G4double waveVectort=std::sqrt(2*9.1093826E-31*1.602176487E-19)/(6.6260755E-34/(2.0*pi));
381                                                << 335   
382   // Parameter for an exponential barrier of p << 336   //Parameter for an exponential barrier of potential (Thèse P68)
383   G4double at=0.5E-10;                         << 337   const G4double at=0.5E-10;
384                                                << 338   
                                                   >> 339   //G4double modif already declared in .hh
385   crossingProbability=0;                          340   crossingProbability=0;
386                                                << 341   
387   G4double kft=waveVectort*std::sqrt(ekint+ene    342   G4double kft=waveVectort*std::sqrt(ekint+energyThreshold)*std::cos(thetaft);
388   G4double kit=waveVectort*std::sqrt(ekinNorma    343   G4double kit=waveVectort*std::sqrt(ekinNormalt); 
389                                                << 344   
390   crossingProbability=1-(std::pow(std::sinh(pi << 345   G4double yy = std::sinh(pi*at*(kit-kft))/std::sinh(pi*at*(kit+kft));
391                                                << 346   crossingProbability = 1 - yy*yy;
392   // First case: the electron crosses the surf << 347   
393   if((aleat<=crossingProbability)&&(ekint> ene << 348   //First case: the electron crosses the surface
394   {                                            << 349   if((aleat<=crossingProbability)&&(ekint>std::abs(energyThreshold)))
395     if (aStep.GetPreStepPoint()->GetMaterial() << 
396      != aStep.GetPostStepPoint()->GetMaterial( << 
397     {                                             350     {
398       aParticleChange.ProposeEnergy(ekint - en << 351       if (pPreStepPoint->GetMaterial() != pPostStepPoint->GetMaterial()) {
399       flag_franchissement_surface = true;      << 352   flag_franchissement_surface = true;
400     }                                          << 353       }
401                                                << 354       
402     thetaft=std::abs(thetaft-thetat);          << 355       thetaft=std::abs(thetaft-thetat);
403                                                << 356       
404     G4ThreeVector zVerst = aStep.GetPostStepPo << 357       G4ThreeVector zVerst = aStep.GetPostStepPoint()->GetMomentumDirection();
405     G4ThreeVector xVerst = zVerst.orthogonal() << 358       G4ThreeVector xVerst = zVerst.orthogonal();
406     G4ThreeVector yVerst = zVerst.cross(xVerst << 359       G4ThreeVector yVerst = zVerst.cross(xVerst);
407                                                << 360 
408     G4double xDirt = std::sqrt(1. - std::cos(t << 361       G4double cost  = std::cos(thetaft); 
409     G4double yDirt = xDirt;                    << 362       G4double xDirt = std::sqrt(1. - cost*cost);
410                                                << 363       G4double yDirt = xDirt;
411     G4ThreeVector zPrimeVerst=((xDirt*xVerst + << 364       
                                                   >> 365       G4ThreeVector zPrimeVerst = xDirt*xVerst + yDirt*yVerst + cost*zVerst;
412                                                   366     
413     aParticleChange.ProposeMomentumDirection(z << 367       aParticleChange.ProposeMomentumDirection(zPrimeVerst.unit());
414   }                                            << 
415   else if ((aleat > crossingProbability) && (e << 
416   {                                            << 
417     flag_reflexion = true;                     << 
418     if (flag_normal)                           << 
419     {                                          << 
420       aParticleChange.ProposeMomentumDirection << 
421     }                                          << 
422     else                                       << 
423     {                                          << 
424       aParticleChange.ProposeMomentumDirection << 
425     }                                          << 
426   }                                            << 
427   else                                         << 
428   {                                            << 
429     if (flag_normal)                           << 
430     {                                          << 
431       aParticleChange.ProposeMomentumDirection << 
432     }                                             368     }
433     else                                       << 369   else if ((aleat > crossingProbability) && (ekint>std::abs(energyThreshold)))
434     {                                             370     {
435       aParticleChange.ProposeMomentumDirection << 371       flag_reflexion = true;
                                                   >> 372       if (flag_normal) { aParticleChange.ProposeMomentumDirection(-oldMomentum.unit()); }
                                                   >> 373       else { aParticleChange.ProposeMomentumDirection(Reflexion(aStep.GetPostStepPoint())); }
                                                   >> 374       
436     }                                             375     }
                                                   >> 376   else {
                                                   >> 377     
                                                   >> 378     if (flag_normal) { aParticleChange.ProposeMomentumDirection(-oldMomentum.unit()); }
                                                   >> 379     else { aParticleChange.ProposeMomentumDirection(Reflexion(aStep.GetPostStepPoint())); }
437     flag_reflexion = true;                        380     flag_reflexion = true;
438   }                                               381   }
                                                   >> 382   
439   return G4VDiscreteProcess::PostStepDoIt(aTra    383   return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
440 }                                                 384 }
441                                                   385 
442 //....oooOO0OOooo........oooOO0OOooo........oo    386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
443                                                   387 
444 G4double G4MicroElecSurface::GetMeanFreePath(c    388 G4double G4MicroElecSurface::GetMeanFreePath(const G4Track&, G4double,
445                                              G << 389                                                  G4ForceCondition* condition)
446 {                                                 390 { 
447   *condition = Forced;                            391   *condition = Forced;
448   return  DBL_MAX;                                392   return  DBL_MAX;   
449 }                                                 393 }
450                                                   394 
451 //....oooOO0OOooo........oooOO0OOooo........oo    395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
452                                                   396 
453 G4double G4MicroElecSurface::GetIncidentAngle(    397 G4double G4MicroElecSurface::GetIncidentAngle() 
454 {                                                 398 {
455   theFacetNormal=theGlobalNormal;                 399   theFacetNormal=theGlobalNormal; 
456                                                   400   
457   G4double PdotN = oldMomentum * theFacetNorma    401   G4double PdotN = oldMomentum * theFacetNormal;
458   G4double magP= oldMomentum.mag();               402   G4double magP= oldMomentum.mag();
459   G4double magN= theFacetNormal.mag();            403   G4double magN= theFacetNormal.mag();
460                                                   404   
461   G4double incidentangle = pi - std::acos(Pdot    405   G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
462                                                   406   
463   return incidentangle;                           407   return incidentangle;
464 }                                                 408 }
465                                                   409 
466 //....oooOO0OOooo........oooOO0OOooo........oo    410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
467                                                   411 
468 G4ThreeVector G4MicroElecSurface::Reflexion(co    412 G4ThreeVector G4MicroElecSurface::Reflexion(const G4StepPoint* PostStepPoint)
469 {                                                 413 {
470   // Normale                                   << 414   //Normale
471   G4double Nx = theGlobalNormal.x();              415   G4double Nx = theGlobalNormal.x();
472   G4double Ny = theGlobalNormal.y();              416   G4double Ny = theGlobalNormal.y();
473   G4double Nz = theGlobalNormal.z();              417   G4double Nz = theGlobalNormal.z();
474                                                   418   
475   // PostStepPoint                             << 419   //PostStepPoint
476   G4double PSx = PostStepPoint->GetPosition().    420   G4double PSx = PostStepPoint->GetPosition().x();
477   G4double PSy = PostStepPoint->GetPosition().    421   G4double PSy = PostStepPoint->GetPosition().y();
478   G4double PSz = PostStepPoint->GetPosition().    422   G4double PSz = PostStepPoint->GetPosition().z();
479                                                   423   
480   // P(alpha,beta,gamma) - PostStep avec trans << 424   //P(alpha,beta,gamma) - PostStep avec translation momentum
481   G4double alpha = PSx + oldMomentum.x();         425   G4double alpha = PSx + oldMomentum.x();
482   G4double beta = PSy + oldMomentum.y();          426   G4double beta = PSy + oldMomentum.y();
483   G4double gamma = PSz + oldMomentum.z();         427   G4double gamma = PSz + oldMomentum.z();
484   G4double r = theGlobalNormal.mag();             428   G4double r = theGlobalNormal.mag();
485   G4double x, y, z, d, A, B, PM2x, PM2y, PM2z;    429   G4double x, y, z, d, A, B, PM2x, PM2y, PM2z;
486   d = -(Nx*PSx + Ny*PSy + Nz*PSz);                430   d = -(Nx*PSx + Ny*PSy + Nz*PSz);
487                                                   431   
488   if (Ny == 0 && Nx == 0)                      << 432   if (Ny == 0 && Nx == 0) {
489   {                                            << 
490     gamma = -gamma;                               433     gamma = -gamma;
491   }                                               434   }  
492   else                                         << 435   else {
493   {                                            << 436     if (Ny == 0) {
494     if (Ny == 0)                               << 437       
495     {                                          << 
496       A = (Nz*Nz*alpha) + (Nx*Nx*PSx) + (Nx*Nz    438       A = (Nz*Nz*alpha) + (Nx*Nx*PSx) + (Nx*Nz*(PSz - gamma));
497       B = r*r;                                    439       B = r*r;
498                                                   440       
499       // M(x,y,z) - Projection de P sur la sur << 441       //M(x,y,z) - Projection de P sur la surface
500       x = A / B;                                  442       x = A / B;
501       y = beta;                                   443       y = beta;
502       z = (x - alpha)*(Nz / Nx) + gamma;          444       z = (x - alpha)*(Nz / Nx) + gamma;
503     }                                             445     }
504     else                                       << 446     
505     {                                          << 447     else {
                                                   >> 448       
506       A = (r*r) / Ny;                             449       A = (r*r) / Ny;
507       B = (beta / Ny)*(Nx*Nx + Nz*Nz) - (Nx*al    450       B = (beta / Ny)*(Nx*Nx + Nz*Nz) - (Nx*alpha + Nz*gamma + d);
508                                                   451       
509       //M(x,y,z) - Projection de P sur la surf    452       //M(x,y,z) - Projection de P sur la surface
510       y = B / A;                                  453       y = B / A;
511       x = (y - beta)*(Nx / Ny) + alpha;           454       x = (y - beta)*(Nx / Ny) + alpha;
512       z = (y - beta)*(Nz / Ny) + gamma;           455       z = (y - beta)*(Nz / Ny) + gamma;
513     }                                             456     }
514                                                   457     
515     // Vecteur 2*PM                            << 458     //Vecteur 2*PM
516     PM2x = 2 * (x - alpha);  PM2y = 2 * (y - b    459     PM2x = 2 * (x - alpha);  PM2y = 2 * (y - beta);  PM2z = 2 * (z - gamma);
517                                                   460     
518     // Nouveau point P                         << 461     //Nouveau point P
519     alpha += PM2x; beta += PM2y; gamma += PM2z    462     alpha += PM2x; beta += PM2y; gamma += PM2z;
520                                                   463     
521   }                                               464   }
522                                                   465   
523   G4ThreeVector newMomentum = G4ThreeVector(al    466   G4ThreeVector newMomentum = G4ThreeVector(alpha-PSx,beta-PSy,gamma-PSz);  
524   return newMomentum.unit();                      467   return newMomentum.unit();
525 }                                                 468 }
526                                                   469 
527 //....oooOO0OOooo........oooOO0OOooo........oo    470 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
528                                                   471 
529 G4MicroElecSurfaceStatus G4MicroElecSurface::G    472 G4MicroElecSurfaceStatus G4MicroElecSurface::GetStatus() const 
530 {                                                 473 { 
531   return theStatus;                               474   return theStatus; 
532 }                                                 475 } 
533                                                   476 
534 //....oooOO0OOooo........oooOO0OOooo........oo    477 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
535                                                   478 
536                                                   479   
537 void G4MicroElecSurface::SetFlagFranchissement    480 void G4MicroElecSurface::SetFlagFranchissement() 
538 {                                                 481 { 
539   flag_franchissement_surface = false;            482   flag_franchissement_surface = false; 
540 }                                                 483 } 
541                                                   484  
542 //....oooOO0OOooo........oooOO0OOooo........oo    485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
543                                                   486 
544                                                   487