| 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 // Optical Photon Boundary Process Class Imple
28 //////////////////////////////////////////////
29 //
30 // File: G4OpBoundaryProcess.cc
31 // Description: Discrete Process -- reflection
32 // optical in
33 // Version: 1.1
34 // Created: 1997-06-18
35 // Modified: 1998-05-25 - Correct parallel
36 // (thanks to: Stefa
37 // 1998-05-28 - NULL Rindex point
38 // (thanks to: Stefa
39 // 1998-06-11 - delete *sint1 in
40 // (thanks to: Giova
41 // 1998-06-19 - move from GetLoca
42 // method: GetLocalE
43 // the surface norma
44 // 1998-11-07 - NULL OpticalSurfa
45 // comparison not sh
46 // remove sin1, sin2
47 // (thanks to Stefan
48 // 1999-10-10 - Accommodate chang
49 // changing logic in
50 // 2001-10-18 - avoid Linux (gcc-
51 // might be used uni
52 // moved E2_perp, E2
53 // 2003-11-27 - Modified line 168
54 // G4OpticalSurface
55 // 2004-02-02 - Set theStatus = U
56 // 2005-07-28 - add G4ProcessType
57 // 2006-11-04 - add capability of
58 // off a metal surfa
59 // of refraction - T
60 // Hauptman (Dept. o
61 // 2009-11-10 - add capability of
62 // with Look-Up-Tabl
63 // optical reflectan
64 // treatments - Than
65 // William Moses (La
66 // 2013-06-01 - add the capabilit
67 // of a dichronic fi
68 // 2017-02-24 - add capability of
69 // with Look-Up-Tabl
70 //
71 // Author: Peter Gumplinger
72 // adopted from work by Werner Keil - April
73 //
74 //////////////////////////////////////////////
75
76 #include "G4OpBoundaryProcess.hh"
77
78 #include "G4ios.hh"
79 #include "G4GeometryTolerance.hh"
80 #include "G4LogicalBorderSurface.hh"
81 #include "G4LogicalSkinSurface.hh"
82 #include "G4OpProcessSubType.hh"
83 #include "G4OpticalParameters.hh"
84 #include "G4ParallelWorldProcess.hh"
85 #include "G4PhysicalConstants.hh"
86 #include "G4SystemOfUnits.hh"
87 #include "G4TransportationManager.hh"
88 #include "G4VSensitiveDetector.hh"
89
90 //....oooOO0OOooo........oooOO0OOooo........oo
91 G4OpBoundaryProcess::G4OpBoundaryProcess(const
92 G4Pro
93 : G4VDiscreteProcess(processName, ptype)
94 {
95 Initialise();
96
97 if(verboseLevel > 0)
98 {
99 G4cout << GetProcessName() << " is created
100 }
101 SetProcessSubType(fOpBoundary);
102
103 fStatus = Undefined;
104 fModel = glisur;
105 fFinish = polished;
106 fReflectivity = 1.;
107 fEfficiency = 0.;
108 fTransmittance = 0.;
109 fSurfaceRoughness = 0.;
110 fProb_sl = 0.;
111 fProb_ss = 0.;
112 fProb_bs = 0.;
113
114 fRealRIndexMPV = nullptr;
115 fImagRIndexMPV = nullptr;
116 fMaterial1 = nullptr;
117 fMaterial2 = nullptr;
118 fOpticalSurface = nullptr;
119 fCarTolerance = G4GeometryTolerance::GetIn
120
121 f_iTE = f_iTM = 0;
122 fPhotonMomentum = 0.;
123 fRindex1 = fRindex2 = 1.;
124 fSint1 = 0.;
125 fDichroicVector = nullptr;
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oo
129 G4OpBoundaryProcess::~G4OpBoundaryProcess() =
130
131 //....oooOO0OOooo........oooOO0OOooo........oo
132 void G4OpBoundaryProcess::PreparePhysicsTable(
133 {
134 Initialise();
135 }
136
137 //....oooOO0OOooo........oooOO0OOooo........oo
138 void G4OpBoundaryProcess::Initialise()
139 {
140 G4OpticalParameters* params = G4OpticalParam
141 SetInvokeSD(params->GetBoundaryInvokeSD());
142 SetVerboseLevel(params->GetBoundaryVerboseLe
143 }
144
145 //....oooOO0OOooo........oooOO0OOooo........oo
146 G4VParticleChange* G4OpBoundaryProcess::PostSt
147
148 {
149 fStatus = Undefined;
150 aParticleChange.Initialize(aTrack);
151 aParticleChange.ProposeVelocity(aTrack.GetVe
152
153 // Get hyperStep from G4ParallelWorldProces
154 // NOTE: PostSetpDoIt of this process to be
155 // G4ParallelWorldProcess!
156 const G4Step* pStep = &aStep;
157 const G4Step* hStep = G4ParallelWorldProcess
158 if(hStep != nullptr)
159 pStep = hStep;
160
161 if(pStep->GetPostStepPoint()->GetStepStatus(
162 {
163 fMaterial1 = pStep->GetPreStepPoint()->Get
164 fMaterial2 = pStep->GetPostStepPoint()->Ge
165 }
166 else
167 {
168 fStatus = NotAtBoundary;
169 if(verboseLevel > 1)
170 BoundaryProcessVerbose();
171 return G4VDiscreteProcess::PostStepDoIt(aT
172 }
173
174 G4VPhysicalVolume* thePrePV = pStep->GetPre
175 G4VPhysicalVolume* thePostPV = pStep->GetPos
176
177 if(verboseLevel > 1)
178 {
179 G4cout << " Photon at Boundary! " << G4end
180 if(thePrePV != nullptr)
181 G4cout << " thePrePV: " << thePrePV->Ge
182 if(thePostPV != nullptr)
183 G4cout << " thePostPV: " << thePostPV->G
184 }
185
186 G4double stepLength = aTrack.GetStepLength()
187 if(stepLength <= fCarTolerance)
188 {
189 fStatus = StepTooSmall;
190 if(verboseLevel > 1)
191 BoundaryProcessVerbose();
192
193 G4MaterialPropertyVector* groupvel = nullp
194 G4MaterialPropertiesTable* aMPT = fMateria
195 if(aMPT != nullptr)
196 {
197 groupvel = aMPT->GetProperty(kGROUPVEL);
198 }
199
200 if(groupvel != nullptr)
201 {
202 aParticleChange.ProposeVelocity(
203 groupvel->Value(fPhotonMomentum, idx_g
204 }
205 return G4VDiscreteProcess::PostStepDoIt(aT
206 }
207 else if (stepLength <= 10.*fCarTolerance &&
208 { // see bug 2510
209 ++fNumSmallStepWarnings;
210 if(verboseLevel > 0)
211 {
212 G4ExceptionDescription ed;
213 ed << "G4OpBoundaryProcess: "
214 << "Opticalphoton step length: " << s
215 << "This is larger than the threshold
216 "to set status StepTooSmall." << G
217 << "Boundary scattering may be incorr
218 if(fNumSmallStepWarnings == 10)
219 {
220 ed << G4endl << "*** Step size warning
221 }
222 G4Exception("G4OpBoundaryProcess", "OpBo
223 }
224 }
225
226 const G4DynamicParticle* aParticle = aTrack.
227
228 fPhotonMomentum = aParticle->GetTotalMoment
229 fOldMomentum = aParticle->GetMomentumDir
230 fOldPolarization = aParticle->GetPolarizatio
231
232 if(verboseLevel > 1)
233 {
234 G4cout << " Old Momentum Direction: " << f
235 << " Old Polarization: " << f
236 }
237
238 G4ThreeVector theGlobalPoint = pStep->GetPos
239 G4bool valid;
240
241 // ID of Navigator which limits step
242 G4int hNavId = G4ParallelWorldProcess::GetHy
243 auto iNav = G4TransportationManager::GetT
244 ->GetActiveNavigatorsIterator(
245 fGlobalNormal = (iNav[hNavId])->GetGlobalExi
246
247 if(valid)
248 {
249 fGlobalNormal = -fGlobalNormal;
250 }
251 else
252 {
253 G4ExceptionDescription ed;
254 ed << " G4OpBoundaryProcess/PostStepDoIt()
255 << " The Navigator reports that it retu
256 G4Exception(
257 "G4OpBoundaryProcess::PostStepDoIt", "Op
258 "Invalid Surface Normal - Geometry must
259 }
260
261 if(fOldMomentum * fGlobalNormal > 0.0)
262 {
263 #ifdef G4OPTICAL_DEBUG
264 G4ExceptionDescription ed;
265 ed << " G4OpBoundaryProcess/PostStepDoIt()
266 "wrong direction. "
267 << G4endl
268 << " The momentum of the photon arriv
269 << " must exit the volume cross in th
270 << " So it MUST have dot < 0 with the
271 "volume (globalNormal)."
272 << G4endl << " >> The dot product of
273 << fOldMomentum * fGlobalNormal << G4en
274 << " Old Momentum (during step)
275 << " Global Normal (Exiting New Vol
276 << G4endl;
277 G4Exception("G4OpBoundaryProcess::PostStep
278 EventMustBeAborted, // Or Jus
279 // repeat
280 ed,
281 "Invalid Surface Normal - Geom
282 "normal pointing in the right
283 #else
284 fGlobalNormal = -fGlobalNormal;
285 #endif
286 }
287
288 G4MaterialPropertyVector* rIndexMPV = nullpt
289 G4MaterialPropertiesTable* MPT = fMaterial1-
290 if(MPT != nullptr)
291 {
292 rIndexMPV = MPT->GetProperty(kRINDEX);
293 }
294 if(rIndexMPV != nullptr)
295 {
296 fRindex1 = rIndexMPV->Value(fPhotonMomentu
297 }
298 else
299 {
300 fStatus = NoRINDEX;
301 if(verboseLevel > 1)
302 BoundaryProcessVerbose();
303 aParticleChange.ProposeLocalEnergyDeposit(
304 aParticleChange.ProposeTrackStatus(fStopAn
305 return G4VDiscreteProcess::PostStepDoIt(aT
306 }
307
308 fReflectivity = 1.;
309 fEfficiency = 0.;
310 fTransmittance = 0.;
311 fSurfaceRoughness = 0.;
312 fModel = glisur;
313 fFinish = polished;
314 G4SurfaceType type = dielectric_dielectric;
315
316 rIndexMPV = nullptr;
317 fOpticalSurface = nullptr;
318
319 G4LogicalSurface* surface =
320 G4LogicalBorderSurface::GetSurface(thePreP
321 if(surface == nullptr)
322 {
323 if(thePostPV->GetMotherLogical() == thePre
324 {
325 surface = G4LogicalSkinSurface::GetSurfa
326 if(surface == nullptr)
327 {
328 surface =
329 G4LogicalSkinSurface::GetSurface(the
330 }
331 }
332 else
333 {
334 surface = G4LogicalSkinSurface::GetSurfa
335 if(surface == nullptr)
336 {
337 surface =
338 G4LogicalSkinSurface::GetSurface(the
339 }
340 }
341 }
342
343 if(surface != nullptr)
344 {
345 fOpticalSurface =
346 dynamic_cast<G4OpticalSurface*>(surface-
347 }
348 if(fOpticalSurface != nullptr)
349 {
350 type = fOpticalSurface->GetType();
351 fModel = fOpticalSurface->GetModel();
352 fFinish = fOpticalSurface->GetFinish();
353
354 G4MaterialPropertiesTable* sMPT =
355 fOpticalSurface->GetMaterialPropertiesTa
356 if(sMPT != nullptr)
357 {
358 if(fFinish == polishedbackpainted || fFi
359 {
360 rIndexMPV = sMPT->GetProperty(kRINDEX)
361 if(rIndexMPV != nullptr)
362 {
363 fRindex2 = rIndexMPV->Value(fPhotonM
364 }
365 else
366 {
367 fStatus = NoRINDEX;
368 if(verboseLevel > 1)
369 BoundaryProcessVerbose();
370 aParticleChange.ProposeLocalEnergyDe
371 aParticleChange.ProposeTrackStatus(f
372 return G4VDiscreteProcess::PostStepD
373 }
374 }
375
376 fRealRIndexMPV = sMPT->GetProperty(kREAL
377 fImagRIndexMPV = sMPT->GetProperty(kIMAG
378 f_iTE = f_iTM = 1;
379
380 G4MaterialPropertyVector* pp;
381 if((pp = sMPT->GetProperty(kREFLECTIVITY
382 {
383 fReflectivity = pp->Value(fPhotonMomen
384 }
385 else if(fRealRIndexMPV && fImagRIndexMPV
386 {
387 CalculateReflectivity();
388 }
389
390 if((pp = sMPT->GetProperty(kEFFICIENCY))
391 {
392 fEfficiency = pp->Value(fPhotonMomentu
393 }
394 if((pp = sMPT->GetProperty(kTRANSMITTANC
395 {
396 fTransmittance = pp->Value(fPhotonMome
397 }
398 if(sMPT->ConstPropertyExists(kSURFACEROU
399 {
400 fSurfaceRoughness = sMPT->GetConstProp
401 }
402
403 if(fModel == unified)
404 {
405 fProb_sl = (pp = sMPT->GetProperty(kSP
406 ? pp->Value(fPhotonMoment
407 : 0.;
408 fProb_ss = (pp = sMPT->GetProperty(kSP
409 ? pp->Value(fPhotonMoment
410 : 0.;
411 fProb_bs = (pp = sMPT->GetProperty(kBA
412 ? pp->Value(fPhotonMoment
413 : 0.;
414 }
415 } // end of if(sMPT)
416 else if(fFinish == polishedbackpainted ||
417 {
418 aParticleChange.ProposeLocalEnergyDeposi
419 aParticleChange.ProposeTrackStatus(fStop
420 return G4VDiscreteProcess::PostStepDoIt(
421 }
422 } // end of if(fOpticalSurface)
423
424 // DIELECTRIC-DIELECTRIC
425 if(type == dielectric_dielectric)
426 {
427 if(fFinish == polished || fFinish == groun
428 {
429 if(fMaterial1 == fMaterial2)
430 {
431 fStatus = SameMaterial;
432 if(verboseLevel > 1)
433 BoundaryProcessVerbose();
434 return G4VDiscreteProcess::PostStepDoI
435 }
436 MPT = fMaterial2->GetMaterialPrope
437 rIndexMPV = nullptr;
438 if(MPT != nullptr)
439 {
440 rIndexMPV = MPT->GetProperty(kRINDEX);
441 }
442 if(rIndexMPV != nullptr)
443 {
444 fRindex2 = rIndexMPV->Value(fPhotonMom
445 }
446 else
447 {
448 fStatus = NoRINDEX;
449 if(verboseLevel > 1)
450 BoundaryProcessVerbose();
451 aParticleChange.ProposeLocalEnergyDepo
452 aParticleChange.ProposeTrackStatus(fSt
453 return G4VDiscreteProcess::PostStepDoI
454 }
455 }
456 if(fFinish == polishedbackpainted || fFini
457 {
458 DielectricDielectric();
459 }
460 else
461 {
462 G4double rand = G4UniformRand();
463 if(rand > fReflectivity + fTransmittance
464 {
465 DoAbsorption();
466 }
467 else if(rand > fReflectivity)
468 {
469 fStatus = Transmission;
470 fNewMomentum = fOldMomentum;
471 fNewPolarization = fOldPolarization;
472 }
473 else
474 {
475 if(fFinish == polishedfrontpainted)
476 {
477 DoReflection();
478 }
479 else if(fFinish == groundfrontpainted)
480 {
481 fStatus = LambertianReflection;
482 DoReflection();
483 }
484 else
485 {
486 DielectricDielectric();
487 }
488 }
489 }
490 }
491 else if(type == dielectric_metal)
492 {
493 DielectricMetal();
494 }
495 else if(type == dielectric_LUT)
496 {
497 DielectricLUT();
498 }
499 else if(type == dielectric_LUTDAVIS)
500 {
501 DielectricLUTDAVIS();
502 }
503 else if(type == dielectric_dichroic)
504 {
505 DielectricDichroic();
506 }
507 else if(type == coated)
508 {
509 CoatedDielectricDielectric();
510 }
511 else
512 {
513 if(fNumBdryTypeWarnings <= 10)
514 {
515 ++fNumBdryTypeWarnings;
516 if(verboseLevel > 0)
517 {
518 G4ExceptionDescription ed;
519 ed << " PostStepDoIt(): Illegal bounda
520 if(fNumBdryTypeWarnings == 10)
521 {
522 ed << "** Boundary type warnings sto
523 }
524 G4Exception("G4OpBoundaryProcess", "Op
525 }
526 }
527 return G4VDiscreteProcess::PostStepDoIt(aT
528 }
529
530 fNewMomentum = fNewMomentum.unit();
531 fNewPolarization = fNewPolarization.unit();
532
533 if(verboseLevel > 1)
534 {
535 G4cout << " New Momentum Direction: " << f
536 << " New Polarization: " << f
537 BoundaryProcessVerbose();
538 }
539
540 aParticleChange.ProposeMomentumDirection(fNe
541 aParticleChange.ProposePolarization(fNewPola
542
543 if(fStatus == FresnelRefraction || fStatus =
544 {
545 // not all surface types check that fMater
546 G4MaterialPropertiesTable* aMPT = fMateria
547 G4MaterialPropertyVector* groupvel = nullp
548 if(aMPT != nullptr)
549 {
550 groupvel = aMPT->GetProperty(kGROUPVEL);
551 }
552 if(groupvel != nullptr)
553 {
554 aParticleChange.ProposeVelocity(
555 groupvel->Value(fPhotonMomentum, idx_g
556 }
557 }
558
559 if(fStatus == Detection && fInvokeSD)
560 InvokeSD(pStep);
561 return G4VDiscreteProcess::PostStepDoIt(aTra
562 }
563
564 //....oooOO0OOooo........oooOO0OOooo........oo
565 void G4OpBoundaryProcess::BoundaryProcessVerbo
566 {
567 G4cout << " *** ";
568 if(fStatus == Undefined)
569 G4cout << "Undefined";
570 else if(fStatus == Transmission)
571 G4cout << "Transmission";
572 else if(fStatus == FresnelRefraction)
573 G4cout << "FresnelRefraction";
574 else if(fStatus == FresnelReflection)
575 G4cout << "FresnelReflection";
576 else if(fStatus == TotalInternalReflection)
577 G4cout << "TotalInternalReflection";
578 else if(fStatus == LambertianReflection)
579 G4cout << "LambertianReflection";
580 else if(fStatus == LobeReflection)
581 G4cout << "LobeReflection";
582 else if(fStatus == SpikeReflection)
583 G4cout << "SpikeReflection";
584 else if(fStatus == BackScattering)
585 G4cout << "BackScattering";
586 else if(fStatus == PolishedLumirrorAirReflec
587 G4cout << "PolishedLumirrorAirReflection";
588 else if(fStatus == PolishedLumirrorGlueRefle
589 G4cout << "PolishedLumirrorGlueReflection"
590 else if(fStatus == PolishedAirReflection)
591 G4cout << "PolishedAirReflection";
592 else if(fStatus == PolishedTeflonAirReflecti
593 G4cout << "PolishedTeflonAirReflection";
594 else if(fStatus == PolishedTiOAirReflection)
595 G4cout << "PolishedTiOAirReflection";
596 else if(fStatus == PolishedTyvekAirReflectio
597 G4cout << "PolishedTyvekAirReflection";
598 else if(fStatus == PolishedVM2000AirReflecti
599 G4cout << "PolishedVM2000AirReflection";
600 else if(fStatus == PolishedVM2000GlueReflect
601 G4cout << "PolishedVM2000GlueReflection";
602 else if(fStatus == EtchedLumirrorAirReflecti
603 G4cout << "EtchedLumirrorAirReflection";
604 else if(fStatus == EtchedLumirrorGlueReflect
605 G4cout << "EtchedLumirrorGlueReflection";
606 else if(fStatus == EtchedAirReflection)
607 G4cout << "EtchedAirReflection";
608 else if(fStatus == EtchedTeflonAirReflection
609 G4cout << "EtchedTeflonAirReflection";
610 else if(fStatus == EtchedTiOAirReflection)
611 G4cout << "EtchedTiOAirReflection";
612 else if(fStatus == EtchedTyvekAirReflection)
613 G4cout << "EtchedTyvekAirReflection";
614 else if(fStatus == EtchedVM2000AirReflection
615 G4cout << "EtchedVM2000AirReflection";
616 else if(fStatus == EtchedVM2000GlueReflectio
617 G4cout << "EtchedVM2000GlueReflection";
618 else if(fStatus == GroundLumirrorAirReflecti
619 G4cout << "GroundLumirrorAirReflection";
620 else if(fStatus == GroundLumirrorGlueReflect
621 G4cout << "GroundLumirrorGlueReflection";
622 else if(fStatus == GroundAirReflection)
623 G4cout << "GroundAirReflection";
624 else if(fStatus == GroundTeflonAirReflection
625 G4cout << "GroundTeflonAirReflection";
626 else if(fStatus == GroundTiOAirReflection)
627 G4cout << "GroundTiOAirReflection";
628 else if(fStatus == GroundTyvekAirReflection)
629 G4cout << "GroundTyvekAirReflection";
630 else if(fStatus == GroundVM2000AirReflection
631 G4cout << "GroundVM2000AirReflection";
632 else if(fStatus == GroundVM2000GlueReflectio
633 G4cout << "GroundVM2000GlueReflection";
634 else if(fStatus == Absorption)
635 G4cout << "Absorption";
636 else if(fStatus == Detection)
637 G4cout << "Detection";
638 else if(fStatus == NotAtBoundary)
639 G4cout << "NotAtBoundary";
640 else if(fStatus == SameMaterial)
641 G4cout << "SameMaterial";
642 else if(fStatus == StepTooSmall)
643 G4cout << "StepTooSmall";
644 else if(fStatus == NoRINDEX)
645 G4cout << "NoRINDEX";
646 else if(fStatus == Dichroic)
647 G4cout << "Dichroic Transmission";
648 else if(fStatus == CoatedDielectricReflectio
649 G4cout << "Coated Dielectric Reflection";
650 else if(fStatus == CoatedDielectricRefractio
651 G4cout << "Coated Dielectric Refraction";
652 else if(fStatus == CoatedDielectricFrustrate
653 G4cout << "Coated Dielectric Frustrated Tr
654
655 G4cout << " ***" << G4endl;
656 }
657
658 //....oooOO0OOooo........oooOO0OOooo........oo
659 G4ThreeVector G4OpBoundaryProcess::GetFacetNor
660 const G4ThreeVector& momentum, const G4Three
661 {
662 G4ThreeVector facetNormal;
663 if(fModel == unified || fModel == LUT || fMo
664 {
665 /* This function codes alpha to a random v
666 distribution p(alpha) = g(alpha; 0, sigma_
667 for alpha > 0 and alpha < 90, where g(alph
668 gaussian distribution with mean 0 and stan
669
670 G4double sigma_alpha = 0.0;
671 if(fOpticalSurface)
672 sigma_alpha = fOpticalSurface->GetSigmaA
673 if(sigma_alpha == 0.0)
674 {
675 return normal;
676 }
677
678 G4double f_max = std::min(1.0, 4. * sigma_
679 G4double alpha, phi, sinAlpha;
680
681 do
682 { // Loop checking, 13-Aug-2015, Peter Gu
683 do
684 { // Loop checking, 13-Aug-2015, Peter
685 alpha = G4RandGauss::shoot(0.0, sig
686 sinAlpha = std::sin(alpha);
687 } while(G4UniformRand() * f_max > sinAlp
688
689 phi = G4UniformRand() * twopi;
690 facetNormal.set(sinAlpha * std::cos(phi)
691 std::cos(alpha));
692 facetNormal.rotateUz(normal);
693 } while(momentum * facetNormal >= 0.0);
694 }
695 else
696 {
697 G4double polish = 1.0;
698 if(fOpticalSurface)
699 polish = fOpticalSurface->GetPolish();
700 if(polish < 1.0)
701 {
702 do
703 { // Loop checking, 13-Aug-2015, Peter
704 G4ThreeVector smear;
705 do
706 { // Loop checking, 13-Aug-2015, Pete
707 smear.setX(2. * G4UniformRand() - 1.
708 smear.setY(2. * G4UniformRand() - 1.
709 smear.setZ(2. * G4UniformRand() - 1.
710 } while(smear.mag2() > 1.0);
711 facetNormal = normal + (1. - polish) *
712 } while(momentum * facetNormal >= 0.0);
713 facetNormal = facetNormal.unit();
714 }
715 else
716 {
717 facetNormal = normal;
718 }
719 }
720 return facetNormal;
721 }
722
723 //....oooOO0OOooo........oooOO0OOooo........oo
724 void G4OpBoundaryProcess::DielectricMetal()
725 {
726 G4int n = 0;
727 G4double rand;
728 G4ThreeVector A_trans;
729
730 do
731 {
732 ++n;
733 rand = G4UniformRand();
734 if(rand > fReflectivity && n == 1)
735 {
736 if(rand > fReflectivity + fTransmittance
737 {
738 DoAbsorption();
739 }
740 else
741 {
742 fStatus = Transmission;
743 fNewMomentum = fOldMomentum;
744 fNewPolarization = fOldPolarization;
745 }
746 break;
747 }
748 else
749 {
750 if(fRealRIndexMPV && fImagRIndexMPV)
751 {
752 if(n > 1)
753 {
754 CalculateReflectivity();
755 if(!G4BooleanRand(fReflectivity))
756 {
757 DoAbsorption();
758 break;
759 }
760 }
761 }
762 if(fModel == glisur || fFinish == polish
763 {
764 DoReflection();
765 }
766 else
767 {
768 if(n == 1)
769 ChooseReflection();
770 if(fStatus == LambertianReflection)
771 {
772 DoReflection();
773 }
774 else if(fStatus == BackScattering)
775 {
776 fNewMomentum = -fOldMomentum;
777 fNewPolarization = -fOldPolarization
778 }
779 else
780 {
781 if(fStatus == LobeReflection)
782 {
783 if(!fRealRIndexMPV || !fImagRIndex
784 {
785 fFacetNormal = GetFacetNormal(fO
786 }
787 // else
788 // case of complex rindex needs t
789 }
790 fNewMomentum =
791 fOldMomentum - 2. * fOldMomentum *
792
793 if(f_iTE > 0 && f_iTM > 0)
794 {
795 fNewPolarization =
796 -fOldPolarization +
797 (2. * fOldPolarization * fFacetN
798 }
799 else if(f_iTE > 0)
800 {
801 A_trans = (fSint1 > 0.0) ? fOldMom
802 : fOldPol
803 fNewPolarization = -A_trans;
804 }
805 else if(f_iTM > 0)
806 {
807 fNewPolarization =
808 -fNewMomentum.cross(A_trans).uni
809 }
810 }
811 }
812 fOldMomentum = fNewMomentum;
813 fOldPolarization = fNewPolarization;
814 }
815 // Loop checking, 13-Aug-2015, Peter Gumpl
816 } while(fNewMomentum * fGlobalNormal < 0.0);
817 }
818
819 //....oooOO0OOooo........oooOO0OOooo........oo
820 void G4OpBoundaryProcess::DielectricLUT()
821 {
822 G4int thetaIndex, phiIndex;
823 G4double angularDistVal, thetaRad, phiRad;
824 G4ThreeVector perpVectorTheta, perpVectorPhi
825
826 fStatus = G4OpBoundaryProcessStatus(
827 G4int(fFinish) + (G4int(NoRINDEX) - G4int(
828
829 G4int thetaIndexMax = fOpticalSurface->GetTh
830 G4int phiIndexMax = fOpticalSurface->GetPh
831
832 G4double rand;
833
834 do
835 {
836 rand = G4UniformRand();
837 if(rand > fReflectivity)
838 {
839 if(rand > fReflectivity + fTransmittance
840 {
841 DoAbsorption();
842 }
843 else
844 {
845 fStatus = Transmission;
846 fNewMomentum = fOldMomentum;
847 fNewPolarization = fOldPolarization;
848 }
849 break;
850 }
851 else
852 {
853 // Calculate Angle between Normal and Ph
854 G4double anglePhotonToNormal = fOldMomen
855 // Round to closest integer: LBNL model
856 G4int angleIncident = (G4int)std::lrint(
857
858 // Take random angles THETA and PHI,
859 // and see if below Probability - if not
860 do
861 {
862 thetaIndex = (G4int)G4RandFlat::shootI
863 phiIndex = (G4int)G4RandFlat::shootI
864 // Find probability with the new indec
865 angularDistVal = fOpticalSurface->GetA
866 angleIncident, thetaIndex, phiIndex)
867 // Loop checking, 13-Aug-2015, Peter G
868 } while(!G4BooleanRand(angularDistVal));
869
870 thetaRad = G4double(-90 + 4 * thetaIndex
871 phiRad = G4double(-90 + 5 * phiIndex)
872 // Rotate Photon Momentum in Theta, then
873 fNewMomentum = -fOldMomentum;
874
875 perpVectorTheta = fNewMomentum.cross(fGl
876 if(perpVectorTheta.mag() < fCarTolerance
877 {
878 perpVectorTheta = fNewMomentum.orthogo
879 }
880 fNewMomentum =
881 fNewMomentum.rotate(anglePhotonToNorma
882 perpVectorPhi = perpVectorTheta.cross(fN
883 fNewMomentum = fNewMomentum.rotate(-phi
884
885 // Rotate Polarization too:
886 fFacetNormal = (fNewMomentum - fOldM
887 fNewPolarization = -fOldPolarization +
888 (2. * fOldPolarizatio
889 }
890 // Loop checking, 13-Aug-2015, Peter Gumpl
891 } while(fNewMomentum * fGlobalNormal <= 0.0)
892 }
893
894 //....oooOO0OOooo........oooOO0OOooo........oo
895 void G4OpBoundaryProcess::DielectricLUTDAVIS()
896 {
897 G4int angindex, random, angleIncident;
898 G4double reflectivityValue, elevation, azimu
899 G4double anglePhotonToNormal;
900
901 G4int lutbin = fOpticalSurface->GetLUTbins(
902 G4double rand = G4UniformRand();
903
904 G4double sinEl;
905 G4ThreeVector u, vNorm, w;
906
907 do
908 {
909 anglePhotonToNormal = fOldMomentum.angle(-
910
911 // Davis model has 90 reflection bins: rou
912 // don't allow angleIncident to be 90 for
913 angleIncident = std::min(
914 static_cast<G4int>(std::floor(anglePhoto
915 reflectivityValue = fOpticalSurface->GetRe
916
917 if(rand > reflectivityValue)
918 {
919 if(fEfficiency > 0.)
920 {
921 DoAbsorption();
922 break;
923 }
924 else
925 {
926 fStatus = Transmission;
927
928 if(angleIncident <= 0.01)
929 {
930 fNewMomentum = fOldMomentum;
931 break;
932 }
933
934 do
935 {
936 random = (G4int)G4RandFlat::shootInt
937 angindex =
938 (((random * 2) - 1)) + angleIncide
939
940 azimuth =
941 fOpticalSurface->GetAngularDistrib
942 elevation = fOpticalSurface->GetAngu
943 } while(elevation == 0. && azimuth ==
944
945 sinEl = std::sin(elevation);
946 vNorm = (fGlobalNormal.cross(fOldMomen
947 u = vNorm.cross(fGlobalNormal) * (
948 vNorm *= (sinEl * std::sin(azimuth));
949 // fGlobalNormal shouldn't be modified
950 w = (fGlobalNormal *= std::
951 fNewMomentum = u + vNorm + w;
952
953 // Rotate Polarization too:
954 fFacetNormal = (fNewMomentum - fOl
955 fNewPolarization = -fOldPolarization +
956
957 }
958 }
959 else
960 {
961 fStatus = LobeReflection;
962
963 if(angleIncident == 0)
964 {
965 fNewMomentum = -fOldMomentum;
966 break;
967 }
968
969 do
970 {
971 random = (G4int)G4RandFlat::shootInt
972 angindex = (((random * 2) - 1)) + (ang
973
974 azimuth = fOpticalSurface->GetAngularD
975 elevation = fOpticalSurface->GetAngula
976 } while(elevation == 0. && azimuth == 0.
977
978 sinEl = std::sin(elevation);
979 vNorm = (fGlobalNormal.cross(fOldMomentu
980 u = vNorm.cross(fGlobalNormal) * (si
981 vNorm *= (sinEl * std::sin(azimuth));
982 // fGlobalNormal shouldn't be modified h
983 w = (fGlobalNormal *= std::cos(elevation
984
985 fNewMomentum = u + vNorm + w;
986
987 // Rotate Polarization too: (needs revis
988 fNewPolarization = fOldPolarization;
989 }
990 } while(fNewMomentum * fGlobalNormal <= 0.0)
991 }
992
993 //....oooOO0OOooo........oooOO0OOooo........oo
994 void G4OpBoundaryProcess::DielectricDichroic()
995 {
996 // Calculate Angle between Normal and Photon
997 G4double anglePhotonToNormal = fOldMomentum.
998
999 // Round it to closest integer
1000 G4double angleIncident = std::floor(180. /
1001
1002 if(!fDichroicVector)
1003 {
1004 if(fOpticalSurface)
1005 fDichroicVector = fOpticalSurface->GetD
1006 }
1007
1008 if(fDichroicVector)
1009 {
1010 G4double wavelength = h_Planck * c_light
1011 fTransmittance = fDichroicVector->Va
1012 i
1013 perCent;
1014 // G4cout << "wavelength: " << std::flo
1015 // << "nm" <<
1016 // G4cout << "Incident angle: " << angl
1017 // G4cout << "Transmittance: "
1018 // << std::floor(fTransmittance/
1019 }
1020 else
1021 {
1022 G4ExceptionDescription ed;
1023 ed << " G4OpBoundaryProcess/DielectricDic
1024 << " The dichroic surface has no G4Phy
1025 G4Exception("G4OpBoundaryProcess::Dielect
1026 FatalException, ed,
1027 "A dichroic surface must have
1028 }
1029
1030 if(!G4BooleanRand(fTransmittance))
1031 { // Not transmitted, so reflect
1032 if(fModel == glisur || fFinish == polishe
1033 {
1034 DoReflection();
1035 }
1036 else
1037 {
1038 ChooseReflection();
1039 if(fStatus == LambertianReflection)
1040 {
1041 DoReflection();
1042 }
1043 else if(fStatus == BackScattering)
1044 {
1045 fNewMomentum = -fOldMomentum;
1046 fNewPolarization = -fOldPolarization;
1047 }
1048 else
1049 {
1050 G4double PdotN, EdotN;
1051 do
1052 {
1053 if(fStatus == LobeReflection)
1054 {
1055 fFacetNormal = GetFacetNormal(fOl
1056 }
1057 PdotN = fOldMomentum * fFace
1058 fNewMomentum = fOldMomentum - (2. *
1059 // Loop checking, 13-Aug-2015, Pete
1060 } while(fNewMomentum * fGlobalNormal
1061
1062 EdotN = fOldPolarization *
1063 fNewPolarization = -fOldPolarization
1064 }
1065 }
1066 }
1067 else
1068 {
1069 fStatus = Dichroic;
1070 fNewMomentum = fOldMomentum;
1071 fNewPolarization = fOldPolarization;
1072 }
1073 }
1074
1075 //....oooOO0OOooo........oooOO0OOooo........o
1076 void G4OpBoundaryProcess::DielectricDielectri
1077 {
1078 G4bool inside = false;
1079 G4bool swap = false;
1080
1081 if(fFinish == polished)
1082 {
1083 fFacetNormal = fGlobalNormal;
1084 }
1085 else
1086 {
1087 fFacetNormal = GetFacetNormal(fOldMomentu
1088 }
1089 G4double cost1 = -fOldMomentum * fFacetNorm
1090 G4double cost2 = 0.;
1091 G4double sint2 = 0.;
1092
1093 G4bool surfaceRoughnessCriterionPass = true
1094 if(fSurfaceRoughness != 0. && fRindex1 > fR
1095 {
1096 G4double wavelength = h_Pl
1097 G4double surfaceRoughnessCriterion = std:
1098 (4. * pi * fSurfaceRoughness * fRindex1
1099 surfaceRoughnessCriterionPass = G4Boolean
1100 }
1101
1102 leap:
1103
1104 G4bool through = false;
1105 G4bool done = false;
1106
1107 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1108 G4double E1_perp, E1_parl;
1109 G4double s1, s2, E2_perp, E2_parl, E2_total
1110 G4double E2_abs, C_parl, C_perp;
1111 G4double alpha;
1112
1113 do
1114 {
1115 if(through)
1116 {
1117 swap = !swap;
1118 through = false;
1119 fGlobalNormal = -fGlobalNormal;
1120 G4SwapPtr(fMaterial1, fMaterial2);
1121 G4SwapObj(&fRindex1, &fRindex2);
1122 }
1123
1124 if(fFinish == polished)
1125 {
1126 fFacetNormal = fGlobalNormal;
1127 }
1128 else
1129 {
1130 fFacetNormal = GetFacetNormal(fOldMomen
1131 }
1132
1133 cost1 = -fOldMomentum * fFacetNormal;
1134 if(std::abs(cost1) < 1.0 - fCarTolerance)
1135 {
1136 fSint1 = std::sqrt(1. - cost1 * cost1);
1137 sint2 = fSint1 * fRindex1 / fRindex2;
1138 // this isn't a sine as we might expect
1139 }
1140 else
1141 {
1142 fSint1 = 0.0;
1143 sint2 = 0.0;
1144 }
1145
1146 // TOTAL INTERNAL REFLECTION
1147 if(sint2 >= 1.0)
1148 {
1149 swap = false;
1150
1151 fStatus = TotalInternalReflection;
1152 if(!surfaceRoughnessCriterionPass)
1153 fStatus = LambertianReflection;
1154 if(fModel == unified && fFinish != poli
1155 ChooseReflection();
1156 if(fStatus == LambertianReflection)
1157 {
1158 DoReflection();
1159 }
1160 else if(fStatus == BackScattering)
1161 {
1162 fNewMomentum = -fOldMomentum;
1163 fNewPolarization = -fOldPolarization;
1164 }
1165 else
1166 {
1167 fNewMomentum =
1168 fOldMomentum - 2. * fOldMomentum *
1169 fNewPolarization = -fOldPolarization
1170
1171 }
1172 }
1173 // NOT TIR
1174 else if(sint2 < 1.0)
1175 {
1176 // Calculate amplitude for transmission
1177 if(cost1 > 0.0)
1178 {
1179 cost2 = std::sqrt(1. - sint2 * sint2)
1180 }
1181 else
1182 {
1183 cost2 = -std::sqrt(1. - sint2 * sint2
1184 }
1185
1186 if(fSint1 > 0.0)
1187 {
1188 A_trans = (fOldMomentum.cross(fFacetN
1189 E1_perp = fOldPolarization * A_trans;
1190 E1pp = E1_perp * A_trans;
1191 E1pl = fOldPolarization - E1pp;
1192 E1_parl = E1pl.mag();
1193 }
1194 else
1195 {
1196 A_trans = fOldPolarization;
1197 // Here we Follow Jackson's conventio
1198 // component = 1 in case of a ray per
1199 E1_perp = 0.0;
1200 E1_parl = 1.0;
1201 }
1202
1203 s1 = fRindex1 * cost1;
1204 E2_perp = 2. * s1 * E1_perp / (fRindex
1205 E2_parl = 2. * s1 * E1_parl / (fRindex
1206 E2_total = E2_perp * E2_perp + E2_parl
1207 s2 = fRindex2 * cost2 * E2_total;
1208
1209 // D.Sawkey, 24 May 24
1210 // Transmittance has already been taken
1211 // For e.g. specular surfaces, the rati
1212 // reflection should be given by the ma
1213 // TRANSMITTANCE
1214 //if(fTransmittance > 0.)
1215 // transCoeff = fTransmittance;
1216 //else if(cost1 != 0.0)
1217 if(cost1 != 0.0)
1218 transCoeff = s2 / s1;
1219 else
1220 transCoeff = 0.0;
1221
1222 // NOT TIR: REFLECTION
1223 if(!G4BooleanRand(transCoeff))
1224 {
1225 swap = false;
1226 fStatus = FresnelReflection;
1227
1228 if(!surfaceRoughnessCriterionPass)
1229 fStatus = LambertianReflection;
1230 if(fModel == unified && fFinish != po
1231 ChooseReflection();
1232 if(fStatus == LambertianReflection)
1233 {
1234 DoReflection();
1235 }
1236 else if(fStatus == BackScattering)
1237 {
1238 fNewMomentum = -fOldMomentum;
1239 fNewPolarization = -fOldPolarizatio
1240 }
1241 else
1242 {
1243 fNewMomentum =
1244 fOldMomentum - 2. * fOldMomentum
1245 if(fSint1 > 0.0)
1246 { // incident ray oblique
1247 E2_parl = fRindex2 * E2_parl / f
1248 E2_perp = E2_perp - E1_perp;
1249 E2_total = E2_perp * E2_perp + E2
1250 A_paral = (fNewMomentum.cross(A_
1251 E2_abs = std::sqrt(E2_total);
1252 C_parl = E2_parl / E2_abs;
1253 C_perp = E2_perp / E2_abs;
1254
1255 fNewPolarization = C_parl * A_par
1256 }
1257 else
1258 { // incident ray perpendicular
1259 if(fRindex2 > fRindex1)
1260 {
1261 fNewPolarization = -fOldPolariz
1262 }
1263 else
1264 {
1265 fNewPolarization = fOldPolariza
1266 }
1267 }
1268 }
1269 }
1270 // NOT TIR: TRANSMISSION
1271 else
1272 {
1273 inside = !inside;
1274 through = true;
1275 fStatus = FresnelRefraction;
1276
1277 if(fSint1 > 0.0)
1278 { // incident ray oblique
1279 alpha = cost1 - cost2 * (fRi
1280 fNewMomentum = (fOldMomentum + alph
1281 A_paral = (fNewMomentum.cross(
1282 E2_abs = std::sqrt(E2_total);
1283 C_parl = E2_parl / E2_abs;
1284 C_perp = E2_perp / E2_abs;
1285
1286 fNewPolarization = C_parl * A_paral
1287 }
1288 else
1289 { // incident ray perpendicular
1290 fNewMomentum = fOldMomentum;
1291 fNewPolarization = fOldPolarization
1292 }
1293 }
1294 }
1295
1296 fOldMomentum = fNewMomentum.unit();
1297 fOldPolarization = fNewPolarization.unit(
1298
1299 if(fStatus == FresnelRefraction)
1300 {
1301 done = (fNewMomentum * fGlobalNormal <=
1302 }
1303 else
1304 {
1305 done = (fNewMomentum * fGlobalNormal >=
1306 }
1307 // Loop checking, 13-Aug-2015, Peter Gump
1308 } while(!done);
1309
1310 if(inside && !swap)
1311 {
1312 if(fFinish == polishedbackpainted || fFin
1313 {
1314 G4double rand = G4UniformRand();
1315 if(rand > fReflectivity + fTransmittanc
1316 {
1317 DoAbsorption();
1318 }
1319 else if(rand > fReflectivity)
1320 {
1321 fStatus = Transmission;
1322 fNewMomentum = fOldMomentum;
1323 fNewPolarization = fOldPolarization;
1324 }
1325 else
1326 {
1327 if(fStatus != FresnelRefraction)
1328 {
1329 fGlobalNormal = -fGlobalNormal;
1330 }
1331 else
1332 {
1333 swap = !swap;
1334 G4SwapPtr(fMaterial1, fMaterial2);
1335 G4SwapObj(&fRindex1, &fRindex2);
1336 }
1337 if(fFinish == groundbackpainted)
1338 fStatus = LambertianReflection;
1339
1340 DoReflection();
1341
1342 fGlobalNormal = -fGlobalNormal;
1343 fOldMomentum = fNewMomentum;
1344
1345 goto leap;
1346 }
1347 }
1348 }
1349 }
1350
1351 //....oooOO0OOooo........oooOO0OOooo........o
1352 G4double G4OpBoundaryProcess::GetMeanFreePath
1353
1354 {
1355 *condition = Forced;
1356 return DBL_MAX;
1357 }
1358
1359 //....oooOO0OOooo........oooOO0OOooo........o
1360 G4double G4OpBoundaryProcess::GetIncidentAngl
1361 {
1362 return pi - std::acos(fOldMomentum * fFacet
1363 (fOldMomentum.mag() *
1364 }
1365
1366 //....oooOO0OOooo........oooOO0OOooo........o
1367 G4double G4OpBoundaryProcess::GetReflectivity
1368
1369
1370
1371
1372 {
1373 G4complex reflectivity, reflectivity_TE, re
1374 G4complex N1(fRindex1, 0.), N2(realRindex,
1375 G4complex cosPhi;
1376
1377 G4complex u(1., 0.); // unit number 1
1378
1379 G4complex numeratorTE; // E1_perp=1 E1_par
1380 G4complex numeratorTM; // E1_parl=1 E1_per
1381 G4complex denominatorTE, denominatorTM;
1382 G4complex rTM, rTE;
1383
1384 G4MaterialPropertiesTable* MPT = fMaterial1
1385 G4MaterialPropertyVector* ppR = MPT->GetPr
1386 G4MaterialPropertyVector* ppI = MPT->GetPr
1387 if(ppR && ppI)
1388 {
1389 G4double rRindex = ppR->Value(fPhotonMome
1390 G4double iRindex = ppI->Value(fPhotonMome
1391 N1 = G4complex(rRindex, iRi
1392 }
1393
1394 // Following two equations, rTM and rTE, ar
1395 // Optics" written by Fowles
1396 cosPhi = std::sqrt(u - ((std::sin(incidenta
1397 (N1 * N1) / (N2 * N
1398
1399 numeratorTE = N1 * std::cos(incidentangle
1400 denominatorTE = N1 * std::cos(incidentangle
1401 rTE = numeratorTE / denominatorTE
1402
1403 numeratorTM = N2 * std::cos(incidentangle
1404 denominatorTM = N2 * std::cos(incidentangle
1405 rTM = numeratorTM / denominatorTM
1406
1407 // This is my (PG) calculaton for reflectiv
1408 // depending on the fraction of TE and TM p
1409 // when TE polarization, E1_parl=0 and E1_p
1410 // when TM polarization, E1_parl=1 and E1_p
1411
1412 reflectivity_TE = (rTE * conj(rTE)) * (E1_p
1413 (E1_perp * E1_perp + E1_p
1414 reflectivity_TM = (rTM * conj(rTM)) * (E1_p
1415 (E1_perp * E1_perp + E1_p
1416 reflectivity = reflectivity_TE + reflectivi
1417
1418 do
1419 {
1420 if(G4UniformRand() * real(reflectivity) >
1421 {
1422 f_iTE = -1;
1423 }
1424 else
1425 {
1426 f_iTE = 1;
1427 }
1428 if(G4UniformRand() * real(reflectivity) >
1429 {
1430 f_iTM = -1;
1431 }
1432 else
1433 {
1434 f_iTM = 1;
1435 }
1436 // Loop checking, 13-Aug-2015, Peter Gump
1437 } while(f_iTE < 0 && f_iTM < 0);
1438
1439 return real(reflectivity);
1440 }
1441
1442 //....oooOO0OOooo........oooOO0OOooo........o
1443 void G4OpBoundaryProcess::CalculateReflectivi
1444 {
1445 G4double realRindex = fRealRIndexMPV->Value
1446 G4double imaginaryRindex =
1447 fImagRIndexMPV->Value(fPhotonMomentum, id
1448
1449 // calculate FacetNormal
1450 if(fFinish == ground)
1451 {
1452 fFacetNormal = GetFacetNormal(fOldMomentu
1453 }
1454 else
1455 {
1456 fFacetNormal = fGlobalNormal;
1457 }
1458
1459 G4double cost1 = -fOldMomentum * fFacetNorm
1460 if(std::abs(cost1) < 1.0 - fCarTolerance)
1461 {
1462 fSint1 = std::sqrt(1. - cost1 * cost1);
1463 }
1464 else
1465 {
1466 fSint1 = 0.0;
1467 }
1468
1469 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1470 G4double E1_perp, E1_parl;
1471
1472 if(fSint1 > 0.0)
1473 {
1474 A_trans = (fOldMomentum.cross(fFacetNorma
1475 E1_perp = fOldPolarization * A_trans;
1476 E1pp = E1_perp * A_trans;
1477 E1pl = fOldPolarization - E1pp;
1478 E1_parl = E1pl.mag();
1479 }
1480 else
1481 {
1482 A_trans = fOldPolarization;
1483 // Here we Follow Jackson's conventions a
1484 // component = 1 in case of a ray perpend
1485 E1_perp = 0.0;
1486 E1_parl = 1.0;
1487 }
1488
1489 G4double incidentangle = GetIncidentAngle()
1490
1491 // calculate the reflectivity depending on
1492 // polarization and complex refractive
1493 fReflectivity = GetReflectivity(E1_perp, E1
1494 imaginaryRi
1495 }
1496
1497 //....oooOO0OOooo........oooOO0OOooo........o
1498 G4bool G4OpBoundaryProcess::InvokeSD(const G4
1499 {
1500 G4Step aStep = *pStep;
1501 aStep.AddTotalEnergyDeposit(fPhotonMomentum
1502
1503 G4VSensitiveDetector* sd = aStep.GetPostSte
1504 if(sd != nullptr)
1505 return sd->Hit(&aStep);
1506 else
1507 return false;
1508 }
1509
1510 //....oooOO0OOooo........oooOO0OOooo........o
1511 inline void G4OpBoundaryProcess::SetInvokeSD(
1512 {
1513 fInvokeSD = flag;
1514 G4OpticalParameters::Instance()->SetBoundar
1515 }
1516
1517 //....oooOO0OOooo........oooOO0OOooo........o
1518 void G4OpBoundaryProcess::SetVerboseLevel(G4i
1519 {
1520 verboseLevel = verbose;
1521 G4OpticalParameters::Instance()->SetBoundar
1522 }
1523
1524 //....oooOO0OOooo........oooOO0OOooo........o
1525 void G4OpBoundaryProcess::CoatedDielectricDie
1526 {
1527 G4MaterialPropertyVector* pp = nullptr;
1528
1529 G4MaterialPropertiesTable* MPT = fMaterial2
1530 if((pp = MPT->GetProperty(kRINDEX)))
1531 {
1532 fRindex2 = pp->Value(fPhotonMomentum, idx
1533 }
1534
1535 MPT = fOpticalSurface->GetMaterialPropertie
1536 if((pp = MPT->GetProperty(kCOATEDRINDEX)))
1537 {
1538 fCoatedRindex = pp->Value(fPhotonMomentum
1539 }
1540 if(MPT->ConstPropertyExists(kCOATEDTHICKNES
1541 {
1542 fCoatedThickness = MPT->GetConstProperty(
1543 }
1544 if(MPT->ConstPropertyExists(kCOATEDFRUSTRAT
1545 {
1546 fCoatedFrustratedTransmission =
1547 (G4bool)MPT->GetConstProperty(kCOATEDFR
1548 }
1549
1550 G4double sintTL;
1551 G4double wavelength = h_Planck * c_light /
1552 G4double PdotN;
1553 G4double E1_perp, E1_parl;
1554 G4double s1, E2_perp, E2_parl, E2_total, tr
1555 G4double E2_abs, C_parl, C_perp;
1556 G4double alpha;
1557 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1558 //G4bool Inside = false;
1559 //G4bool Swap = false;
1560 G4bool through = false;
1561 G4bool done = false;
1562
1563 do {
1564 if (through)
1565 {
1566 //Swap = !Swap;
1567 through = false;
1568 fGlobalNormal = -fGlobalNormal;
1569 G4SwapPtr(fMaterial1, fMaterial2);
1570 G4SwapObj(&fRindex1, &fRindex2);
1571 }
1572
1573 if(fFinish == polished)
1574 {
1575 fFacetNormal = fGlobalNormal;
1576 }
1577 else
1578 {
1579 fFacetNormal = GetFacetNormal(fOldMomen
1580 }
1581
1582 PdotN = fOldMomentum * fFacetNormal;
1583 G4double cost1 = -PdotN;
1584 G4double sint2, cost2 = 0.;
1585
1586 if (std::abs(cost1) < 1.0 - fCarTolerance
1587 {
1588 fSint1 = std::sqrt(1. - cost1 * cost1);
1589 sint2 = fSint1 * fRindex1 / fRindex2;
1590 sintTL = fSint1 * fRindex1 / fCoatedRin
1591 } else
1592 {
1593 fSint1 = 0.0;
1594 sint2 = 0.0;
1595 sintTL = 0.0;
1596 }
1597
1598 if (fSint1 > 0.0)
1599 {
1600 A_trans = fOldMomentum.cross(fFacetNorm
1601 A_trans = A_trans.unit();
1602 E1_perp = fOldPolarization * A_trans;
1603 E1pp = E1_perp * A_trans;
1604 E1pl = fOldPolarization - E1pp;
1605 E1_parl = E1pl.mag();
1606 }
1607 else
1608 {
1609 A_trans = fOldPolarization;
1610 E1_perp = 0.0;
1611 E1_parl = 1.0;
1612 }
1613
1614 s1 = fRindex1 * cost1;
1615
1616 if (cost1 > 0.0)
1617 {
1618 cost2 = std::sqrt(1. - sint2 * sint2);
1619 }
1620 else
1621 {
1622 cost2 = -std::sqrt(1. - sint2 * sint2);
1623 }
1624
1625 transCoeff = 0.0;
1626
1627 if (sintTL >= 1.0)
1628 { // --> Angle > Angle Limit
1629 //Swap = false;
1630 }
1631 E2_perp = 2. * s1 * E1_perp / (fRindex1 *
1632 E2_parl = 2. * s1 * E1_parl / (fRindex2 *
1633 E2_total = E2_perp * E2_perp + E2_parl *
1634
1635 transCoeff = 1. - GetReflectivityThroughT
1636 sintTL, E1_perp, E1_p
1637 if (!G4BooleanRand(transCoeff))
1638 {
1639 if(verboseLevel > 2)
1640 G4cout << "Reflection from " << fMate
1641 << fMaterial2->GetName() << G4
1642
1643 //Swap = false;
1644
1645 if (sintTL >= 1.0)
1646 {
1647 fStatus = TotalInternalReflection;
1648 }
1649 else
1650 {
1651 fStatus = CoatedDielectricReflection;
1652 }
1653
1654 PdotN = fOldMomentum * fFacetNormal;
1655 fNewMomentum = fOldMomentum - (2. * Pdo
1656
1657 if (fSint1 > 0.0) { // incident ray o
1658
1659 E2_parl = fRindex2 * E2_parl / fRinde
1660 E2_perp = E2_perp - E1_perp;
1661 E2_total = E2_perp * E2_perp + E2_par
1662 A_paral = fNewMomentum.cross(A_trans)
1663 A_paral = A_paral.unit();
1664 E2_abs = std::sqrt(E2_total);
1665 C_parl = E2_parl / E2_abs;
1666 C_perp = E2_perp / E2_abs;
1667
1668 fNewPolarization = C_parl * A_paral +
1669
1670 }
1671 else
1672 { // incident ray perpend
1673 if (fRindex2 > fRindex1)
1674 {
1675 fNewPolarization = -fOldPolarizatio
1676 }
1677 else
1678 {
1679 fNewPolarization = fOldPolarization
1680 }
1681 }
1682
1683 } else { // photon gets transmitted
1684 if (verboseLevel > 2)
1685 G4cout << "Transmission from " << fMa
1686 << fMaterial2->GetName() << G4
1687
1688 //Inside = !Inside;
1689 through = true;
1690
1691 if (fEfficiency > 0.)
1692 {
1693 DoAbsorption();
1694 return;
1695 }
1696 else
1697 {
1698 if (sintTL >= 1.0)
1699 {
1700 fStatus = CoatedDielectricFrustrate
1701 }
1702 else
1703 {
1704 fStatus = CoatedDielectricRefractio
1705 }
1706
1707 if (fSint1 > 0.0) { // incident
1708
1709 alpha = cost1 - cost2 * (fRindex2 /
1710 fNewMomentum = fOldMomentum + alpha
1711 fNewMomentum = fNewMomentum.unit();
1712 A_paral = fNewMomentum.cross(A_tran
1713 A_paral = A_paral.unit();
1714 E2_abs = std::sqrt(E2_total);
1715 C_parl = E2_parl / E2_abs;
1716 C_perp = E2_perp / E2_abs;
1717
1718 fNewPolarization = C_parl * A_paral
1719
1720 }
1721 else
1722 { // incident ray pe
1723 fNewMomentum = fOldMomentum;
1724 fNewPolarization = fOldPolarization
1725 }
1726 }
1727 }
1728
1729 fOldMomentum = fNewMomentum.unit();
1730 fOldPolarization = fNewPolarization.unit(
1731 if ((fStatus == CoatedDielectricFrustrate
1732 (fStatus == CoatedDielectricRefractio
1733 {
1734 done = (fNewMomentum * fGlobalNormal <=
1735 }
1736 else
1737 {
1738 done = (fNewMomentum * fGlobalNormal >=
1739 }
1740
1741 } while (!done);
1742 }
1743
1744 //....oooOO0OOooo........oooOO0OOooo........o
1745 G4double G4OpBoundaryProcess::GetReflectivity
1746 G4double E1_perp,
1747 G4double E1_parl,
1748 G4double wavelength, G4dou
1749 G4complex Reflectivity, Reflectivity_TE, Re
1750 G4double gammaTL, costTL;
1751
1752 G4complex i(0, 1);
1753 G4complex rTM, rTE;
1754 G4complex r1toTL, rTLto2;
1755 G4double k0 = 2 * pi / wavelength;
1756
1757 // Angle > Angle limit
1758 if (sinTL >= 1.0) {
1759 if (fCoatedFrustratedTransmission) { //Fr
1760
1761 if (cost1 > 0.0)
1762 {
1763 gammaTL = std::sqrt(fRindex1 * fRinde
1764 fCoatedRindex * fCoatedRin
1765 }
1766 else
1767 {
1768 gammaTL = -std::sqrt(fRindex1 * fRind
1769 fCoatedRindex * fCoatedRin
1770 }
1771
1772 // TE
1773 r1toTL = (fRindex1 * cost1 - i * gammaT
1774 rTLto2 = (i * gammaTL - fRindex2 * cost
1775 if (cost1 != 0.0)
1776 {
1777 rTE = (r1toTL + rTLto2 * std::exp(-2
1778 (1.0 + r1toTL * rTLto2 * std
1779 }
1780 // TM
1781 r1toTL = (fRindex1 * i * gammaTL - fCoa
1782 (fRindex1 * i * gammaTL + f
1783 rTLto2 = (fCoatedRindex * fCoatedRindex
1784 (fCoatedRindex * fCoatedRin
1785 if (cost1 != 0.0)
1786 {
1787 rTM = (r1toTL + rTLto2 * std::exp(-2
1788 (1.0 + r1toTL * rTLto2 * std
1789 }
1790 }
1791 else
1792 { //Total reflection
1793 return(1.);
1794 }
1795 }
1796
1797 // Angle <= Angle limit
1798 else //if (sinTL < 1.0)
1799 {
1800 if (cost1 > 0.0)
1801 {
1802 costTL = std::sqrt(1. - sinTL * sinTL);
1803 }
1804 else
1805 {
1806 costTL = -std::sqrt(1. - sinTL * sinTL)
1807 }
1808 // TE
1809 r1toTL = (fRindex1 * cost1 - fCoatedRinde
1810 rTLto2 = (fCoatedRindex * costTL - fRinde
1811 if (cost1 != 0.0)
1812 {
1813 rTE = (r1toTL + rTLto2 * std::exp(2.0 *
1814 (1.0 + r1toTL * rTLto2 * std::exp
1815 }
1816 // TM
1817 r1toTL = (fRindex1 * costTL - fCoatedRind
1818 rTLto2 = (fCoatedRindex * cost2 - fRindex
1819 if (cost1 != 0.0)
1820 {
1821 rTM = (r1toTL + rTLto2 * std::exp(2.0 *
1822 (1.0 + r1toTL * rTLto2 * std::exp
1823 }
1824 }
1825
1826 Reflectivity_TE = (rTE * conj(rTE)) * (E1_p
1827 Reflectivity_TM = (rTM * conj(rTM)) * (E1_p
1828 Reflectivity = Reflectivity_TE + Reflectivi
1829
1830 return real(Reflectivity);
1831 }
1832