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