| 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 // Optical Photon Boundary Process Class Implementation
28 ////////////////////////////////////////////////////////////////////////
29 //
30 // File: G4OpBoundaryProcess.cc
31 // Description: Discrete Process -- reflection/refraction at
32 // optical interfaces
33 // Version: 1.1
34 // Created: 1997-06-18
35 // Modified: 1998-05-25 - Correct parallel component of polarization
36 // (thanks to: Stefano Magni + Giovanni Pieri)
37 // 1998-05-28 - NULL Rindex pointer before reuse
38 // (thanks to: Stefano Magni)
39 // 1998-06-11 - delete *sint1 in oblique reflection
40 // (thanks to: Giovanni Pieri)
41 // 1998-06-19 - move from GetLocalExitNormal() to the new
42 // method: GetLocalExitNormal(&valid) to get
43 // the surface normal in all cases
44 // 1998-11-07 - NULL OpticalSurface pointer before use
45 // comparison not sharp for: std::abs(cost1) < 1.0
46 // remove sin1, sin2 in lines 556,567
47 // (thanks to Stefano Magni)
48 // 1999-10-10 - Accommodate changes done in DoAbsorption by
49 // changing logic in DielectricMetal
50 // 2001-10-18 - avoid Linux (gcc-2.95.2) warning about variables
51 // might be used uninitialized in this function
52 // moved E2_perp, E2_parl and E2_total out of 'if'
53 // 2003-11-27 - Modified line 168-9 to reflect changes made to
54 // G4OpticalSurface class ( by Fan Lei)
55 // 2004-02-02 - Set theStatus = Undefined at start of DoIt
56 // 2005-07-28 - add G4ProcessType to constructor
57 // 2006-11-04 - add capability of calculating the reflectivity
58 // off a metal surface by way of a complex index
59 // of refraction - Thanks to Sehwook Lee and John
60 // Hauptman (Dept. of Physics - Iowa State Univ.)
61 // 2009-11-10 - add capability of simulating surface reflections
62 // with Look-Up-Tables (LUT) containing measured
63 // optical reflectance for a variety of surface
64 // treatments - Thanks to Martin Janecek and
65 // William Moses (Lawrence Berkeley National Lab.)
66 // 2013-06-01 - add the capability of simulating the transmission
67 // of a dichronic filter
68 // 2017-02-24 - add capability of simulating surface reflections
69 // with Look-Up-Tables (LUT) developed in DAVIS
70 //
71 // Author: Peter Gumplinger
72 // adopted from work by Werner Keil - April 2/96
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........oooOO0OOooo........oooOO0OOooo......
91 G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
92 G4ProcessType ptype)
93 : G4VDiscreteProcess(processName, ptype)
94 {
95 Initialise();
96
97 if(verboseLevel > 0)
98 {
99 G4cout << GetProcessName() << " is created " << G4endl;
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::GetInstance()->GetSurfaceTolerance();
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........oooOO0OOooo........oooOO0OOooo......
129 G4OpBoundaryProcess::~G4OpBoundaryProcess() = default;
130
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 void G4OpBoundaryProcess::PreparePhysicsTable(const G4ParticleDefinition&)
133 {
134 Initialise();
135 }
136
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
138 void G4OpBoundaryProcess::Initialise()
139 {
140 G4OpticalParameters* params = G4OpticalParameters::Instance();
141 SetInvokeSD(params->GetBoundaryInvokeSD());
142 SetVerboseLevel(params->GetBoundaryVerboseLevel());
143 }
144
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146 G4VParticleChange* G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack,
147 const G4Step& aStep)
148 {
149 fStatus = Undefined;
150 aParticleChange.Initialize(aTrack);
151 aParticleChange.ProposeVelocity(aTrack.GetVelocity());
152
153 // Get hyperStep from G4ParallelWorldProcess
154 // NOTE: PostSetpDoIt of this process to be invoked after
155 // G4ParallelWorldProcess!
156 const G4Step* pStep = &aStep;
157 const G4Step* hStep = G4ParallelWorldProcess::GetHyperStep();
158 if(hStep != nullptr)
159 pStep = hStep;
160
161 if(pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
162 {
163 fMaterial1 = pStep->GetPreStepPoint()->GetMaterial();
164 fMaterial2 = pStep->GetPostStepPoint()->GetMaterial();
165 }
166 else
167 {
168 fStatus = NotAtBoundary;
169 if(verboseLevel > 1)
170 BoundaryProcessVerbose();
171 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
172 }
173
174 G4VPhysicalVolume* thePrePV = pStep->GetPreStepPoint()->GetPhysicalVolume();
175 G4VPhysicalVolume* thePostPV = pStep->GetPostStepPoint()->GetPhysicalVolume();
176
177 if(verboseLevel > 1)
178 {
179 G4cout << " Photon at Boundary! " << G4endl;
180 if(thePrePV != nullptr)
181 G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
182 if(thePostPV != nullptr)
183 G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
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 = nullptr;
194 G4MaterialPropertiesTable* aMPT = fMaterial2->GetMaterialPropertiesTable();
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_groupvel));
204 }
205 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
206 }
207 else if (stepLength <= 10.*fCarTolerance && fNumSmallStepWarnings < 10)
208 { // see bug 2510
209 ++fNumSmallStepWarnings;
210 if(verboseLevel > 0)
211 {
212 G4ExceptionDescription ed;
213 ed << "G4OpBoundaryProcess: "
214 << "Opticalphoton step length: " << stepLength/mm << " mm." << G4endl
215 << "This is larger than the threshold " << fCarTolerance/mm << " mm "
216 "to set status StepTooSmall." << G4endl
217 << "Boundary scattering may be incorrect. ";
218 if(fNumSmallStepWarnings == 10)
219 {
220 ed << G4endl << "*** Step size warnings stopped.";
221 }
222 G4Exception("G4OpBoundaryProcess", "OpBoun06", JustWarning, ed, "");
223 }
224 }
225
226 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
227
228 fPhotonMomentum = aParticle->GetTotalMomentum();
229 fOldMomentum = aParticle->GetMomentumDirection();
230 fOldPolarization = aParticle->GetPolarization();
231
232 if(verboseLevel > 1)
233 {
234 G4cout << " Old Momentum Direction: " << fOldMomentum << G4endl
235 << " Old Polarization: " << fOldPolarization << G4endl;
236 }
237
238 G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
239 G4bool valid;
240
241 // ID of Navigator which limits step
242 G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
243 auto iNav = G4TransportationManager::GetTransportationManager()
244 ->GetActiveNavigatorsIterator();
245 fGlobalNormal = (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint, &valid);
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 returned an invalid normal" << G4endl;
256 G4Exception(
257 "G4OpBoundaryProcess::PostStepDoIt", "OpBoun01", EventMustBeAborted, ed,
258 "Invalid Surface Normal - Geometry must return valid surface normal");
259 }
260
261 if(fOldMomentum * fGlobalNormal > 0.0)
262 {
263 #ifdef G4OPTICAL_DEBUG
264 G4ExceptionDescription ed;
265 ed << " G4OpBoundaryProcess/PostStepDoIt(): fGlobalNormal points in a "
266 "wrong direction. "
267 << G4endl
268 << " The momentum of the photon arriving at interface (oldMomentum)"
269 << " must exit the volume cross in the step. " << G4endl
270 << " So it MUST have dot < 0 with the normal that Exits the new "
271 "volume (globalNormal)."
272 << G4endl << " >> The dot product of oldMomentum and global Normal is "
273 << fOldMomentum * fGlobalNormal << G4endl
274 << " Old Momentum (during step) = " << fOldMomentum << G4endl
275 << " Global Normal (Exiting New Vol) = " << fGlobalNormal << G4endl
276 << G4endl;
277 G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
278 EventMustBeAborted, // Or JustWarning to see if it happens
279 // repeatedly on one ray
280 ed,
281 "Invalid Surface Normal - Geometry must return valid surface "
282 "normal pointing in the right direction");
283 #else
284 fGlobalNormal = -fGlobalNormal;
285 #endif
286 }
287
288 G4MaterialPropertyVector* rIndexMPV = nullptr;
289 G4MaterialPropertiesTable* MPT = fMaterial1->GetMaterialPropertiesTable();
290 if(MPT != nullptr)
291 {
292 rIndexMPV = MPT->GetProperty(kRINDEX);
293 }
294 if(rIndexMPV != nullptr)
295 {
296 fRindex1 = rIndexMPV->Value(fPhotonMomentum, idx_rindex1);
297 }
298 else
299 {
300 fStatus = NoRINDEX;
301 if(verboseLevel > 1)
302 BoundaryProcessVerbose();
303 aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
304 aParticleChange.ProposeTrackStatus(fStopAndKill);
305 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
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(thePrePV, thePostPV);
321 if(surface == nullptr)
322 {
323 if(thePostPV->GetMotherLogical() == thePrePV->GetLogicalVolume())
324 {
325 surface = G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
326 if(surface == nullptr)
327 {
328 surface =
329 G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
330 }
331 }
332 else
333 {
334 surface = G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
335 if(surface == nullptr)
336 {
337 surface =
338 G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
339 }
340 }
341 }
342
343 if(surface != nullptr)
344 {
345 fOpticalSurface =
346 dynamic_cast<G4OpticalSurface*>(surface->GetSurfaceProperty());
347 }
348 if(fOpticalSurface != nullptr)
349 {
350 type = fOpticalSurface->GetType();
351 fModel = fOpticalSurface->GetModel();
352 fFinish = fOpticalSurface->GetFinish();
353
354 G4MaterialPropertiesTable* sMPT =
355 fOpticalSurface->GetMaterialPropertiesTable();
356 if(sMPT != nullptr)
357 {
358 if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
359 {
360 rIndexMPV = sMPT->GetProperty(kRINDEX);
361 if(rIndexMPV != nullptr)
362 {
363 fRindex2 = rIndexMPV->Value(fPhotonMomentum, idx_rindex_surface);
364 }
365 else
366 {
367 fStatus = NoRINDEX;
368 if(verboseLevel > 1)
369 BoundaryProcessVerbose();
370 aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
371 aParticleChange.ProposeTrackStatus(fStopAndKill);
372 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
373 }
374 }
375
376 fRealRIndexMPV = sMPT->GetProperty(kREALRINDEX);
377 fImagRIndexMPV = sMPT->GetProperty(kIMAGINARYRINDEX);
378 f_iTE = f_iTM = 1;
379
380 G4MaterialPropertyVector* pp;
381 if((pp = sMPT->GetProperty(kREFLECTIVITY)))
382 {
383 fReflectivity = pp->Value(fPhotonMomentum, idx_reflect);
384 }
385 else if(fRealRIndexMPV && fImagRIndexMPV)
386 {
387 CalculateReflectivity();
388 }
389
390 if((pp = sMPT->GetProperty(kEFFICIENCY)))
391 {
392 fEfficiency = pp->Value(fPhotonMomentum, idx_eff);
393 }
394 if((pp = sMPT->GetProperty(kTRANSMITTANCE)))
395 {
396 fTransmittance = pp->Value(fPhotonMomentum, idx_trans);
397 }
398 if(sMPT->ConstPropertyExists(kSURFACEROUGHNESS))
399 {
400 fSurfaceRoughness = sMPT->GetConstProperty(kSURFACEROUGHNESS);
401 }
402
403 if(fModel == unified)
404 {
405 fProb_sl = (pp = sMPT->GetProperty(kSPECULARLOBECONSTANT))
406 ? pp->Value(fPhotonMomentum, idx_lobe)
407 : 0.;
408 fProb_ss = (pp = sMPT->GetProperty(kSPECULARSPIKECONSTANT))
409 ? pp->Value(fPhotonMomentum, idx_spike)
410 : 0.;
411 fProb_bs = (pp = sMPT->GetProperty(kBACKSCATTERCONSTANT))
412 ? pp->Value(fPhotonMomentum, idx_back)
413 : 0.;
414 }
415 } // end of if(sMPT)
416 else if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
417 {
418 aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
419 aParticleChange.ProposeTrackStatus(fStopAndKill);
420 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
421 }
422 } // end of if(fOpticalSurface)
423
424 // DIELECTRIC-DIELECTRIC
425 if(type == dielectric_dielectric)
426 {
427 if(fFinish == polished || fFinish == ground)
428 {
429 if(fMaterial1 == fMaterial2)
430 {
431 fStatus = SameMaterial;
432 if(verboseLevel > 1)
433 BoundaryProcessVerbose();
434 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
435 }
436 MPT = fMaterial2->GetMaterialPropertiesTable();
437 rIndexMPV = nullptr;
438 if(MPT != nullptr)
439 {
440 rIndexMPV = MPT->GetProperty(kRINDEX);
441 }
442 if(rIndexMPV != nullptr)
443 {
444 fRindex2 = rIndexMPV->Value(fPhotonMomentum, idx_rindex2);
445 }
446 else
447 {
448 fStatus = NoRINDEX;
449 if(verboseLevel > 1)
450 BoundaryProcessVerbose();
451 aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
452 aParticleChange.ProposeTrackStatus(fStopAndKill);
453 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
454 }
455 }
456 if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
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 boundary type." << G4endl;
520 if(fNumBdryTypeWarnings == 10)
521 {
522 ed << "** Boundary type warnings stopped." << G4endl;
523 }
524 G4Exception("G4OpBoundaryProcess", "OpBoun04", JustWarning, ed);
525 }
526 }
527 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
528 }
529
530 fNewMomentum = fNewMomentum.unit();
531 fNewPolarization = fNewPolarization.unit();
532
533 if(verboseLevel > 1)
534 {
535 G4cout << " New Momentum Direction: " << fNewMomentum << G4endl
536 << " New Polarization: " << fNewPolarization << G4endl;
537 BoundaryProcessVerbose();
538 }
539
540 aParticleChange.ProposeMomentumDirection(fNewMomentum);
541 aParticleChange.ProposePolarization(fNewPolarization);
542
543 if(fStatus == FresnelRefraction || fStatus == Transmission)
544 {
545 // not all surface types check that fMaterial2 has an MPT
546 G4MaterialPropertiesTable* aMPT = fMaterial2->GetMaterialPropertiesTable();
547 G4MaterialPropertyVector* groupvel = nullptr;
548 if(aMPT != nullptr)
549 {
550 groupvel = aMPT->GetProperty(kGROUPVEL);
551 }
552 if(groupvel != nullptr)
553 {
554 aParticleChange.ProposeVelocity(
555 groupvel->Value(fPhotonMomentum, idx_groupvel));
556 }
557 }
558
559 if(fStatus == Detection && fInvokeSD)
560 InvokeSD(pStep);
561 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
562 }
563
564 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
565 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
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 == PolishedLumirrorAirReflection)
587 G4cout << "PolishedLumirrorAirReflection";
588 else if(fStatus == PolishedLumirrorGlueReflection)
589 G4cout << "PolishedLumirrorGlueReflection";
590 else if(fStatus == PolishedAirReflection)
591 G4cout << "PolishedAirReflection";
592 else if(fStatus == PolishedTeflonAirReflection)
593 G4cout << "PolishedTeflonAirReflection";
594 else if(fStatus == PolishedTiOAirReflection)
595 G4cout << "PolishedTiOAirReflection";
596 else if(fStatus == PolishedTyvekAirReflection)
597 G4cout << "PolishedTyvekAirReflection";
598 else if(fStatus == PolishedVM2000AirReflection)
599 G4cout << "PolishedVM2000AirReflection";
600 else if(fStatus == PolishedVM2000GlueReflection)
601 G4cout << "PolishedVM2000GlueReflection";
602 else if(fStatus == EtchedLumirrorAirReflection)
603 G4cout << "EtchedLumirrorAirReflection";
604 else if(fStatus == EtchedLumirrorGlueReflection)
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 == EtchedVM2000GlueReflection)
617 G4cout << "EtchedVM2000GlueReflection";
618 else if(fStatus == GroundLumirrorAirReflection)
619 G4cout << "GroundLumirrorAirReflection";
620 else if(fStatus == GroundLumirrorGlueReflection)
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 == GroundVM2000GlueReflection)
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 == CoatedDielectricReflection)
649 G4cout << "Coated Dielectric Reflection";
650 else if(fStatus == CoatedDielectricRefraction)
651 G4cout << "Coated Dielectric Refraction";
652 else if(fStatus == CoatedDielectricFrustratedTransmission)
653 G4cout << "Coated Dielectric Frustrated Transmission";
654
655 G4cout << " ***" << G4endl;
656 }
657
658 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
659 G4ThreeVector G4OpBoundaryProcess::GetFacetNormal(
660 const G4ThreeVector& momentum, const G4ThreeVector& normal) const
661 {
662 G4ThreeVector facetNormal;
663 if(fModel == unified || fModel == LUT || fModel == DAVIS)
664 {
665 /* This function codes alpha to a random value taken from the
666 distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
667 for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha) is a
668 gaussian distribution with mean 0 and standard deviation sigma_alpha. */
669
670 G4double sigma_alpha = 0.0;
671 if(fOpticalSurface)
672 sigma_alpha = fOpticalSurface->GetSigmaAlpha();
673 if(sigma_alpha == 0.0)
674 {
675 return normal;
676 }
677
678 G4double f_max = std::min(1.0, 4. * sigma_alpha);
679 G4double alpha, phi, sinAlpha;
680
681 do
682 { // Loop checking, 13-Aug-2015, Peter Gumplinger
683 do
684 { // Loop checking, 13-Aug-2015, Peter Gumplinger
685 alpha = G4RandGauss::shoot(0.0, sigma_alpha);
686 sinAlpha = std::sin(alpha);
687 } while(G4UniformRand() * f_max > sinAlpha || alpha >= halfpi);
688
689 phi = G4UniformRand() * twopi;
690 facetNormal.set(sinAlpha * std::cos(phi), sinAlpha * std::sin(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 Gumplinger
704 G4ThreeVector smear;
705 do
706 { // Loop checking, 13-Aug-2015, Peter Gumplinger
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) * smear;
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........oooOO0OOooo........oooOO0OOooo......
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 == polished)
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 || !fImagRIndexMPV)
784 {
785 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
786 }
787 // else
788 // case of complex rindex needs to be implemented
789 }
790 fNewMomentum =
791 fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
792
793 if(f_iTE > 0 && f_iTM > 0)
794 {
795 fNewPolarization =
796 -fOldPolarization +
797 (2. * fOldPolarization * fFacetNormal * fFacetNormal);
798 }
799 else if(f_iTE > 0)
800 {
801 A_trans = (fSint1 > 0.0) ? fOldMomentum.cross(fFacetNormal).unit()
802 : fOldPolarization;
803 fNewPolarization = -A_trans;
804 }
805 else if(f_iTM > 0)
806 {
807 fNewPolarization =
808 -fNewMomentum.cross(A_trans).unit(); // = -A_paral
809 }
810 }
811 }
812 fOldMomentum = fNewMomentum;
813 fOldPolarization = fNewPolarization;
814 }
815 // Loop checking, 13-Aug-2015, Peter Gumplinger
816 } while(fNewMomentum * fGlobalNormal < 0.0);
817 }
818
819 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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(groundbackpainted)));
828
829 G4int thetaIndexMax = fOpticalSurface->GetThetaIndexMax();
830 G4int phiIndexMax = fOpticalSurface->GetPhiIndexMax();
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 Photon Momentum
854 G4double anglePhotonToNormal = fOldMomentum.angle(-fGlobalNormal);
855 // Round to closest integer: LBNL model array has 91 values
856 G4int angleIncident = (G4int)std::lrint(anglePhotonToNormal / CLHEP::deg);
857
858 // Take random angles THETA and PHI,
859 // and see if below Probability - if not - Redo
860 do
861 {
862 thetaIndex = (G4int)G4RandFlat::shootInt(thetaIndexMax - 1);
863 phiIndex = (G4int)G4RandFlat::shootInt(phiIndexMax - 1);
864 // Find probability with the new indeces from LUT
865 angularDistVal = fOpticalSurface->GetAngularDistributionValue(
866 angleIncident, thetaIndex, phiIndex);
867 // Loop checking, 13-Aug-2015, Peter Gumplinger
868 } while(!G4BooleanRand(angularDistVal));
869
870 thetaRad = G4double(-90 + 4 * thetaIndex) * pi / 180.;
871 phiRad = G4double(-90 + 5 * phiIndex) * pi / 180.;
872 // Rotate Photon Momentum in Theta, then in Phi
873 fNewMomentum = -fOldMomentum;
874
875 perpVectorTheta = fNewMomentum.cross(fGlobalNormal);
876 if(perpVectorTheta.mag() < fCarTolerance)
877 {
878 perpVectorTheta = fNewMomentum.orthogonal();
879 }
880 fNewMomentum =
881 fNewMomentum.rotate(anglePhotonToNormal - thetaRad, perpVectorTheta);
882 perpVectorPhi = perpVectorTheta.cross(fNewMomentum);
883 fNewMomentum = fNewMomentum.rotate(-phiRad, perpVectorPhi);
884
885 // Rotate Polarization too:
886 fFacetNormal = (fNewMomentum - fOldMomentum).unit();
887 fNewPolarization = -fOldPolarization +
888 (2. * fOldPolarization * fFacetNormal * fFacetNormal);
889 }
890 // Loop checking, 13-Aug-2015, Peter Gumplinger
891 } while(fNewMomentum * fGlobalNormal <= 0.0);
892 }
893
894 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
895 void G4OpBoundaryProcess::DielectricLUTDAVIS()
896 {
897 G4int angindex, random, angleIncident;
898 G4double reflectivityValue, elevation, azimuth;
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(-fGlobalNormal);
910
911 // Davis model has 90 reflection bins: round down
912 // don't allow angleIncident to be 90 for anglePhotonToNormal close to 90
913 angleIncident = std::min(
914 static_cast<G4int>(std::floor(anglePhotonToNormal / CLHEP::deg)), 89);
915 reflectivityValue = fOpticalSurface->GetReflectivityLUTValue(angleIncident);
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(1, lutbin + 1);
937 angindex =
938 (((random * 2) - 1)) + angleIncident * lutbin * 2 + 3640000;
939
940 azimuth =
941 fOpticalSurface->GetAngularDistributionValueLUT(angindex - 1);
942 elevation = fOpticalSurface->GetAngularDistributionValueLUT(angindex);
943 } while(elevation == 0. && azimuth == 0.);
944
945 sinEl = std::sin(elevation);
946 vNorm = (fGlobalNormal.cross(fOldMomentum)).unit();
947 u = vNorm.cross(fGlobalNormal) * (sinEl * std::cos(azimuth));
948 vNorm *= (sinEl * std::sin(azimuth));
949 // fGlobalNormal shouldn't be modified here
950 w = (fGlobalNormal *= std::cos(elevation));
951 fNewMomentum = u + vNorm + w;
952
953 // Rotate Polarization too:
954 fFacetNormal = (fNewMomentum - fOldMomentum).unit();
955 fNewPolarization = -fOldPolarization + (2. * fOldPolarization *
956 fFacetNormal * fFacetNormal);
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(1, lutbin + 1);
972 angindex = (((random * 2) - 1)) + (angleIncident - 1) * lutbin * 2;
973
974 azimuth = fOpticalSurface->GetAngularDistributionValueLUT(angindex - 1);
975 elevation = fOpticalSurface->GetAngularDistributionValueLUT(angindex);
976 } while(elevation == 0. && azimuth == 0.);
977
978 sinEl = std::sin(elevation);
979 vNorm = (fGlobalNormal.cross(fOldMomentum)).unit();
980 u = vNorm.cross(fGlobalNormal) * (sinEl * std::cos(azimuth));
981 vNorm *= (sinEl * std::sin(azimuth));
982 // fGlobalNormal shouldn't be modified here
983 w = (fGlobalNormal *= std::cos(elevation));
984
985 fNewMomentum = u + vNorm + w;
986
987 // Rotate Polarization too: (needs revision)
988 fNewPolarization = fOldPolarization;
989 }
990 } while(fNewMomentum * fGlobalNormal <= 0.0);
991 }
992
993 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
994 void G4OpBoundaryProcess::DielectricDichroic()
995 {
996 // Calculate Angle between Normal and Photon Momentum
997 G4double anglePhotonToNormal = fOldMomentum.angle(-fGlobalNormal);
998
999 // Round it to closest integer
1000 G4double angleIncident = std::floor(180. / pi * anglePhotonToNormal + 0.5);
1001
1002 if(!fDichroicVector)
1003 {
1004 if(fOpticalSurface)
1005 fDichroicVector = fOpticalSurface->GetDichroicVector();
1006 }
1007
1008 if(fDichroicVector)
1009 {
1010 G4double wavelength = h_Planck * c_light / fPhotonMomentum;
1011 fTransmittance = fDichroicVector->Value(wavelength / nm, angleIncident,
1012 idx_dichroicX, idx_dichroicY) *
1013 perCent;
1014 // G4cout << "wavelength: " << std::floor(wavelength/nm)
1015 // << "nm" << G4endl;
1016 // G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
1017 // G4cout << "Transmittance: "
1018 // << std::floor(fTransmittance/perCent) << "%" << G4endl;
1019 }
1020 else
1021 {
1022 G4ExceptionDescription ed;
1023 ed << " G4OpBoundaryProcess/DielectricDichroic(): "
1024 << " The dichroic surface has no G4Physics2DVector" << G4endl;
1025 G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
1026 FatalException, ed,
1027 "A dichroic surface must have an associated G4Physics2DVector");
1028 }
1029
1030 if(!G4BooleanRand(fTransmittance))
1031 { // Not transmitted, so reflect
1032 if(fModel == glisur || fFinish == polished)
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(fOldMomentum, fGlobalNormal);
1056 }
1057 PdotN = fOldMomentum * fFacetNormal;
1058 fNewMomentum = fOldMomentum - (2. * PdotN) * fFacetNormal;
1059 // Loop checking, 13-Aug-2015, Peter Gumplinger
1060 } while(fNewMomentum * fGlobalNormal <= 0.0);
1061
1062 EdotN = fOldPolarization * fFacetNormal;
1063 fNewPolarization = -fOldPolarization + (2. * EdotN) * fFacetNormal;
1064 }
1065 }
1066 }
1067 else
1068 {
1069 fStatus = Dichroic;
1070 fNewMomentum = fOldMomentum;
1071 fNewPolarization = fOldPolarization;
1072 }
1073 }
1074
1075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1076 void G4OpBoundaryProcess::DielectricDielectric()
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(fOldMomentum, fGlobalNormal);
1088 }
1089 G4double cost1 = -fOldMomentum * fFacetNormal;
1090 G4double cost2 = 0.;
1091 G4double sint2 = 0.;
1092
1093 G4bool surfaceRoughnessCriterionPass = true;
1094 if(fSurfaceRoughness != 0. && fRindex1 > fRindex2)
1095 {
1096 G4double wavelength = h_Planck * c_light / fPhotonMomentum;
1097 G4double surfaceRoughnessCriterion = std::exp(-std::pow(
1098 (4. * pi * fSurfaceRoughness * fRindex1 * cost1 / wavelength), 2));
1099 surfaceRoughnessCriterionPass = G4BooleanRand(surfaceRoughnessCriterion);
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, transCoeff;
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(fOldMomentum, fGlobalNormal);
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; // *** Snell's Law ***
1138 // this isn't a sine as we might expect from the name; can be > 1
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 != polished)
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 * fFacetNormal * fFacetNormal;
1169 fNewPolarization = -fOldPolarization + (2. * fOldPolarization *
1170 fFacetNormal * fFacetNormal);
1171 }
1172 }
1173 // NOT TIR
1174 else if(sint2 < 1.0)
1175 {
1176 // Calculate amplitude for transmission (Q = P x N)
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(fFacetNormal)).unit();
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 conventions and set the parallel
1198 // component = 1 in case of a ray perpendicular to the surface
1199 E1_perp = 0.0;
1200 E1_parl = 1.0;
1201 }
1202
1203 s1 = fRindex1 * cost1;
1204 E2_perp = 2. * s1 * E1_perp / (fRindex1 * cost1 + fRindex2 * cost2);
1205 E2_parl = 2. * s1 * E1_parl / (fRindex2 * cost1 + fRindex1 * cost2);
1206 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1207 s2 = fRindex2 * cost2 * E2_total;
1208
1209 // D.Sawkey, 24 May 24
1210 // Transmittance has already been taken into account in PostStepDoIt.
1211 // For e.g. specular surfaces, the ratio of Fresnel refraction to
1212 // reflection should be given by the math, not material property
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 != polished)
1231 ChooseReflection();
1232 if(fStatus == LambertianReflection)
1233 {
1234 DoReflection();
1235 }
1236 else if(fStatus == BackScattering)
1237 {
1238 fNewMomentum = -fOldMomentum;
1239 fNewPolarization = -fOldPolarization;
1240 }
1241 else
1242 {
1243 fNewMomentum =
1244 fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
1245 if(fSint1 > 0.0)
1246 { // incident ray oblique
1247 E2_parl = fRindex2 * E2_parl / fRindex1 - E1_parl;
1248 E2_perp = E2_perp - E1_perp;
1249 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1250 A_paral = (fNewMomentum.cross(A_trans)).unit();
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_paral + C_perp * A_trans;
1256 }
1257 else
1258 { // incident ray perpendicular
1259 if(fRindex2 > fRindex1)
1260 {
1261 fNewPolarization = -fOldPolarization;
1262 }
1263 else
1264 {
1265 fNewPolarization = fOldPolarization;
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 * (fRindex2 / fRindex1);
1280 fNewMomentum = (fOldMomentum + alpha * fFacetNormal).unit();
1281 A_paral = (fNewMomentum.cross(A_trans)).unit();
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 + C_perp * A_trans;
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 <= 0.0);
1302 }
1303 else
1304 {
1305 done = (fNewMomentum * fGlobalNormal >= -fCarTolerance);
1306 }
1307 // Loop checking, 13-Aug-2015, Peter Gumplinger
1308 } while(!done);
1309
1310 if(inside && !swap)
1311 {
1312 if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
1313 {
1314 G4double rand = G4UniformRand();
1315 if(rand > fReflectivity + fTransmittance)
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........oooOO0OOooo........oooOO0OOooo......
1352 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track&, G4double,
1353 G4ForceCondition* condition)
1354 {
1355 *condition = Forced;
1356 return DBL_MAX;
1357 }
1358
1359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1360 G4double G4OpBoundaryProcess::GetIncidentAngle()
1361 {
1362 return pi - std::acos(fOldMomentum * fFacetNormal /
1363 (fOldMomentum.mag() * fFacetNormal.mag()));
1364 }
1365
1366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1367 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1368 G4double E1_parl,
1369 G4double incidentangle,
1370 G4double realRindex,
1371 G4double imaginaryRindex)
1372 {
1373 G4complex reflectivity, reflectivity_TE, reflectivity_TM;
1374 G4complex N1(fRindex1, 0.), N2(realRindex, imaginaryRindex);
1375 G4complex cosPhi;
1376
1377 G4complex u(1., 0.); // unit number 1
1378
1379 G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1380 G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1381 G4complex denominatorTE, denominatorTM;
1382 G4complex rTM, rTE;
1383
1384 G4MaterialPropertiesTable* MPT = fMaterial1->GetMaterialPropertiesTable();
1385 G4MaterialPropertyVector* ppR = MPT->GetProperty(kREALRINDEX);
1386 G4MaterialPropertyVector* ppI = MPT->GetProperty(kIMAGINARYRINDEX);
1387 if(ppR && ppI)
1388 {
1389 G4double rRindex = ppR->Value(fPhotonMomentum, idx_rrindex);
1390 G4double iRindex = ppI->Value(fPhotonMomentum, idx_irindex);
1391 N1 = G4complex(rRindex, iRindex);
1392 }
1393
1394 // Following two equations, rTM and rTE, are from: "Introduction To Modern
1395 // Optics" written by Fowles
1396 cosPhi = std::sqrt(u - ((std::sin(incidentangle) * std::sin(incidentangle)) *
1397 (N1 * N1) / (N2 * N2)));
1398
1399 numeratorTE = N1 * std::cos(incidentangle) - N2 * cosPhi;
1400 denominatorTE = N1 * std::cos(incidentangle) + N2 * cosPhi;
1401 rTE = numeratorTE / denominatorTE;
1402
1403 numeratorTM = N2 * std::cos(incidentangle) - N1 * cosPhi;
1404 denominatorTM = N2 * std::cos(incidentangle) + N1 * cosPhi;
1405 rTM = numeratorTM / denominatorTM;
1406
1407 // This is my (PG) calculaton for reflectivity on a metallic surface
1408 // depending on the fraction of TE and TM polarization
1409 // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1410 // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1411
1412 reflectivity_TE = (rTE * conj(rTE)) * (E1_perp * E1_perp) /
1413 (E1_perp * E1_perp + E1_parl * E1_parl);
1414 reflectivity_TM = (rTM * conj(rTM)) * (E1_parl * E1_parl) /
1415 (E1_perp * E1_perp + E1_parl * E1_parl);
1416 reflectivity = reflectivity_TE + reflectivity_TM;
1417
1418 do
1419 {
1420 if(G4UniformRand() * real(reflectivity) > real(reflectivity_TE))
1421 {
1422 f_iTE = -1;
1423 }
1424 else
1425 {
1426 f_iTE = 1;
1427 }
1428 if(G4UniformRand() * real(reflectivity) > real(reflectivity_TM))
1429 {
1430 f_iTM = -1;
1431 }
1432 else
1433 {
1434 f_iTM = 1;
1435 }
1436 // Loop checking, 13-Aug-2015, Peter Gumplinger
1437 } while(f_iTE < 0 && f_iTM < 0);
1438
1439 return real(reflectivity);
1440 }
1441
1442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1443 void G4OpBoundaryProcess::CalculateReflectivity()
1444 {
1445 G4double realRindex = fRealRIndexMPV->Value(fPhotonMomentum, idx_rrindex);
1446 G4double imaginaryRindex =
1447 fImagRIndexMPV->Value(fPhotonMomentum, idx_irindex);
1448
1449 // calculate FacetNormal
1450 if(fFinish == ground)
1451 {
1452 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1453 }
1454 else
1455 {
1456 fFacetNormal = fGlobalNormal;
1457 }
1458
1459 G4double cost1 = -fOldMomentum * fFacetNormal;
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(fFacetNormal)).unit();
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 and we set the parallel
1484 // component = 1 in case of a ray perpendicular to the surface
1485 E1_perp = 0.0;
1486 E1_parl = 1.0;
1487 }
1488
1489 G4double incidentangle = GetIncidentAngle();
1490
1491 // calculate the reflectivity depending on incident angle,
1492 // polarization and complex refractive
1493 fReflectivity = GetReflectivity(E1_perp, E1_parl, incidentangle, realRindex,
1494 imaginaryRindex);
1495 }
1496
1497 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1498 G4bool G4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
1499 {
1500 G4Step aStep = *pStep;
1501 aStep.AddTotalEnergyDeposit(fPhotonMomentum);
1502
1503 G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
1504 if(sd != nullptr)
1505 return sd->Hit(&aStep);
1506 else
1507 return false;
1508 }
1509
1510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1511 inline void G4OpBoundaryProcess::SetInvokeSD(G4bool flag)
1512 {
1513 fInvokeSD = flag;
1514 G4OpticalParameters::Instance()->SetBoundaryInvokeSD(fInvokeSD);
1515 }
1516
1517 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1518 void G4OpBoundaryProcess::SetVerboseLevel(G4int verbose)
1519 {
1520 verboseLevel = verbose;
1521 G4OpticalParameters::Instance()->SetBoundaryVerboseLevel(verboseLevel);
1522 }
1523
1524 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1525 void G4OpBoundaryProcess::CoatedDielectricDielectric()
1526 {
1527 G4MaterialPropertyVector* pp = nullptr;
1528
1529 G4MaterialPropertiesTable* MPT = fMaterial2->GetMaterialPropertiesTable();
1530 if((pp = MPT->GetProperty(kRINDEX)))
1531 {
1532 fRindex2 = pp->Value(fPhotonMomentum, idx_rindex2);
1533 }
1534
1535 MPT = fOpticalSurface->GetMaterialPropertiesTable();
1536 if((pp = MPT->GetProperty(kCOATEDRINDEX)))
1537 {
1538 fCoatedRindex = pp->Value(fPhotonMomentum, idx_coatedrindex);
1539 }
1540 if(MPT->ConstPropertyExists(kCOATEDTHICKNESS))
1541 {
1542 fCoatedThickness = MPT->GetConstProperty(kCOATEDTHICKNESS);
1543 }
1544 if(MPT->ConstPropertyExists(kCOATEDFRUSTRATEDTRANSMISSION))
1545 {
1546 fCoatedFrustratedTransmission =
1547 (G4bool)MPT->GetConstProperty(kCOATEDFRUSTRATEDTRANSMISSION);
1548 }
1549
1550 G4double sintTL;
1551 G4double wavelength = h_Planck * c_light / fPhotonMomentum;
1552 G4double PdotN;
1553 G4double E1_perp, E1_parl;
1554 G4double s1, E2_perp, E2_parl, E2_total, transCoeff;
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(fOldMomentum, fGlobalNormal);
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 / fCoatedRindex;
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(fFacetNormal);
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 * cost1 + fRindex2 * cost2);
1632 E2_parl = 2. * s1 * E1_parl / (fRindex2 * cost1 + fRindex1 * cost2);
1633 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1634
1635 transCoeff = 1. - GetReflectivityThroughThinLayer(
1636 sintTL, E1_perp, E1_parl, wavelength, cost1, cost2);
1637 if (!G4BooleanRand(transCoeff))
1638 {
1639 if(verboseLevel > 2)
1640 G4cout << "Reflection from " << fMaterial1->GetName() << " to "
1641 << fMaterial2->GetName() << G4endl;
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. * PdotN) * fFacetNormal;
1656
1657 if (fSint1 > 0.0) { // incident ray oblique
1658
1659 E2_parl = fRindex2 * E2_parl / fRindex1 - E1_parl;
1660 E2_perp = E2_perp - E1_perp;
1661 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
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 + C_perp * A_trans;
1669
1670 }
1671 else
1672 { // incident ray perpendicular
1673 if (fRindex2 > fRindex1)
1674 {
1675 fNewPolarization = -fOldPolarization;
1676 }
1677 else
1678 {
1679 fNewPolarization = fOldPolarization;
1680 }
1681 }
1682
1683 } else { // photon gets transmitted
1684 if (verboseLevel > 2)
1685 G4cout << "Transmission from " << fMaterial1->GetName() << " to "
1686 << fMaterial2->GetName() << G4endl;
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 = CoatedDielectricFrustratedTransmission;
1701 }
1702 else
1703 {
1704 fStatus = CoatedDielectricRefraction;
1705 }
1706
1707 if (fSint1 > 0.0) { // incident ray oblique
1708
1709 alpha = cost1 - cost2 * (fRindex2 / fRindex1);
1710 fNewMomentum = fOldMomentum + alpha * fFacetNormal;
1711 fNewMomentum = fNewMomentum.unit();
1712 A_paral = fNewMomentum.cross(A_trans);
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 + C_perp * A_trans;
1719
1720 }
1721 else
1722 { // incident ray perpendicular
1723 fNewMomentum = fOldMomentum;
1724 fNewPolarization = fOldPolarization;
1725 }
1726 }
1727 }
1728
1729 fOldMomentum = fNewMomentum.unit();
1730 fOldPolarization = fNewPolarization.unit();
1731 if ((fStatus == CoatedDielectricFrustratedTransmission) ||
1732 (fStatus == CoatedDielectricRefraction))
1733 {
1734 done = (fNewMomentum * fGlobalNormal <= 0.0);
1735 }
1736 else
1737 {
1738 done = (fNewMomentum * fGlobalNormal >= -fCarTolerance);
1739 }
1740
1741 } while (!done);
1742 }
1743
1744 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1745 G4double G4OpBoundaryProcess::GetReflectivityThroughThinLayer(G4double sinTL,
1746 G4double E1_perp,
1747 G4double E1_parl,
1748 G4double wavelength, G4double cost1, G4double cost2) {
1749 G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
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) { //Frustrated transmission
1760
1761 if (cost1 > 0.0)
1762 {
1763 gammaTL = std::sqrt(fRindex1 * fRindex1 * fSint1 * fSint1 -
1764 fCoatedRindex * fCoatedRindex);
1765 }
1766 else
1767 {
1768 gammaTL = -std::sqrt(fRindex1 * fRindex1 * fSint1 * fSint1 -
1769 fCoatedRindex * fCoatedRindex);
1770 }
1771
1772 // TE
1773 r1toTL = (fRindex1 * cost1 - i * gammaTL) / (fRindex1 * cost1 + i * gammaTL);
1774 rTLto2 = (i * gammaTL - fRindex2 * cost2) / (i * gammaTL + fRindex2 * cost2);
1775 if (cost1 != 0.0)
1776 {
1777 rTE = (r1toTL + rTLto2 * std::exp(-2 * k0 * fCoatedThickness * gammaTL)) /
1778 (1.0 + r1toTL * rTLto2 * std::exp(-2 * k0 * fCoatedThickness * gammaTL));
1779 }
1780 // TM
1781 r1toTL = (fRindex1 * i * gammaTL - fCoatedRindex * fCoatedRindex * cost1) /
1782 (fRindex1 * i * gammaTL + fCoatedRindex * fCoatedRindex * cost1);
1783 rTLto2 = (fCoatedRindex * fCoatedRindex * cost2 - fRindex2 * i * gammaTL) /
1784 (fCoatedRindex * fCoatedRindex * cost2 + fRindex2 * i * gammaTL);
1785 if (cost1 != 0.0)
1786 {
1787 rTM = (r1toTL + rTLto2 * std::exp(-2 * k0 * fCoatedThickness * gammaTL)) /
1788 (1.0 + r1toTL * rTLto2 * std::exp(-2 * k0 * fCoatedThickness * gammaTL));
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 - fCoatedRindex * costTL) / (fRindex1 * cost1 + fCoatedRindex * costTL);
1810 rTLto2 = (fCoatedRindex * costTL - fRindex2 * cost2) / (fCoatedRindex * costTL + fRindex2 * cost2);
1811 if (cost1 != 0.0)
1812 {
1813 rTE = (r1toTL + rTLto2 * std::exp(2.0 * i * k0 * fCoatedRindex * fCoatedThickness * costTL)) /
1814 (1.0 + r1toTL * rTLto2 * std::exp(2.0 * i * k0 * fCoatedRindex * fCoatedThickness * costTL));
1815 }
1816 // TM
1817 r1toTL = (fRindex1 * costTL - fCoatedRindex * cost1) / (fRindex1 * costTL + fCoatedRindex * cost1);
1818 rTLto2 = (fCoatedRindex * cost2 - fRindex2 * costTL) / (fCoatedRindex * cost2 + fRindex2 * costTL);
1819 if (cost1 != 0.0)
1820 {
1821 rTM = (r1toTL + rTLto2 * std::exp(2.0 * i * k0 * fCoatedRindex * fCoatedThickness * costTL)) /
1822 (1.0 + r1toTL * rTLto2 * std::exp(2.0 * i * k0 * fCoatedRindex * fCoatedThickness * costTL));
1823 }
1824 }
1825
1826 Reflectivity_TE = (rTE * conj(rTE)) * (E1_perp * E1_perp) / (E1_perp * E1_perp + E1_parl * E1_parl);
1827 Reflectivity_TM = (rTM * conj(rTM)) * (E1_parl * E1_parl) / (E1_perp * E1_perp + E1_parl * E1_parl);
1828 Reflectivity = Reflectivity_TE + Reflectivity_TM;
1829
1830 return real(Reflectivity);
1831 }
1832