| 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 type)
93 : G4VDiscreteProcess(processName, type)
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() {}
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 if(aTrack.GetStepLength() <= fCarTolerance)
187 {
188 fStatus = StepTooSmall;
189 if(verboseLevel > 1)
190 BoundaryProcessVerbose();
191
192 G4MaterialPropertyVector* groupvel =
193 fMaterial2->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
194 if(groupvel != nullptr)
195 {
196 aParticleChange.ProposeVelocity(
197 groupvel->Value(fPhotonMomentum, idx_groupvel));
198 }
199 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
200 }
201
202 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
203
204 fPhotonMomentum = aParticle->GetTotalMomentum();
205 fOldMomentum = aParticle->GetMomentumDirection();
206 fOldPolarization = aParticle->GetPolarization();
207
208 if(verboseLevel > 1)
209 {
210 G4cout << " Old Momentum Direction: " << fOldMomentum << G4endl
211 << " Old Polarization: " << fOldPolarization << G4endl;
212 }
213
214 G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
215 G4bool valid;
216
217 // ID of Navigator which limits step
218 G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
219 auto iNav = G4TransportationManager::GetTransportationManager()
220 ->GetActiveNavigatorsIterator();
221 fGlobalNormal = (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint, &valid);
222
223 if(valid)
224 {
225 fGlobalNormal = -fGlobalNormal;
226 }
227 else
228 {
229 G4ExceptionDescription ed;
230 ed << " G4OpBoundaryProcess/PostStepDoIt(): "
231 << " The Navigator reports that it returned an invalid normal" << G4endl;
232 G4Exception(
233 "G4OpBoundaryProcess::PostStepDoIt", "OpBoun01", EventMustBeAborted, ed,
234 "Invalid Surface Normal - Geometry must return valid surface normal");
235 }
236
237 if(fOldMomentum * fGlobalNormal > 0.0)
238 {
239 #ifdef G4OPTICAL_DEBUG
240 G4ExceptionDescription ed;
241 ed << " G4OpBoundaryProcess/PostStepDoIt(): fGlobalNormal points in a "
242 "wrong direction. "
243 << G4endl
244 << " The momentum of the photon arriving at interface (oldMomentum)"
245 << " must exit the volume cross in the step. " << G4endl
246 << " So it MUST have dot < 0 with the normal that Exits the new "
247 "volume (globalNormal)."
248 << G4endl << " >> The dot product of oldMomentum and global Normal is "
249 << fOldMomentum * fGlobalNormal << G4endl
250 << " Old Momentum (during step) = " << fOldMomentum << G4endl
251 << " Global Normal (Exiting New Vol) = " << fGlobalNormal << G4endl
252 << G4endl;
253 G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
254 EventMustBeAborted, // Or JustWarning to see if it happens
255 // repeatedly on one ray
256 ed,
257 "Invalid Surface Normal - Geometry must return valid surface "
258 "normal pointing in the right direction");
259 #else
260 fGlobalNormal = -fGlobalNormal;
261 #endif
262 }
263
264 G4MaterialPropertyVector* rIndexMPV = nullptr;
265 G4MaterialPropertiesTable* MPT = fMaterial1->GetMaterialPropertiesTable();
266 if(MPT != nullptr)
267 {
268 rIndexMPV = MPT->GetProperty(kRINDEX);
269 }
270 if(rIndexMPV != nullptr)
271 {
272 fRindex1 = rIndexMPV->Value(fPhotonMomentum, idx_rindex1);
273 }
274 else
275 {
276 fStatus = NoRINDEX;
277 if(verboseLevel > 1)
278 BoundaryProcessVerbose();
279 aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
280 aParticleChange.ProposeTrackStatus(fStopAndKill);
281 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
282 }
283
284 fReflectivity = 1.;
285 fEfficiency = 0.;
286 fTransmittance = 0.;
287 fSurfaceRoughness = 0.;
288 fModel = glisur;
289 fFinish = polished;
290 G4SurfaceType type = dielectric_dielectric;
291
292 rIndexMPV = nullptr;
293 fOpticalSurface = nullptr;
294
295 G4LogicalSurface* surface =
296 G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
297 if(surface == nullptr)
298 {
299 if(thePostPV->GetMotherLogical() == thePrePV->GetLogicalVolume())
300 {
301 surface = G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
302 if(surface == nullptr)
303 {
304 surface =
305 G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
306 }
307 }
308 else
309 {
310 surface = G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
311 if(surface == nullptr)
312 {
313 surface =
314 G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
315 }
316 }
317 }
318
319 if(surface != nullptr)
320 {
321 fOpticalSurface =
322 dynamic_cast<G4OpticalSurface*>(surface->GetSurfaceProperty());
323 }
324 if(fOpticalSurface != nullptr)
325 {
326 type = fOpticalSurface->GetType();
327 fModel = fOpticalSurface->GetModel();
328 fFinish = fOpticalSurface->GetFinish();
329
330 G4MaterialPropertiesTable* sMPT =
331 fOpticalSurface->GetMaterialPropertiesTable();
332 if(sMPT != nullptr)
333 {
334 if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
335 {
336 rIndexMPV = sMPT->GetProperty(kRINDEX);
337 if(rIndexMPV != nullptr)
338 {
339 fRindex2 = rIndexMPV->Value(fPhotonMomentum, idx_rindex_surface);
340 }
341 else
342 {
343 fStatus = NoRINDEX;
344 if(verboseLevel > 1)
345 BoundaryProcessVerbose();
346 aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
347 aParticleChange.ProposeTrackStatus(fStopAndKill);
348 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
349 }
350 }
351
352 fRealRIndexMPV = sMPT->GetProperty(kREALRINDEX);
353 fImagRIndexMPV = sMPT->GetProperty(kIMAGINARYRINDEX);
354 f_iTE = f_iTM = 1;
355
356 G4MaterialPropertyVector* pp;
357 if((pp = sMPT->GetProperty(kREFLECTIVITY)))
358 {
359 fReflectivity = pp->Value(fPhotonMomentum, idx_reflect);
360 }
361 else if(fRealRIndexMPV && fImagRIndexMPV)
362 {
363 CalculateReflectivity();
364 }
365
366 if((pp = sMPT->GetProperty(kEFFICIENCY)))
367 {
368 fEfficiency = pp->Value(fPhotonMomentum, idx_eff);
369 }
370 if((pp = sMPT->GetProperty(kTRANSMITTANCE)))
371 {
372 fTransmittance = pp->Value(fPhotonMomentum, idx_trans);
373 }
374 if(sMPT->ConstPropertyExists(kSURFACEROUGHNESS))
375 {
376 fSurfaceRoughness = sMPT->GetConstProperty(kSURFACEROUGHNESS);
377 }
378
379 if(fModel == unified)
380 {
381 fProb_sl = (pp = sMPT->GetProperty(kSPECULARLOBECONSTANT))
382 ? pp->Value(fPhotonMomentum, idx_lobe)
383 : 0.;
384 fProb_ss = (pp = sMPT->GetProperty(kSPECULARSPIKECONSTANT))
385 ? pp->Value(fPhotonMomentum, idx_spike)
386 : 0.;
387 fProb_bs = (pp = sMPT->GetProperty(kBACKSCATTERCONSTANT))
388 ? pp->Value(fPhotonMomentum, idx_back)
389 : 0.;
390 }
391 } // end of if(sMPT)
392 else if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
393 {
394 aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
395 aParticleChange.ProposeTrackStatus(fStopAndKill);
396 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
397 }
398 } // end of if(fOpticalSurface)
399
400 // DIELECTRIC-DIELECTRIC
401 if(type == dielectric_dielectric)
402 {
403 if(fFinish == polished || fFinish == ground)
404 {
405 if(fMaterial1 == fMaterial2)
406 {
407 fStatus = SameMaterial;
408 if(verboseLevel > 1)
409 BoundaryProcessVerbose();
410 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
411 }
412 MPT = fMaterial2->GetMaterialPropertiesTable();
413 if(MPT != nullptr)
414 {
415 rIndexMPV = MPT->GetProperty(kRINDEX);
416 }
417 if(rIndexMPV != nullptr)
418 {
419 fRindex2 = rIndexMPV->Value(fPhotonMomentum, idx_rindex2);
420 }
421 else
422 {
423 fStatus = NoRINDEX;
424 if(verboseLevel > 1)
425 BoundaryProcessVerbose();
426 aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
427 aParticleChange.ProposeTrackStatus(fStopAndKill);
428 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
429 }
430 }
431 if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
432 {
433 DielectricDielectric();
434 }
435 else
436 {
437 G4double rand = G4UniformRand();
438 if(rand > fReflectivity + fTransmittance)
439 {
440 DoAbsorption();
441 }
442 else if(rand > fReflectivity)
443 {
444 fStatus = Transmission;
445 fNewMomentum = fOldMomentum;
446 fNewPolarization = fOldPolarization;
447 }
448 else
449 {
450 if(fFinish == polishedfrontpainted)
451 {
452 DoReflection();
453 }
454 else if(fFinish == groundfrontpainted)
455 {
456 fStatus = LambertianReflection;
457 DoReflection();
458 }
459 else
460 {
461 DielectricDielectric();
462 }
463 }
464 }
465 }
466 else if(type == dielectric_metal)
467 {
468 DielectricMetal();
469 }
470 else if(type == dielectric_LUT)
471 {
472 DielectricLUT();
473 }
474 else if(type == dielectric_LUTDAVIS)
475 {
476 DielectricLUTDAVIS();
477 }
478 else if(type == dielectric_dichroic)
479 {
480 DielectricDichroic();
481 }
482 else
483 {
484 G4ExceptionDescription ed;
485 ed << " PostStepDoIt(): Illegal boundary type." << G4endl;
486 G4Exception("G4OpBoundaryProcess", "OpBoun04", JustWarning, ed);
487 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
488 }
489
490 fNewMomentum = fNewMomentum.unit();
491 fNewPolarization = fNewPolarization.unit();
492
493 if(verboseLevel > 1)
494 {
495 G4cout << " New Momentum Direction: " << fNewMomentum << G4endl
496 << " New Polarization: " << fNewPolarization << G4endl;
497 BoundaryProcessVerbose();
498 }
499
500 aParticleChange.ProposeMomentumDirection(fNewMomentum);
501 aParticleChange.ProposePolarization(fNewPolarization);
502
503 if(fStatus == FresnelRefraction || fStatus == Transmission)
504 {
505 G4MaterialPropertyVector* groupvel =
506 fMaterial2->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
507 if(groupvel != nullptr)
508 {
509 aParticleChange.ProposeVelocity(
510 groupvel->Value(fPhotonMomentum, idx_groupvel));
511 }
512 }
513
514 if(fStatus == Detection && fInvokeSD)
515 InvokeSD(pStep);
516 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
517 }
518
519 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
520 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
521 {
522 G4cout << " *** ";
523 if(fStatus == Undefined)
524 G4cout << "Undefined";
525 else if(fStatus == Transmission)
526 G4cout << "Transmission";
527 else if(fStatus == FresnelRefraction)
528 G4cout << "FresnelRefraction";
529 else if(fStatus == FresnelReflection)
530 G4cout << "FresnelReflection";
531 else if(fStatus == TotalInternalReflection)
532 G4cout << "TotalInternalReflection";
533 else if(fStatus == LambertianReflection)
534 G4cout << "LambertianReflection";
535 else if(fStatus == LobeReflection)
536 G4cout << "LobeReflection";
537 else if(fStatus == SpikeReflection)
538 G4cout << "SpikeReflection";
539 else if(fStatus == BackScattering)
540 G4cout << "BackScattering";
541 else if(fStatus == PolishedLumirrorAirReflection)
542 G4cout << "PolishedLumirrorAirReflection";
543 else if(fStatus == PolishedLumirrorGlueReflection)
544 G4cout << "PolishedLumirrorGlueReflection";
545 else if(fStatus == PolishedAirReflection)
546 G4cout << "PolishedAirReflection";
547 else if(fStatus == PolishedTeflonAirReflection)
548 G4cout << "PolishedTeflonAirReflection";
549 else if(fStatus == PolishedTiOAirReflection)
550 G4cout << "PolishedTiOAirReflection";
551 else if(fStatus == PolishedTyvekAirReflection)
552 G4cout << "PolishedTyvekAirReflection";
553 else if(fStatus == PolishedVM2000AirReflection)
554 G4cout << "PolishedVM2000AirReflection";
555 else if(fStatus == PolishedVM2000GlueReflection)
556 G4cout << "PolishedVM2000GlueReflection";
557 else if(fStatus == EtchedLumirrorAirReflection)
558 G4cout << "EtchedLumirrorAirReflection";
559 else if(fStatus == EtchedLumirrorGlueReflection)
560 G4cout << "EtchedLumirrorGlueReflection";
561 else if(fStatus == EtchedAirReflection)
562 G4cout << "EtchedAirReflection";
563 else if(fStatus == EtchedTeflonAirReflection)
564 G4cout << "EtchedTeflonAirReflection";
565 else if(fStatus == EtchedTiOAirReflection)
566 G4cout << "EtchedTiOAirReflection";
567 else if(fStatus == EtchedTyvekAirReflection)
568 G4cout << "EtchedTyvekAirReflection";
569 else if(fStatus == EtchedVM2000AirReflection)
570 G4cout << "EtchedVM2000AirReflection";
571 else if(fStatus == EtchedVM2000GlueReflection)
572 G4cout << "EtchedVM2000GlueReflection";
573 else if(fStatus == GroundLumirrorAirReflection)
574 G4cout << "GroundLumirrorAirReflection";
575 else if(fStatus == GroundLumirrorGlueReflection)
576 G4cout << "GroundLumirrorGlueReflection";
577 else if(fStatus == GroundAirReflection)
578 G4cout << "GroundAirReflection";
579 else if(fStatus == GroundTeflonAirReflection)
580 G4cout << "GroundTeflonAirReflection";
581 else if(fStatus == GroundTiOAirReflection)
582 G4cout << "GroundTiOAirReflection";
583 else if(fStatus == GroundTyvekAirReflection)
584 G4cout << "GroundTyvekAirReflection";
585 else if(fStatus == GroundVM2000AirReflection)
586 G4cout << "GroundVM2000AirReflection";
587 else if(fStatus == GroundVM2000GlueReflection)
588 G4cout << "GroundVM2000GlueReflection";
589 else if(fStatus == Absorption)
590 G4cout << "Absorption";
591 else if(fStatus == Detection)
592 G4cout << "Detection";
593 else if(fStatus == NotAtBoundary)
594 G4cout << "NotAtBoundary";
595 else if(fStatus == SameMaterial)
596 G4cout << "SameMaterial";
597 else if(fStatus == StepTooSmall)
598 G4cout << "StepTooSmall";
599 else if(fStatus == NoRINDEX)
600 G4cout << "NoRINDEX";
601 else if(fStatus == Dichroic)
602 G4cout << "Dichroic Transmission";
603 G4cout << " ***" << G4endl;
604 }
605
606 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
607 G4ThreeVector G4OpBoundaryProcess::GetFacetNormal(
608 const G4ThreeVector& momentum, const G4ThreeVector& normal) const
609 {
610 G4ThreeVector facetNormal;
611 if(fModel == unified || fModel == LUT || fModel == DAVIS)
612 {
613 /* This function codes alpha to a random value taken from the
614 distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
615 for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha) is a
616 gaussian distribution with mean 0 and standard deviation sigma_alpha. */
617
618 G4double sigma_alpha = 0.0;
619 if(fOpticalSurface)
620 sigma_alpha = fOpticalSurface->GetSigmaAlpha();
621 if(sigma_alpha == 0.0)
622 {
623 return normal;
624 }
625
626 G4double f_max = std::min(1.0, 4. * sigma_alpha);
627 G4double alpha, phi, sinAlpha;
628
629 do
630 { // Loop checking, 13-Aug-2015, Peter Gumplinger
631 do
632 { // Loop checking, 13-Aug-2015, Peter Gumplinger
633 alpha = G4RandGauss::shoot(0.0, sigma_alpha);
634 sinAlpha = std::sin(alpha);
635 } while(G4UniformRand() * f_max > sinAlpha || alpha >= halfpi);
636
637 phi = G4UniformRand() * twopi;
638 facetNormal.set(sinAlpha * std::cos(phi), sinAlpha * std::sin(phi),
639 std::cos(alpha));
640 facetNormal.rotateUz(normal);
641 } while(momentum * facetNormal >= 0.0);
642 }
643 else
644 {
645 G4double polish = 1.0;
646 if(fOpticalSurface)
647 polish = fOpticalSurface->GetPolish();
648 if(polish < 1.0)
649 {
650 do
651 { // Loop checking, 13-Aug-2015, Peter Gumplinger
652 G4ThreeVector smear;
653 do
654 { // Loop checking, 13-Aug-2015, Peter Gumplinger
655 smear.setX(2. * G4UniformRand() - 1.);
656 smear.setY(2. * G4UniformRand() - 1.);
657 smear.setZ(2. * G4UniformRand() - 1.);
658 } while(smear.mag2() > 1.0);
659 facetNormal = normal + (1. - polish) * smear;
660 } while(momentum * facetNormal >= 0.0);
661 facetNormal = facetNormal.unit();
662 }
663 else
664 {
665 facetNormal = normal;
666 }
667 }
668 return facetNormal;
669 }
670
671 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
672 void G4OpBoundaryProcess::DielectricMetal()
673 {
674 G4int n = 0;
675 G4double rand;
676 G4ThreeVector A_trans;
677
678 do
679 {
680 ++n;
681 rand = G4UniformRand();
682 if(rand > fReflectivity && n == 1)
683 {
684 if(rand > fReflectivity + fTransmittance)
685 {
686 DoAbsorption();
687 }
688 else
689 {
690 fStatus = Transmission;
691 fNewMomentum = fOldMomentum;
692 fNewPolarization = fOldPolarization;
693 }
694 break;
695 }
696 else
697 {
698 if(fRealRIndexMPV && fImagRIndexMPV)
699 {
700 if(n > 1)
701 {
702 CalculateReflectivity();
703 if(!G4BooleanRand(fReflectivity))
704 {
705 DoAbsorption();
706 break;
707 }
708 }
709 }
710 if(fModel == glisur || fFinish == polished)
711 {
712 DoReflection();
713 }
714 else
715 {
716 if(n == 1)
717 ChooseReflection();
718 if(fStatus == LambertianReflection)
719 {
720 DoReflection();
721 }
722 else if(fStatus == BackScattering)
723 {
724 fNewMomentum = -fOldMomentum;
725 fNewPolarization = -fOldPolarization;
726 }
727 else
728 {
729 if(fStatus == LobeReflection)
730 {
731 if(!fRealRIndexMPV || !fImagRIndexMPV)
732 {
733 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
734 }
735 //else
736 // case of complex rindex needs to be implemented
737 }
738 fNewMomentum =
739 fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
740
741 if(f_iTE > 0 && f_iTM > 0)
742 {
743 fNewPolarization =
744 -fOldPolarization +
745 (2. * fOldPolarization * fFacetNormal * fFacetNormal);
746 }
747 else if(f_iTE > 0)
748 {
749 A_trans = (fSint1 > 0.0) ? fOldMomentum.cross(fFacetNormal).unit()
750 : fOldPolarization;
751 fNewPolarization = -A_trans;
752 }
753 else if(f_iTM > 0)
754 {
755 fNewPolarization =
756 -fNewMomentum.cross(A_trans).unit(); // = -A_paral
757 }
758 }
759 }
760 fOldMomentum = fNewMomentum;
761 fOldPolarization = fNewPolarization;
762 }
763 // Loop checking, 13-Aug-2015, Peter Gumplinger
764 } while(fNewMomentum * fGlobalNormal < 0.0);
765 }
766
767 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
768 void G4OpBoundaryProcess::DielectricLUT()
769 {
770 G4int thetaIndex, phiIndex;
771 G4double angularDistVal, thetaRad, phiRad;
772 G4ThreeVector perpVectorTheta, perpVectorPhi;
773
774 fStatus = G4OpBoundaryProcessStatus(
775 G4int(fFinish) + (G4int(NoRINDEX) - G4int(groundbackpainted)));
776
777 G4int thetaIndexMax = fOpticalSurface->GetThetaIndexMax();
778 G4int phiIndexMax = fOpticalSurface->GetPhiIndexMax();
779
780 G4double rand;
781
782 do
783 {
784 rand = G4UniformRand();
785 if(rand > fReflectivity)
786 {
787 if(rand > fReflectivity + fTransmittance)
788 {
789 DoAbsorption();
790 }
791 else
792 {
793 fStatus = Transmission;
794 fNewMomentum = fOldMomentum;
795 fNewPolarization = fOldPolarization;
796 }
797 break;
798 }
799 else
800 {
801 // Calculate Angle between Normal and Photon Momentum
802 G4double anglePhotonToNormal = fOldMomentum.angle(-fGlobalNormal);
803 // Round to closest integer: LBNL model array has 91 values
804 G4int angleIncident = G4lrint(anglePhotonToNormal / CLHEP::deg);
805
806 // Take random angles THETA and PHI,
807 // and see if below Probability - if not - Redo
808 do
809 {
810 thetaIndex = G4RandFlat::shootInt(thetaIndexMax - 1);
811 phiIndex = G4RandFlat::shootInt(phiIndexMax - 1);
812 // Find probability with the new indeces from LUT
813 angularDistVal = fOpticalSurface->GetAngularDistributionValue(
814 angleIncident, thetaIndex, phiIndex);
815 // Loop checking, 13-Aug-2015, Peter Gumplinger
816 } while(!G4BooleanRand(angularDistVal));
817
818 thetaRad = G4double(-90 + 4 * thetaIndex) * pi / 180.;
819 phiRad = G4double(-90 + 5 * phiIndex) * pi / 180.;
820 // Rotate Photon Momentum in Theta, then in Phi
821 fNewMomentum = -fOldMomentum;
822
823 perpVectorTheta = fNewMomentum.cross(fGlobalNormal);
824 if(perpVectorTheta.mag() < fCarTolerance)
825 {
826 perpVectorTheta = fNewMomentum.orthogonal();
827 }
828 fNewMomentum =
829 fNewMomentum.rotate(anglePhotonToNormal - thetaRad, perpVectorTheta);
830 perpVectorPhi = perpVectorTheta.cross(fNewMomentum);
831 fNewMomentum = fNewMomentum.rotate(-phiRad, perpVectorPhi);
832
833 // Rotate Polarization too:
834 fFacetNormal = (fNewMomentum - fOldMomentum).unit();
835 fNewPolarization = -fOldPolarization +
836 (2. * fOldPolarization * fFacetNormal * fFacetNormal);
837 }
838 // Loop checking, 13-Aug-2015, Peter Gumplinger
839 } while(fNewMomentum * fGlobalNormal <= 0.0);
840 }
841
842 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
843 void G4OpBoundaryProcess::DielectricLUTDAVIS()
844 {
845 G4int angindex, random, angleIncident;
846 G4double reflectivityValue, elevation, azimuth;
847 G4double anglePhotonToNormal;
848
849 G4int lutbin = fOpticalSurface->GetLUTbins();
850 G4double rand = G4UniformRand();
851
852 G4double sinEl;
853 G4ThreeVector u, vNorm, w;
854
855 do
856 {
857 anglePhotonToNormal = fOldMomentum.angle(-fGlobalNormal);
858
859 // Davis model has 90 reflection bins: round down
860 angleIncident = G4lint(anglePhotonToNormal / CLHEP::deg);
861 reflectivityValue = fOpticalSurface->GetReflectivityLUTValue(angleIncident);
862
863 if(rand > reflectivityValue)
864 {
865 if(fEfficiency > 0.)
866 {
867 DoAbsorption();
868 break;
869 }
870 else
871 {
872 fStatus = Transmission;
873
874 if(angleIncident <= 0.01)
875 {
876 fNewMomentum = fOldMomentum;
877 break;
878 }
879
880 do
881 {
882 random = G4RandFlat::shootInt(1, lutbin + 1);
883 angindex =
884 (((random * 2) - 1)) + angleIncident * lutbin * 2 + 3640000;
885
886 azimuth =
887 fOpticalSurface->GetAngularDistributionValueLUT(angindex - 1);
888 elevation = fOpticalSurface->GetAngularDistributionValueLUT(angindex);
889 } while(elevation == 0. && azimuth == 0.);
890
891 sinEl = std::sin(elevation);
892 vNorm = (fGlobalNormal.cross(fOldMomentum)).unit();
893 u = vNorm.cross(fGlobalNormal) * (sinEl * std::cos(azimuth));
894 vNorm *= (sinEl * std::sin(azimuth));
895 // fGlobalNormal shouldn't be modified here
896 w = (fGlobalNormal *= std::cos(elevation));
897 fNewMomentum = u + vNorm + w;
898
899 // Rotate Polarization too:
900 fFacetNormal = (fNewMomentum - fOldMomentum).unit();
901 fNewPolarization = -fOldPolarization + (2. * fOldPolarization *
902 fFacetNormal * fFacetNormal);
903 }
904 }
905 else
906 {
907 fStatus = LobeReflection;
908
909 if(angleIncident == 0)
910 {
911 fNewMomentum = -fOldMomentum;
912 break;
913 }
914
915 do
916 {
917 random = G4RandFlat::shootInt(1, lutbin + 1);
918 angindex = (((random * 2) - 1)) + (angleIncident - 1) * lutbin * 2;
919
920 azimuth = fOpticalSurface->GetAngularDistributionValueLUT(angindex - 1);
921 elevation = fOpticalSurface->GetAngularDistributionValueLUT(angindex);
922 } while(elevation == 0. && azimuth == 0.);
923
924 sinEl = std::sin(elevation);
925 vNorm = (fGlobalNormal.cross(fOldMomentum)).unit();
926 u = vNorm.cross(fGlobalNormal) * (sinEl * std::cos(azimuth));
927 vNorm *= (sinEl * std::sin(azimuth));
928 // fGlobalNormal shouldn't be modified here
929 w = (fGlobalNormal *= std::cos(elevation));
930
931 fNewMomentum = u + vNorm + w;
932
933 // Rotate Polarization too: (needs revision)
934 fNewPolarization = fOldPolarization;
935 }
936 } while(fNewMomentum * fGlobalNormal <= 0.0);
937 }
938
939 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
940 void G4OpBoundaryProcess::DielectricDichroic()
941 {
942 // Calculate Angle between Normal and Photon Momentum
943 G4double anglePhotonToNormal = fOldMomentum.angle(-fGlobalNormal);
944
945 // Round it to closest integer
946 G4double angleIncident = std::floor(180. / pi * anglePhotonToNormal + 0.5);
947
948 if(!fDichroicVector)
949 {
950 if(fOpticalSurface)
951 fDichroicVector = fOpticalSurface->GetDichroicVector();
952 }
953
954 if(fDichroicVector)
955 {
956 G4double wavelength = h_Planck * c_light / fPhotonMomentum;
957 fTransmittance = fDichroicVector->Value(wavelength / nm, angleIncident,
958 idx_dichroicX, idx_dichroicY) *
959 perCent;
960 // G4cout << "wavelength: " << std::floor(wavelength/nm)
961 // << "nm" << G4endl;
962 // G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
963 // G4cout << "Transmittance: "
964 // << std::floor(fTransmittance/perCent) << "%" << G4endl;
965 }
966 else
967 {
968 G4ExceptionDescription ed;
969 ed << " G4OpBoundaryProcess/DielectricDichroic(): "
970 << " The dichroic surface has no G4Physics2DVector" << G4endl;
971 G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
972 FatalException, ed,
973 "A dichroic surface must have an associated G4Physics2DVector");
974 }
975
976 if(!G4BooleanRand(fTransmittance))
977 { // Not transmitted, so reflect
978 if(fModel == glisur || fFinish == polished)
979 {
980 DoReflection();
981 }
982 else
983 {
984 ChooseReflection();
985 if(fStatus == LambertianReflection)
986 {
987 DoReflection();
988 }
989 else if(fStatus == BackScattering)
990 {
991 fNewMomentum = -fOldMomentum;
992 fNewPolarization = -fOldPolarization;
993 }
994 else
995 {
996 G4double PdotN, EdotN;
997 do
998 {
999 if(fStatus == LobeReflection)
1000 {
1001 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1002 }
1003 PdotN = fOldMomentum * fFacetNormal;
1004 fNewMomentum = fOldMomentum - (2. * PdotN) * fFacetNormal;
1005 // Loop checking, 13-Aug-2015, Peter Gumplinger
1006 } while(fNewMomentum * fGlobalNormal <= 0.0);
1007
1008 EdotN = fOldPolarization * fFacetNormal;
1009 fNewPolarization = -fOldPolarization + (2. * EdotN) * fFacetNormal;
1010 }
1011 }
1012 }
1013 else
1014 {
1015 fStatus = Dichroic;
1016 fNewMomentum = fOldMomentum;
1017 fNewPolarization = fOldPolarization;
1018 }
1019 }
1020
1021 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1022 void G4OpBoundaryProcess::DielectricDielectric()
1023 {
1024 G4bool inside = false;
1025 G4bool swap = false;
1026
1027 if(fFinish == polished)
1028 {
1029 fFacetNormal = fGlobalNormal;
1030 }
1031 else
1032 {
1033 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1034 }
1035 G4double cost1 = -fOldMomentum * fFacetNormal;
1036 G4double cost2 = 0.;
1037 G4double sint2 = 0.;
1038
1039 G4bool surfaceRoughnessCriterionPass = true;
1040 if(fSurfaceRoughness != 0. && fRindex1 > fRindex2)
1041 {
1042 G4double wavelength = h_Planck * c_light / fPhotonMomentum;
1043 G4double surfaceRoughnessCriterion = std::exp(-std::pow(
1044 (4. * pi * fSurfaceRoughness * fRindex1 * cost1 / wavelength), 2));
1045 surfaceRoughnessCriterionPass = G4BooleanRand(surfaceRoughnessCriterion);
1046 }
1047
1048 leap:
1049
1050 G4bool through = false;
1051 G4bool done = false;
1052
1053 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1054 G4double E1_perp, E1_parl;
1055 G4double s1, s2, E2_perp, E2_parl, E2_total, transCoeff;
1056 G4double E2_abs, C_parl, C_perp;
1057 G4double alpha;
1058
1059 do
1060 {
1061 if(through)
1062 {
1063 swap = !swap;
1064 through = false;
1065 fGlobalNormal = -fGlobalNormal;
1066 G4SwapPtr(fMaterial1, fMaterial2);
1067 G4SwapObj(&fRindex1, &fRindex2);
1068 }
1069
1070 if(fFinish == polished)
1071 {
1072 fFacetNormal = fGlobalNormal;
1073 }
1074 else
1075 {
1076 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1077 }
1078
1079 cost1 = -fOldMomentum * fFacetNormal;
1080 if(std::abs(cost1) < 1.0 - fCarTolerance)
1081 {
1082 fSint1 = std::sqrt(1. - cost1 * cost1);
1083 sint2 = fSint1 * fRindex1 / fRindex2; // *** Snell's Law ***
1084 // this isn't a sine as we might expect from the name; can be > 1
1085 }
1086 else
1087 {
1088 fSint1 = 0.0;
1089 sint2 = 0.0;
1090 }
1091
1092 // TOTAL INTERNAL REFLECTION
1093 if(sint2 >= 1.0)
1094 {
1095 swap = false;
1096
1097 fStatus = TotalInternalReflection;
1098 if(!surfaceRoughnessCriterionPass)
1099 fStatus = LambertianReflection;
1100 if(fModel == unified && fFinish != polished)
1101 ChooseReflection();
1102 if(fStatus == LambertianReflection)
1103 {
1104 DoReflection();
1105 }
1106 else if(fStatus == BackScattering)
1107 {
1108 fNewMomentum = -fOldMomentum;
1109 fNewPolarization = -fOldPolarization;
1110 }
1111 else
1112 {
1113 fNewMomentum =
1114 fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
1115 fNewPolarization = -fOldPolarization + (2. * fOldPolarization *
1116 fFacetNormal * fFacetNormal);
1117 }
1118 }
1119 // NOT TIR
1120 else if(sint2 < 1.0)
1121 {
1122 // Calculate amplitude for transmission (Q = P x N)
1123 if(cost1 > 0.0)
1124 {
1125 cost2 = std::sqrt(1. - sint2 * sint2);
1126 }
1127 else
1128 {
1129 cost2 = -std::sqrt(1. - sint2 * sint2);
1130 }
1131
1132 if(fSint1 > 0.0)
1133 {
1134 A_trans = (fOldMomentum.cross(fFacetNormal)).unit();
1135 E1_perp = fOldPolarization * A_trans;
1136 E1pp = E1_perp * A_trans;
1137 E1pl = fOldPolarization - E1pp;
1138 E1_parl = E1pl.mag();
1139 }
1140 else
1141 {
1142 A_trans = fOldPolarization;
1143 // Here we Follow Jackson's conventions and set the parallel
1144 // component = 1 in case of a ray perpendicular to the surface
1145 E1_perp = 0.0;
1146 E1_parl = 1.0;
1147 }
1148
1149 s1 = fRindex1 * cost1;
1150 E2_perp = 2. * s1 * E1_perp / (fRindex1 * cost1 + fRindex2 * cost2);
1151 E2_parl = 2. * s1 * E1_parl / (fRindex2 * cost1 + fRindex1 * cost2);
1152 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1153 s2 = fRindex2 * cost2 * E2_total;
1154
1155 if(fTransmittance > 0.)
1156 transCoeff = fTransmittance;
1157 else if(cost1 != 0.0)
1158 transCoeff = s2 / s1;
1159 else
1160 transCoeff = 0.0;
1161
1162 // NOT TIR: REFLECTION
1163 if(!G4BooleanRand(transCoeff))
1164 {
1165 swap = false;
1166 fStatus = FresnelReflection;
1167
1168 if(!surfaceRoughnessCriterionPass)
1169 fStatus = LambertianReflection;
1170 if(fModel == unified && fFinish != polished)
1171 ChooseReflection();
1172 if(fStatus == LambertianReflection)
1173 {
1174 DoReflection();
1175 }
1176 else if(fStatus == BackScattering)
1177 {
1178 fNewMomentum = -fOldMomentum;
1179 fNewPolarization = -fOldPolarization;
1180 }
1181 else
1182 {
1183 fNewMomentum =
1184 fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
1185 if(fSint1 > 0.0)
1186 { // incident ray oblique
1187 E2_parl = fRindex2 * E2_parl / fRindex1 - E1_parl;
1188 E2_perp = E2_perp - E1_perp;
1189 E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1190 A_paral = (fNewMomentum.cross(A_trans)).unit();
1191 E2_abs = std::sqrt(E2_total);
1192 C_parl = E2_parl / E2_abs;
1193 C_perp = E2_perp / E2_abs;
1194
1195 fNewPolarization = C_parl * A_paral + C_perp * A_trans;
1196 }
1197 else
1198 { // incident ray perpendicular
1199 if(fRindex2 > fRindex1)
1200 {
1201 fNewPolarization = -fOldPolarization;
1202 }
1203 else
1204 {
1205 fNewPolarization = fOldPolarization;
1206 }
1207 }
1208 }
1209 }
1210 // NOT TIR: TRANSMISSION
1211 else
1212 {
1213 inside = !inside;
1214 through = true;
1215 fStatus = FresnelRefraction;
1216
1217 if(fSint1 > 0.0)
1218 { // incident ray oblique
1219 alpha = cost1 - cost2 * (fRindex2 / fRindex1);
1220 fNewMomentum = (fOldMomentum + alpha * fFacetNormal).unit();
1221 A_paral = (fNewMomentum.cross(A_trans)).unit();
1222 E2_abs = std::sqrt(E2_total);
1223 C_parl = E2_parl / E2_abs;
1224 C_perp = E2_perp / E2_abs;
1225
1226 fNewPolarization = C_parl * A_paral + C_perp * A_trans;
1227 }
1228 else
1229 { // incident ray perpendicular
1230 fNewMomentum = fOldMomentum;
1231 fNewPolarization = fOldPolarization;
1232 }
1233 }
1234 }
1235
1236 fOldMomentum = fNewMomentum.unit();
1237 fOldPolarization = fNewPolarization.unit();
1238
1239 if(fStatus == FresnelRefraction)
1240 {
1241 done = (fNewMomentum * fGlobalNormal <= 0.0);
1242 }
1243 else
1244 {
1245 done = (fNewMomentum * fGlobalNormal >= -fCarTolerance);
1246 }
1247 // Loop checking, 13-Aug-2015, Peter Gumplinger
1248 } while(!done);
1249
1250 if(inside && !swap)
1251 {
1252 if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
1253 {
1254 G4double rand = G4UniformRand();
1255 if(rand > fReflectivity + fTransmittance)
1256 {
1257 DoAbsorption();
1258 }
1259 else if(rand > fReflectivity)
1260 {
1261 fStatus = Transmission;
1262 fNewMomentum = fOldMomentum;
1263 fNewPolarization = fOldPolarization;
1264 }
1265 else
1266 {
1267 if(fStatus != FresnelRefraction)
1268 {
1269 fGlobalNormal = -fGlobalNormal;
1270 }
1271 else
1272 {
1273 swap = !swap;
1274 G4SwapPtr(fMaterial1, fMaterial2);
1275 G4SwapObj(&fRindex1, &fRindex2);
1276 }
1277 if(fFinish == groundbackpainted)
1278 fStatus = LambertianReflection;
1279
1280 DoReflection();
1281
1282 fGlobalNormal = -fGlobalNormal;
1283 fOldMomentum = fNewMomentum;
1284
1285 goto leap;
1286 }
1287 }
1288 }
1289 }
1290
1291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1292 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track&, G4double,
1293 G4ForceCondition* condition)
1294 {
1295 *condition = Forced;
1296 return DBL_MAX;
1297 }
1298
1299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1300 G4double G4OpBoundaryProcess::GetIncidentAngle()
1301 {
1302 return pi - std::acos(fOldMomentum * fFacetNormal /
1303 (fOldMomentum.mag() * fFacetNormal.mag()));
1304 }
1305
1306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1307 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1308 G4double E1_parl,
1309 G4double incidentangle,
1310 G4double realRindex,
1311 G4double imaginaryRindex)
1312 {
1313 G4complex reflectivity, reflectivity_TE, reflectivity_TM;
1314 G4complex N1(fRindex1, 0.), N2(realRindex, imaginaryRindex);
1315 G4complex cosPhi;
1316
1317 G4complex u(1., 0.); // unit number 1
1318
1319 G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1320 G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1321 G4complex denominatorTE, denominatorTM;
1322 G4complex rTM, rTE;
1323
1324 G4MaterialPropertiesTable* MPT = fMaterial1->GetMaterialPropertiesTable();
1325 G4MaterialPropertyVector* ppR = MPT->GetProperty(kREALRINDEX);
1326 G4MaterialPropertyVector* ppI = MPT->GetProperty(kIMAGINARYRINDEX);
1327 if(ppR && ppI)
1328 {
1329 G4double rRindex = ppR->Value(fPhotonMomentum, idx_rrindex);
1330 G4double iRindex = ppI->Value(fPhotonMomentum, idx_irindex);
1331 N1 = G4complex(rRindex, iRindex);
1332 }
1333
1334 // Following two equations, rTM and rTE, are from: "Introduction To Modern
1335 // Optics" written by Fowles
1336 cosPhi = std::sqrt(u - ((std::sin(incidentangle) * std::sin(incidentangle)) *
1337 (N1 * N1) / (N2 * N2)));
1338
1339 numeratorTE = N1 * std::cos(incidentangle) - N2 * cosPhi;
1340 denominatorTE = N1 * std::cos(incidentangle) + N2 * cosPhi;
1341 rTE = numeratorTE / denominatorTE;
1342
1343 numeratorTM = N2 * std::cos(incidentangle) - N1 * cosPhi;
1344 denominatorTM = N2 * std::cos(incidentangle) + N1 * cosPhi;
1345 rTM = numeratorTM / denominatorTM;
1346
1347 // This is my (PG) calculaton for reflectivity on a metallic surface
1348 // depending on the fraction of TE and TM polarization
1349 // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1350 // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1351
1352 reflectivity_TE = (rTE * conj(rTE)) * (E1_perp * E1_perp) /
1353 (E1_perp * E1_perp + E1_parl * E1_parl);
1354 reflectivity_TM = (rTM * conj(rTM)) * (E1_parl * E1_parl) /
1355 (E1_perp * E1_perp + E1_parl * E1_parl);
1356 reflectivity = reflectivity_TE + reflectivity_TM;
1357
1358 do
1359 {
1360 if(G4UniformRand() * real(reflectivity) > real(reflectivity_TE))
1361 {
1362 f_iTE = -1;
1363 }
1364 else
1365 {
1366 f_iTE = 1;
1367 }
1368 if(G4UniformRand() * real(reflectivity) > real(reflectivity_TM))
1369 {
1370 f_iTM = -1;
1371 }
1372 else
1373 {
1374 f_iTM = 1;
1375 }
1376 // Loop checking, 13-Aug-2015, Peter Gumplinger
1377 } while(f_iTE < 0 && f_iTM < 0);
1378
1379 return real(reflectivity);
1380 }
1381
1382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1383
1384 void G4OpBoundaryProcess::CalculateReflectivity()
1385 {
1386 G4double realRindex = fRealRIndexMPV->Value(fPhotonMomentum, idx_rrindex);
1387 G4double imaginaryRindex =
1388 fImagRIndexMPV->Value(fPhotonMomentum, idx_irindex);
1389
1390 // calculate FacetNormal
1391 if(fFinish == ground)
1392 {
1393 fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1394 }
1395 else
1396 {
1397 fFacetNormal = fGlobalNormal;
1398 }
1399
1400 G4double cost1 = -fOldMomentum * fFacetNormal;
1401 if(std::abs(cost1) < 1.0 - fCarTolerance)
1402 {
1403 fSint1 = std::sqrt(1. - cost1 * cost1);
1404 }
1405 else
1406 {
1407 fSint1 = 0.0;
1408 }
1409
1410 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1411 G4double E1_perp, E1_parl;
1412
1413 if(fSint1 > 0.0)
1414 {
1415 A_trans = (fOldMomentum.cross(fFacetNormal)).unit();
1416 E1_perp = fOldPolarization * A_trans;
1417 E1pp = E1_perp * A_trans;
1418 E1pl = fOldPolarization - E1pp;
1419 E1_parl = E1pl.mag();
1420 }
1421 else
1422 {
1423 A_trans = fOldPolarization;
1424 // Here we Follow Jackson's conventions and we set the parallel
1425 // component = 1 in case of a ray perpendicular to the surface
1426 E1_perp = 0.0;
1427 E1_parl = 1.0;
1428 }
1429
1430 G4double incidentangle = GetIncidentAngle();
1431
1432 // calculate the reflectivity depending on incident angle,
1433 // polarization and complex refractive
1434 fReflectivity = GetReflectivity(E1_perp, E1_parl, incidentangle, realRindex,
1435 imaginaryRindex);
1436 }
1437
1438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1439 G4bool G4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
1440 {
1441 G4Step aStep = *pStep;
1442 aStep.AddTotalEnergyDeposit(fPhotonMomentum);
1443
1444 G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
1445 if(sd != nullptr)
1446 return sd->Hit(&aStep);
1447 else
1448 return false;
1449 }
1450
1451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1452 inline void G4OpBoundaryProcess::SetInvokeSD(G4bool flag)
1453 {
1454 fInvokeSD = flag;
1455 G4OpticalParameters::Instance()->SetBoundaryInvokeSD(fInvokeSD);
1456 }
1457
1458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1459 void G4OpBoundaryProcess::SetVerboseLevel(G4int verbose)
1460 {
1461 verboseLevel = verbose;
1462 G4OpticalParameters::Instance()->SetBoundaryVerboseLevel(verboseLevel);
1463 }
1464