| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // G4MicroElecSurface.cc,
28 // 2020/05/20 P. Caron, C. Ingui
29 // Q. Gibaru is with
30 // M. Raine and D. La
31 //
32 // A part of this work has been funded by the
33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France
34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 T
35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CED
36 //
37 // Based on the following publications
38 //
39 // - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine,
40 // Geant4 physics processes for microdosimet
41 // Extension of MicroElec to very low energi
42 // NIM B, 2020, in review.
43 //
44 // - Modele de transport d'electrons a basse e
45 // applications spatiales (OSMOSEE, GEANT4),
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........oo
58
59 G4MicroElecSurface::G4MicroElecSurface(const G
60 : G4VDiscreteProcess(processName, type),
61 oldMomentum(0.,0.,0.), previousMomentum(0.
62 theGlobalNormal(0.,0.,0.), theFacetNormal(
63 {
64 if ( verboseLevel > 0)
65 {
66 G4cout << GetProcessName() << " is created
67 }
68
69 isInitialised=false;
70 SetProcessSubType(fSurfaceReflection);
71
72 theStatus = UndefinedSurf;
73 material1 = nullptr;
74 material2 = nullptr;
75
76 kCarTolerance = G4GeometryTolerance::GetInst
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 =
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oo
88
89 G4MicroElecSurface::~G4MicroElecSurface()
90 {}
91
92 //....oooOO0OOooo........oooOO0OOooo........oo
93
94 G4bool
95 G4MicroElecSurface::IsApplicable(const G4Parti
96 {
97 return ( aParticleType.GetPDGEncoding() == 1
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101
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
133 {
134 if (isInitialised) { return; }
135
136 G4ProductionCutsTable* theCoupleTable =
137 G4ProductionCutsTable::GetProductionCutsTa
138 G4int numOfCouples = (G4int)theCoupleTable->
139 G4cout << "G4MicroElecSurface::Initialise: N
140 << numOfCouples << G4endl;
141
142 for (G4int i = 0; i < numOfCouples; ++i)
143 {
144 const G4Material* material =
145 theCoupleTable->GetMaterialCutsCouple(i)
146
147 G4cout << "G4Surface, Material " << i + 1
148 << numOfCouples << " : " << materia
149 if (material->GetName() == "Vacuum")
150 {
151 tableWF[material->GetName()] = 0;
152 continue;
153 }
154 G4String mat = material->GetName();
155 G4MicroElecMaterialStructure str = G4Micro
156 tableWF[mat] = str.GetWorkFunction();
157 }
158 isInitialised = true;
159 }
160
161 //....oooOO0OOooo........oooOO0OOooo........oo
162
163 G4VParticleChange* G4MicroElecSurface::PostSte
164 {
165
166 if (!isInitialised) { Initialise(); }
167 theStatus = UndefinedSurf;
168
169 // Definition of the parameters for the part
170 aParticleChange.Initialize(aTrack);
171 aParticleChange.ProposeVelocity(aTrack.GetVe
172
173 G4StepPoint* pPreStepPoint = aStep.GetPreSt
174 G4StepPoint* pPostStepPoint = aStep.GetPostS
175
176 material1 = pPreStepPoint -> GetMaterial();
177 material2 = pPostStepPoint -> GetMaterial();
178
179 theStatus = UndefinedSurf;
180
181 const G4DynamicParticle* aParticle = aTrack.
182
183 theParticleMomentum = aParticle->GetTotalMom
184 previousMomentum = oldMomentum;
185 oldMomentum = aParticle->GetMomentumDirectio
186
187
188 // Fisrt case: not a boundary
189 if (pPostStepPoint->GetStepStatus() != fGeom
190 || pPostStepPoint->GetPhysicalVolume()->
191 {
192 theStatus = NotAtBoundarySurf;
193 flag_franchissement_surface = false;
194 flag_reflexion = false;
195 return G4VDiscreteProcess::PostStepDoIt(aT
196 }
197 theStatus = UndefinedSurf;
198
199 // Third case: same material
200 if (material1 == material2)
201 {
202 theStatus = SameMaterialSurf;
203 return G4VDiscreteProcess::PostStepDoIt(aT
204 }
205 if (verboseLevel > 3)
206 {
207 G4cout << G4endl << " Electron at Boundary
208 G4VPhysicalVolume* thePrePV = pPreStepPoin
209 G4VPhysicalVolume* thePostPV = pPostStepPo
210 if (thePrePV) G4cout << " thePrePV: " <<
211 if (thePostPV) G4cout << " thePostPV: " <<
212 G4cout << " Old Momentum Direction: " << o
213 }
214
215 // Definition of the parameters for the surf
216 G4ThreeVector theGlobalPoint = pPostStepPoin
217
218 G4Navigator* theNavigator = G4Transportation
219 GetTransportationManager()->GetNavigatorFo
220
221 G4bool valid;
222 theGlobalNormal = theNavigator->GetGlobalExi
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 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 {
264 G4String str = "Material ";
265 str += pPostStepPoint->GetMaterial()->Ge
266 G4Exception("G4Surface::G4Surface", "em0
267 return nullptr;
268 }
269 else if (preStepWF == tableWF.end())
270 {
271 G4String str = "Material ";
272 str += pPreStepPoint->GetMaterial()->Get
273 G4Exception("G4Surface::G4Surface", "em0
274 return nullptr;
275 }
276 else
277 {
278 G4double thresholdNew_surface = postStep
279 G4double thresholdOld_surface = preStepW
280 energyThreshold_surface = thresholdNew_s
281 }
282
283 if (flag_franchissement_surface == true)
284 {
285 aParticleChange.ProposeEnergy(aStep.GetP
286 flag_franchissement_surface = false;
287 }
288 if (flag_reflexion == true && flag_normal
289 {
290 aParticleChange.ProposeMomentumDirection
291 flag_reflexion = false;
292 flag_normal = false;
293 }
294 return G4VDiscreteProcess::PostStepDoIt(aT
295 }
296
297 flag_normal = (theGlobalNormal == G4ThreeVec
298 || theGlobalNormal == G4ThreeVec
299
300 G4LogicalSurface* Surface = nullptr;
301
302 Surface = G4LogicalBorderSurface::GetSurface
303 (pPreStepPoint ->GetPhysical
304 pPostStepPoint->GetPhysical
305
306 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 {
321 Surface = G4LogicalSkinSurface::GetSurfa
322 (pPreStepPoint->GetPhysicalVol
323
324 if(Surface == nullptr)
325 Surface = G4LogicalSkinSurface::GetSurfa
326 (pPostStepPoint->GetPhysicalVo
327 }
328 }
329
330 // Definition of the parameters for the surf
331 G4VPhysicalVolume* thePrePV = pPreStepPoint-
332 G4VPhysicalVolume* thePostPV = pPostStepPoin
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->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 {
360 G4double thresholdNew = postStepWF->seco
361 G4double thresholdOld = preStepWF->secon
362
363 energyThreshold = thresholdNew - thresho
364 energyDelta = thresholdOld- thresholdNew
365 }
366 }
367
368 ekint=aStep.GetPostStepPoint()->GetKineticEn
369 thetat= GetIncidentAngle(); //angle d'incide
370 G4double ekinNormalt=ekint*std::cos(thetat)*
371
372 thetaft=std::asin(std::sqrt(ekint/(ekint+ene
373 if(std::sqrt(ekint/(ekint+energyThreshold))*
374 {
375 thetaft=std::asin(1.0);
376 }
377
378 G4double aleat=G4UniformRand();
379
380 G4double waveVectort=std::sqrt(2*9.1093826E-
381
382 // Parameter for an exponential barrier of p
383 G4double at=0.5E-10;
384
385 crossingProbability=0;
386
387 G4double kft=waveVectort*std::sqrt(ekint+ene
388 G4double kit=waveVectort*std::sqrt(ekinNorma
389
390 crossingProbability=1-(std::pow(std::sinh(pi
391
392 // First case: the electron crosses the surf
393 if((aleat<=crossingProbability)&&(ekint> ene
394 {
395 if (aStep.GetPreStepPoint()->GetMaterial()
396 != aStep.GetPostStepPoint()->GetMaterial(
397 {
398 aParticleChange.ProposeEnergy(ekint - en
399 flag_franchissement_surface = true;
400 }
401
402 thetaft=std::abs(thetaft-thetat);
403
404 G4ThreeVector zVerst = aStep.GetPostStepPo
405 G4ThreeVector xVerst = zVerst.orthogonal()
406 G4ThreeVector yVerst = zVerst.cross(xVerst
407
408 G4double xDirt = std::sqrt(1. - std::cos(t
409 G4double yDirt = xDirt;
410
411 G4ThreeVector zPrimeVerst=((xDirt*xVerst +
412
413 aParticleChange.ProposeMomentumDirection(z
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 }
433 else
434 {
435 aParticleChange.ProposeMomentumDirection
436 }
437 flag_reflexion = true;
438 }
439 return G4VDiscreteProcess::PostStepDoIt(aTra
440 }
441
442 //....oooOO0OOooo........oooOO0OOooo........oo
443
444 G4double G4MicroElecSurface::GetMeanFreePath(c
445 G
446 {
447 *condition = Forced;
448 return DBL_MAX;
449 }
450
451 //....oooOO0OOooo........oooOO0OOooo........oo
452
453 G4double G4MicroElecSurface::GetIncidentAngle(
454 {
455 theFacetNormal=theGlobalNormal;
456
457 G4double PdotN = oldMomentum * theFacetNorma
458 G4double magP= oldMomentum.mag();
459 G4double magN= theFacetNormal.mag();
460
461 G4double incidentangle = pi - std::acos(Pdot
462
463 return incidentangle;
464 }
465
466 //....oooOO0OOooo........oooOO0OOooo........oo
467
468 G4ThreeVector G4MicroElecSurface::Reflexion(co
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().
477 G4double PSy = PostStepPoint->GetPosition().
478 G4double PSz = PostStepPoint->GetPosition().
479
480 // P(alpha,beta,gamma) - PostStep avec trans
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
497 B = r*r;
498
499 // M(x,y,z) - Projection de P sur la sur
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*al
508
509 //M(x,y,z) - Projection de P sur la surf
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 - b
517
518 // Nouveau point P
519 alpha += PM2x; beta += PM2y; gamma += PM2z
520
521 }
522
523 G4ThreeVector newMomentum = G4ThreeVector(al
524 return newMomentum.unit();
525 }
526
527 //....oooOO0OOooo........oooOO0OOooo........oo
528
529 G4MicroElecSurfaceStatus G4MicroElecSurface::G
530 {
531 return theStatus;
532 }
533
534 //....oooOO0OOooo........oooOO0OOooo........oo
535
536
537 void G4MicroElecSurface::SetFlagFranchissement
538 {
539 flag_franchissement_surface = false;
540 }
541
542 //....oooOO0OOooo........oooOO0OOooo........oo
543
544