| 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 if(fTransmittance > 0.)
1210 transCoeff = fTransmittance;
1211 else if(cost1 != 0.0)
1212 transCoeff = s2 / s1;
1213 else
1214 transCoeff = 0.0;
1215
1216 // NOT TIR: REFLECTION
1217 if(!G4BooleanRand(transCoeff))
1218 {
1219 swap = false;
1220 fStatus = FresnelReflection;
1221
1222 if(!surfaceRoughnessCriterionPass)
1223 fStatus = LambertianReflection;
1224 if(fModel == unified && fFinish != polished)
1225 ChooseReflection();
1226 if(fStatus == LambertianReflection)
1227 {
1228 DoReflection();
1229 }
1230 else if(fStatus == BackScattering)
1231 {
1232 fNewMomentum = -fOldMomentum;
1233 fNewPolarization = -fOldPolarization;
1234 }
1235 else
1236 {
1237 fNewMomentum =
1238 fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
1239 if(fSint1 > 0.0)
1240 { // incident ray oblique
1241 E2_parl = fRindex2 * E2_parl / fRindex1 - E1_parl;
1242 E2_perp = E2_perp - E1_perp;
1243 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1244 A_paral = (fNewMomentum.cross(A_trans)).unit();
1245 E2_abs = std::sqrt(E2_total);
1246 C_parl = E2_parl / E2_abs;
1247 C_perp = E2_perp / E2_abs;
1248
1249 fNewPolarization = C_parl * A_paral + C_perp * A_trans;
1250 }
1251 else
1252 { // incident ray perpendicular
1253 if(fRindex2 > fRindex1)
1254 {
1255 fNewPolarization = -fOldPolarization;
1256 }
1257 else
1258 {
1259 fNewPolarization = fOldPolarization;
1260 }
1261 }
1262 }
1263 }
1264 // NOT TIR: TRANSMISSION
1265 else
1266 {
1267 inside = !inside;
1268 through = true;
1269 fStatus = FresnelRefraction;
1270
1271 if(fSint1 > 0.0)
1272 { // incident ray oblique
1273 alpha = cost1 - cost2 * (fRindex2 / fRindex1);
1274 fNewMomentum = (fOldMomentum + alpha * fFacetNormal).unit();
1275 A_paral = (fNewMomentum.cross(A_trans)).unit();
1276 E2_abs = std::sqrt(E2_total);
1277 C_parl = E2_parl / E2_abs;
1278 C_perp = E2_perp / E2_abs;
1279
1280 fNewPolarization = C_parl * A_paral + C_perp * A_trans;
1281 }
1282 else
1283 { // incident ray perpendicular
1284 fNewMomentum = fOldMomentum;
1285 fNewPolarization = fOldPolarization;
1286 }
1287 }
1288 }
1289
1290 fOldMomentum = fNewMomentum.unit();
1291 fOldPolarization = fNewPolarization.unit();
1292
1293 if(fStatus == FresnelRefraction)
1294 {
1295 done = (fNewMomentum * fGlobalNormal <= 0.0);
1296 }
1297 else
1298 {
1299 done = (fNewMomentum * fGlobalNormal >= -fCarTolerance);
1300 }
1301 // Loop checking, 13-Aug-2015, Peter Gumplinger
1302 } while(!done);
1303
1304 if(inside && !swap)
1305 {
1306 if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
1307 {
1308 G4double rand = G4UniformRand();
1309 if(rand > fReflectivity + fTransmittance)
1310 {
1311 DoAbsorption();
1312 }
1313 else if(rand > fReflectivity)
1314 {
1315 fStatus = Transmission;
1316 fNewMomentum = fOldMomentum;
1317 fNewPolarization = fOldPolarization;
1318 }
1319 else
1320 {
1321 if(fStatus != FresnelRefraction)
1322 {
1323 fGlobalNormal = -fGlobalNormal;
1324 }
1325 else
1326 {
1327 swap = !swap;
1328 G4SwapPtr(fMaterial1, fMaterial2);
1329 G4SwapObj(&fRindex1, &fRindex2);
1330 }
1331 if(fFinish == groundbackpainted)
1332 fStatus = LambertianReflection;
1333
1334 DoReflection();
1335
1336 fGlobalNormal = -fGlobalNormal;
1337 fOldMomentum = fNewMomentum;
1338
1339 goto leap;
1340 }
1341 }
1342 }
1343 }
1344
1345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1346 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track&, G4double,
1347 G4ForceCondition* condition)
1348 {
1349 *condition = Forced;
1350 return DBL_MAX;
1351 }
1352
1353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1354 G4double G4OpBoundaryProcess::GetIncidentAngle()
1355 {
1356 return pi - std::acos(fOldMomentum * fFacetNormal /
1357 (fOldMomentum.mag() * fFacetNormal.mag()));
1358 }
1359
1360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1361 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1362 G4double E1_parl,
1363 G4double incidentangle,
1364 G4double realRindex,
1365 G4double imaginaryRindex)
1366 {
1367 G4complex reflectivity, reflectivity_TE, reflectivity_TM;
1368 G4complex N1(fRindex1, 0.), N2(realRindex, imaginaryRindex);
1369 G4complex cosPhi;
1370
1371 G4complex u(1., 0.); // unit number 1
1372
1373 G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1374 G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1375 G4complex denominatorTE, denominatorTM;
1376 G4complex rTM, rTE;
1377
1378 G4MaterialPropertiesTable* MPT = fMaterial1->GetMaterialPropertiesTable();
1379 G4MaterialPropertyVector* ppR = MPT->GetProperty(kREALRINDEX);
1380 G4MaterialPropertyVector* ppI = MPT->GetProperty(kIMAGINARYRINDEX);
1381 if(ppR && ppI)
1382 {
1383 G4double rRindex = ppR->Value(fPhotonMomentum, idx_rrindex);
1384 G4double iRindex = ppI->Value(fPhotonMomentum, idx_irindex);
1385 N1 = G4complex(rRindex, iRindex);
1386 }
1387
1388 // Following two equations, rTM and rTE, are from: "Introduction To Modern
1389 // Optics" written by Fowles
1390 cosPhi = std::sqrt(u - ((std::sin(incidentangle) * std::sin(incidentangle)) *
1391 (N1 * N1) / (N2 * N2)));
1392
1393 numeratorTE = N1 * std::cos(incidentangle) - N2 * cosPhi;
1394 denominatorTE = N1 * std::cos(incidentangle) + N2 * cosPhi;
1395 rTE = numeratorTE / denominatorTE;
1396
1397 numeratorTM = N2 * std::cos(incidentangle) - N1 * cosPhi;
1398 denominatorTM = N2 * std::cos(incidentangle) + N1 * cosPhi;
1399 rTM = numeratorTM / denominatorTM;
1400
1401 // This is my (PG) calculaton for reflectivity on a metallic surface
1402 // depending on the fraction of TE and TM polarization
1403 // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1404 // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1405
1406 reflectivity_TE = (rTE * conj(rTE)) * (E1_perp * E1_perp) /
1407 (E1_perp * E1_perp + E1_parl * E1_parl);
1408 reflectivity_TM = (rTM * conj(rTM)) * (E1_parl * E1_parl) /
1409 (E1_perp * E1_perp + E1_parl * E1_parl);
1410 reflectivity = reflectivity_TE + reflectivity_TM;
1411
1412 do
1413 {
1414 if(G4UniformRand() * real(reflectivity) > real(reflectivity_TE))
1415 {
1416 f_iTE = -1;
1417 }
1418 else
1419 {
1420 f_iTE = 1;
1421 }
1422 if(G4UniformRand() * real(reflectivity) > real(reflectivity_TM))
1423 {
1424 f_iTM = -1;
1425 }
1426 else
1427 {
1428 f_iTM = 1;
1429 }
1430 // Loop checking, 13-Aug-2015, Peter Gumplinger
1431 } while(f_iTE < 0 && f_iTM < 0);
1432
1433 return real(reflectivity);
1434 }
1435
1436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1437 void G4OpBoundaryProcess::CalculateReflectivity()
1438 {
1439 G4double realRindex = fRealRIndexMPV->Value(fPhotonMomentum, idx_rrindex);
1440 G4double imaginaryRindex =
1441 fImagRIndexMPV->Value(fPhotonMomentum, idx_irindex);
1442
1443 // calculate FacetNormal
1444 if(fFinish == ground)
1445 {
1446 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1447 }
1448 else
1449 {
1450 fFacetNormal = fGlobalNormal;
1451 }
1452
1453 G4double cost1 = -fOldMomentum * fFacetNormal;
1454 if(std::abs(cost1) < 1.0 - fCarTolerance)
1455 {
1456 fSint1 = std::sqrt(1. - cost1 * cost1);
1457 }
1458 else
1459 {
1460 fSint1 = 0.0;
1461 }
1462
1463 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1464 G4double E1_perp, E1_parl;
1465
1466 if(fSint1 > 0.0)
1467 {
1468 A_trans = (fOldMomentum.cross(fFacetNormal)).unit();
1469 E1_perp = fOldPolarization * A_trans;
1470 E1pp = E1_perp * A_trans;
1471 E1pl = fOldPolarization - E1pp;
1472 E1_parl = E1pl.mag();
1473 }
1474 else
1475 {
1476 A_trans = fOldPolarization;
1477 // Here we Follow Jackson's conventions and we set the parallel
1478 // component = 1 in case of a ray perpendicular to the surface
1479 E1_perp = 0.0;
1480 E1_parl = 1.0;
1481 }
1482
1483 G4double incidentangle = GetIncidentAngle();
1484
1485 // calculate the reflectivity depending on incident angle,
1486 // polarization and complex refractive
1487 fReflectivity = GetReflectivity(E1_perp, E1_parl, incidentangle, realRindex,
1488 imaginaryRindex);
1489 }
1490
1491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1492 G4bool G4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
1493 {
1494 G4Step aStep = *pStep;
1495 aStep.AddTotalEnergyDeposit(fPhotonMomentum);
1496
1497 G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
1498 if(sd != nullptr)
1499 return sd->Hit(&aStep);
1500 else
1501 return false;
1502 }
1503
1504 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1505 inline void G4OpBoundaryProcess::SetInvokeSD(G4bool flag)
1506 {
1507 fInvokeSD = flag;
1508 G4OpticalParameters::Instance()->SetBoundaryInvokeSD(fInvokeSD);
1509 }
1510
1511 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1512 void G4OpBoundaryProcess::SetVerboseLevel(G4int verbose)
1513 {
1514 verboseLevel = verbose;
1515 G4OpticalParameters::Instance()->SetBoundaryVerboseLevel(verboseLevel);
1516 }
1517
1518 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1519 void G4OpBoundaryProcess::CoatedDielectricDielectric()
1520 {
1521 G4MaterialPropertyVector* pp = nullptr;
1522
1523 G4MaterialPropertiesTable* MPT = fMaterial2->GetMaterialPropertiesTable();
1524 if((pp = MPT->GetProperty(kRINDEX)))
1525 {
1526 fRindex2 = pp->Value(fPhotonMomentum, idx_rindex2);
1527 }
1528
1529 MPT = fOpticalSurface->GetMaterialPropertiesTable();
1530 if((pp = MPT->GetProperty(kCOATEDRINDEX)))
1531 {
1532 fCoatedRindex = pp->Value(fPhotonMomentum, idx_coatedrindex);
1533 }
1534 if(MPT->ConstPropertyExists(kCOATEDTHICKNESS))
1535 {
1536 fCoatedThickness = MPT->GetConstProperty(kCOATEDTHICKNESS);
1537 }
1538 if(MPT->ConstPropertyExists(kCOATEDFRUSTRATEDTRANSMISSION))
1539 {
1540 fCoatedFrustratedTransmission =
1541 (G4bool)MPT->GetConstProperty(kCOATEDFRUSTRATEDTRANSMISSION);
1542 }
1543
1544 G4double sintTL;
1545 G4double wavelength = h_Planck * c_light / fPhotonMomentum;
1546 G4double PdotN;
1547 G4double E1_perp, E1_parl;
1548 G4double s1, E2_perp, E2_parl, E2_total, transCoeff;
1549 G4double E2_abs, C_parl, C_perp;
1550 G4double alpha;
1551 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1552 //G4bool Inside = false;
1553 //G4bool Swap = false;
1554 G4bool through = false;
1555 G4bool done = false;
1556
1557 do {
1558 if (through)
1559 {
1560 //Swap = !Swap;
1561 through = false;
1562 fGlobalNormal = -fGlobalNormal;
1563 G4SwapPtr(fMaterial1, fMaterial2);
1564 G4SwapObj(&fRindex1, &fRindex2);
1565 }
1566
1567 if(fFinish == polished)
1568 {
1569 fFacetNormal = fGlobalNormal;
1570 }
1571 else
1572 {
1573 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1574 }
1575
1576 PdotN = fOldMomentum * fFacetNormal;
1577 G4double cost1 = -PdotN;
1578 G4double sint2, cost2 = 0.;
1579
1580 if (std::abs(cost1) < 1.0 - fCarTolerance)
1581 {
1582 fSint1 = std::sqrt(1. - cost1 * cost1);
1583 sint2 = fSint1 * fRindex1 / fRindex2;
1584 sintTL = fSint1 * fRindex1 / fCoatedRindex;
1585 } else
1586 {
1587 fSint1 = 0.0;
1588 sint2 = 0.0;
1589 sintTL = 0.0;
1590 }
1591
1592 if (fSint1 > 0.0)
1593 {
1594 A_trans = fOldMomentum.cross(fFacetNormal);
1595 A_trans = A_trans.unit();
1596 E1_perp = fOldPolarization * A_trans;
1597 E1pp = E1_perp * A_trans;
1598 E1pl = fOldPolarization - E1pp;
1599 E1_parl = E1pl.mag();
1600 }
1601 else
1602 {
1603 A_trans = fOldPolarization;
1604 E1_perp = 0.0;
1605 E1_parl = 1.0;
1606 }
1607
1608 s1 = fRindex1 * cost1;
1609
1610 if (cost1 > 0.0)
1611 {
1612 cost2 = std::sqrt(1. - sint2 * sint2);
1613 }
1614 else
1615 {
1616 cost2 = -std::sqrt(1. - sint2 * sint2);
1617 }
1618
1619 transCoeff = 0.0;
1620
1621 if (sintTL >= 1.0)
1622 { // --> Angle > Angle Limit
1623 //Swap = false;
1624 }
1625 E2_perp = 2. * s1 * E1_perp / (fRindex1 * cost1 + fRindex2 * cost2);
1626 E2_parl = 2. * s1 * E1_parl / (fRindex2 * cost1 + fRindex1 * cost2);
1627 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1628
1629 transCoeff = 1. - GetReflectivityThroughThinLayer(
1630 sintTL, E1_perp, E1_parl, wavelength, cost1, cost2);
1631 if (!G4BooleanRand(transCoeff))
1632 {
1633 if(verboseLevel > 2)
1634 G4cout << "Reflection from " << fMaterial1->GetName() << " to "
1635 << fMaterial2->GetName() << G4endl;
1636
1637 //Swap = false;
1638
1639 if (sintTL >= 1.0)
1640 {
1641 fStatus = TotalInternalReflection;
1642 }
1643 else
1644 {
1645 fStatus = CoatedDielectricReflection;
1646 }
1647
1648 PdotN = fOldMomentum * fFacetNormal;
1649 fNewMomentum = fOldMomentum - (2. * PdotN) * fFacetNormal;
1650
1651 if (fSint1 > 0.0) { // incident ray oblique
1652
1653 E2_parl = fRindex2 * E2_parl / fRindex1 - E1_parl;
1654 E2_perp = E2_perp - E1_perp;
1655 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1656 A_paral = fNewMomentum.cross(A_trans);
1657 A_paral = A_paral.unit();
1658 E2_abs = std::sqrt(E2_total);
1659 C_parl = E2_parl / E2_abs;
1660 C_perp = E2_perp / E2_abs;
1661
1662 fNewPolarization = C_parl * A_paral + C_perp * A_trans;
1663
1664 }
1665 else
1666 { // incident ray perpendicular
1667 if (fRindex2 > fRindex1)
1668 {
1669 fNewPolarization = -fOldPolarization;
1670 }
1671 else
1672 {
1673 fNewPolarization = fOldPolarization;
1674 }
1675 }
1676
1677 } else { // photon gets transmitted
1678 if (verboseLevel > 2)
1679 G4cout << "Transmission from " << fMaterial1->GetName() << " to "
1680 << fMaterial2->GetName() << G4endl;
1681
1682 //Inside = !Inside;
1683 through = true;
1684
1685 if (fEfficiency > 0.)
1686 {
1687 DoAbsorption();
1688 return;
1689 }
1690 else
1691 {
1692 if (sintTL >= 1.0)
1693 {
1694 fStatus = CoatedDielectricFrustratedTransmission;
1695 }
1696 else
1697 {
1698 fStatus = CoatedDielectricRefraction;
1699 }
1700
1701 if (fSint1 > 0.0) { // incident ray oblique
1702
1703 alpha = cost1 - cost2 * (fRindex2 / fRindex1);
1704 fNewMomentum = fOldMomentum + alpha * fFacetNormal;
1705 fNewMomentum = fNewMomentum.unit();
1706 A_paral = fNewMomentum.cross(A_trans);
1707 A_paral = A_paral.unit();
1708 E2_abs = std::sqrt(E2_total);
1709 C_parl = E2_parl / E2_abs;
1710 C_perp = E2_perp / E2_abs;
1711
1712 fNewPolarization = C_parl * A_paral + C_perp * A_trans;
1713
1714 }
1715 else
1716 { // incident ray perpendicular
1717 fNewMomentum = fOldMomentum;
1718 fNewPolarization = fOldPolarization;
1719 }
1720 }
1721 }
1722
1723 fOldMomentum = fNewMomentum.unit();
1724 fOldPolarization = fNewPolarization.unit();
1725 if ((fStatus == CoatedDielectricFrustratedTransmission) ||
1726 (fStatus == CoatedDielectricRefraction))
1727 {
1728 done = (fNewMomentum * fGlobalNormal <= 0.0);
1729 }
1730 else
1731 {
1732 done = (fNewMomentum * fGlobalNormal >= -fCarTolerance);
1733 }
1734
1735 } while (!done);
1736 }
1737
1738 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1739 G4double G4OpBoundaryProcess::GetReflectivityThroughThinLayer(G4double sinTL,
1740 G4double E1_perp,
1741 G4double E1_parl,
1742 G4double wavelength, G4double cost1, G4double cost2) {
1743 G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1744 G4double gammaTL, costTL;
1745
1746 G4complex i(0, 1);
1747 G4complex rTM, rTE;
1748 G4complex r1toTL, rTLto2;
1749 G4double k0 = 2 * pi / wavelength;
1750
1751 // Angle > Angle limit
1752 if (sinTL >= 1.0) {
1753 if (fCoatedFrustratedTransmission) { //Frustrated transmission
1754
1755 if (cost1 > 0.0)
1756 {
1757 gammaTL = std::sqrt(fRindex1 * fRindex1 * fSint1 * fSint1 -
1758 fCoatedRindex * fCoatedRindex);
1759 }
1760 else
1761 {
1762 gammaTL = -std::sqrt(fRindex1 * fRindex1 * fSint1 * fSint1 -
1763 fCoatedRindex * fCoatedRindex);
1764 }
1765
1766 // TE
1767 r1toTL = (fRindex1 * cost1 - i * gammaTL) / (fRindex1 * cost1 + i * gammaTL);
1768 rTLto2 = (i * gammaTL - fRindex2 * cost2) / (i * gammaTL + fRindex2 * cost2);
1769 if (cost1 != 0.0)
1770 {
1771 rTE = (r1toTL + rTLto2 * std::exp(-2 * k0 * fCoatedThickness * gammaTL)) /
1772 (1.0 + r1toTL * rTLto2 * std::exp(-2 * k0 * fCoatedThickness * gammaTL));
1773 }
1774 // TM
1775 r1toTL = (fRindex1 * i * gammaTL - fCoatedRindex * fCoatedRindex * cost1) /
1776 (fRindex1 * i * gammaTL + fCoatedRindex * fCoatedRindex * cost1);
1777 rTLto2 = (fCoatedRindex * fCoatedRindex * cost2 - fRindex2 * i * gammaTL) /
1778 (fCoatedRindex * fCoatedRindex * cost2 + fRindex2 * i * gammaTL);
1779 if (cost1 != 0.0)
1780 {
1781 rTM = (r1toTL + rTLto2 * std::exp(-2 * k0 * fCoatedThickness * gammaTL)) /
1782 (1.0 + r1toTL * rTLto2 * std::exp(-2 * k0 * fCoatedThickness * gammaTL));
1783 }
1784 }
1785 else
1786 { //Total reflection
1787 return(1.);
1788 }
1789 }
1790
1791 // Angle <= Angle limit
1792 else //if (sinTL < 1.0)
1793 {
1794 if (cost1 > 0.0)
1795 {
1796 costTL = std::sqrt(1. - sinTL * sinTL);
1797 }
1798 else
1799 {
1800 costTL = -std::sqrt(1. - sinTL * sinTL);
1801 }
1802 // TE
1803 r1toTL = (fRindex1 * cost1 - fCoatedRindex * costTL) / (fRindex1 * cost1 + fCoatedRindex * costTL);
1804 rTLto2 = (fCoatedRindex * costTL - fRindex2 * cost2) / (fCoatedRindex * costTL + fRindex2 * cost2);
1805 if (cost1 != 0.0)
1806 {
1807 rTE = (r1toTL + rTLto2 * std::exp(2.0 * i * k0 * fCoatedRindex * fCoatedThickness * costTL)) /
1808 (1.0 + r1toTL * rTLto2 * std::exp(2.0 * i * k0 * fCoatedRindex * fCoatedThickness * costTL));
1809 }
1810 // TM
1811 r1toTL = (fRindex1 * costTL - fCoatedRindex * cost1) / (fRindex1 * costTL + fCoatedRindex * cost1);
1812 rTLto2 = (fCoatedRindex * cost2 - fRindex2 * costTL) / (fCoatedRindex * cost2 + fRindex2 * costTL);
1813 if (cost1 != 0.0)
1814 {
1815 rTM = (r1toTL + rTLto2 * std::exp(2.0 * i * k0 * fCoatedRindex * fCoatedThickness * costTL)) /
1816 (1.0 + r1toTL * rTLto2 * std::exp(2.0 * i * k0 * fCoatedRindex * fCoatedThickness * costTL));
1817 }
1818 }
1819
1820 Reflectivity_TE = (rTE * conj(rTE)) * (E1_perp * E1_perp) / (E1_perp * E1_perp + E1_parl * E1_parl);
1821 Reflectivity_TM = (rTM * conj(rTM)) * (E1_parl * E1_parl) / (E1_perp * E1_perp + E1_parl * E1_parl);
1822 Reflectivity = Reflectivity_TE + Reflectivity_TM;
1823
1824 return real(Reflectivity);
1825 }
1826