| Geant4 Cross Reference |
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 = 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 { 193 {
249 flag_reflexion = false; << 194 theGlobalNormal = -theGlobalNormal;
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 { <<
264 G4String str = "Material "; <<
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 { 266 {
312 Surface = G4LogicalSkinSurface::GetSurfa << 267 G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
313 (pPostStepPoint->GetPhysicalVo << 268 ->GetMotherLogical() ==
314 << 269 pPreStepPoint->GetPhysicalVolume()
315 if(Surface == nullptr) << 270 ->GetLogicalVolume());
316 Surface = G4LogicalSkinSurface::GetSurfa << 271 if(enteredDaughter)
317 (pPreStepPoint->GetPhysicalVol << 272 {
318 } << 273 Surface = G4LogicalSkinSurface::GetSurface
319 else << 274 (pPostStepPoint->GetPhysicalVolume()->
320 { << 275 GetLogicalVolume());
321 Surface = G4LogicalSkinSurface::GetSurfa << 276
322 (pPreStepPoint->GetPhysicalVol << 277 if(Surface == nullptr)
323 << 278 Surface = G4LogicalSkinSurface::GetSurface
324 if(Surface == nullptr) << 279 (pPreStepPoint->GetPhysicalVolume()->
325 Surface = G4LogicalSkinSurface::GetSurfa << 280 GetLogicalVolume());
326 (pPostStepPoint->GetPhysicalVo << 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