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