source: trunk/source/processes/optical/src/G4OpBoundaryProcess.cc@ 1196

Last change on this file since 1196 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 38.7 KB
Line 
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//
67// Author: Peter Gumplinger
68// adopted from work by Werner Keil - April 2/96
69// mail: gum@triumf.ca
70//
71////////////////////////////////////////////////////////////////////////
72
73#include "G4ios.hh"
74#include "G4OpProcessSubType.hh"
75
76#include "G4OpBoundaryProcess.hh"
77#include "G4GeometryTolerance.hh"
78
79/////////////////////////
80// Class Implementation
81/////////////////////////
82
83 //////////////
84 // Operators
85 //////////////
86
87// G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
88// {
89// }
90
91 /////////////////
92 // Constructors
93 /////////////////
94
95G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
96 G4ProcessType type)
97 : G4VDiscreteProcess(processName, type)
98{
99 if ( verboseLevel > 0) {
100 G4cout << GetProcessName() << " is created " << G4endl;
101 }
102
103 SetProcessSubType(fOpBoundary);
104
105 theStatus = Undefined;
106 theModel = glisur;
107 theFinish = polished;
108 theReflectivity = 1.;
109 theEfficiency = 0.;
110
111 prob_sl = 0.;
112 prob_ss = 0.;
113 prob_bs = 0.;
114
115 kCarTolerance = G4GeometryTolerance::GetInstance()
116 ->GetSurfaceTolerance();
117
118}
119
120// G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
121// {
122// }
123
124 ////////////////
125 // Destructors
126 ////////////////
127
128G4OpBoundaryProcess::~G4OpBoundaryProcess(){}
129
130 ////////////
131 // Methods
132 ////////////
133
134// PostStepDoIt
135// ------------
136//
137G4VParticleChange*
138G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
139{
140 theStatus = Undefined;
141
142 aParticleChange.Initialize(aTrack);
143
144 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
145 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
146
147 if ( verboseLevel > 0 ) {
148 G4cout << " Photon at Boundary! " << G4endl;
149 G4VPhysicalVolume* thePrePV = pPreStepPoint->GetPhysicalVolume();
150 G4VPhysicalVolume* thePostPV = pPostStepPoint->GetPhysicalVolume();
151 if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
152 if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
153 }
154
155 if (pPostStepPoint->GetStepStatus() != fGeomBoundary){
156 theStatus = NotAtBoundary;
157 if ( verboseLevel > 0) BoundaryProcessVerbose();
158 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
159 }
160 if (aTrack.GetStepLength()<=kCarTolerance/2){
161 theStatus = StepTooSmall;
162 if ( verboseLevel > 0) BoundaryProcessVerbose();
163 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
164 }
165
166 Material1 = pPreStepPoint -> GetMaterial();
167 Material2 = pPostStepPoint -> GetMaterial();
168
169 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
170
171 thePhotonMomentum = aParticle->GetTotalMomentum();
172 OldMomentum = aParticle->GetMomentumDirection();
173 OldPolarization = aParticle->GetPolarization();
174
175 if ( verboseLevel > 0 ) {
176 G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
177 G4cout << " Old Polarization: " << OldPolarization << G4endl;
178 }
179
180 G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
181
182 G4Navigator* theNavigator =
183 G4TransportationManager::GetTransportationManager()->
184 GetNavigatorForTracking();
185
186 G4ThreeVector theLocalPoint = theNavigator->
187 GetGlobalToLocalTransform().
188 TransformPoint(theGlobalPoint);
189
190 G4ThreeVector theLocalNormal; // Normal points back into volume
191
192 G4bool valid;
193 theLocalNormal = theNavigator->GetLocalExitNormal(&valid);
194
195 if (valid) {
196 theLocalNormal = -theLocalNormal;
197 }
198 else {
199 G4cerr << " G4OpBoundaryProcess/PostStepDoIt(): "
200 << " The Navigator reports that it returned an invalid normal"
201 << G4endl;
202 G4Exception("G4OpBoundaryProcess::PostStepDoIt",
203 "Invalid Surface Normal",
204 EventMustBeAborted,
205 "Geometry must return valid surface normal");
206 }
207
208 theGlobalNormal = theNavigator->GetLocalToGlobalTransform().
209 TransformAxis(theLocalNormal);
210
211 if (OldMomentum * theGlobalNormal > 0.0) {
212#ifdef G4DEBUG_OPTICAL
213 G4cerr << " G4OpBoundaryProcess/PostStepDoIt(): "
214 << " theGlobalNormal points the wrong direction "
215 << G4endl;
216#endif
217 theGlobalNormal = -theGlobalNormal;
218 }
219
220 G4MaterialPropertiesTable* aMaterialPropertiesTable;
221 G4MaterialPropertyVector* Rindex;
222
223 aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
224 if (aMaterialPropertiesTable) {
225 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
226 }
227 else {
228 theStatus = NoRINDEX;
229 if ( verboseLevel > 0) BoundaryProcessVerbose();
230 aParticleChange.ProposeTrackStatus(fStopAndKill);
231 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
232 }
233
234 if (Rindex) {
235 Rindex1 = Rindex->GetProperty(thePhotonMomentum);
236 }
237 else {
238 theStatus = NoRINDEX;
239 if ( verboseLevel > 0) BoundaryProcessVerbose();
240 aParticleChange.ProposeTrackStatus(fStopAndKill);
241 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
242 }
243
244 theReflectivity = 1.;
245 theEfficiency = 0.;
246
247 theModel = glisur;
248 theFinish = polished;
249
250 G4SurfaceType type = dielectric_dielectric;
251
252 Rindex = NULL;
253 OpticalSurface = NULL;
254
255 G4LogicalSurface* Surface = NULL;
256
257 Surface = G4LogicalBorderSurface::GetSurface
258 (pPreStepPoint ->GetPhysicalVolume(),
259 pPostStepPoint->GetPhysicalVolume());
260
261 if (Surface == NULL){
262 G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
263 ->GetMotherLogical() ==
264 pPreStepPoint->GetPhysicalVolume()
265 ->GetLogicalVolume());
266 if(enteredDaughter){
267 Surface = G4LogicalSkinSurface::GetSurface
268 (pPostStepPoint->GetPhysicalVolume()->
269 GetLogicalVolume());
270 if(Surface == NULL)
271 Surface = G4LogicalSkinSurface::GetSurface
272 (pPreStepPoint->GetPhysicalVolume()->
273 GetLogicalVolume());
274 }
275 else {
276 Surface = G4LogicalSkinSurface::GetSurface
277 (pPreStepPoint->GetPhysicalVolume()->
278 GetLogicalVolume());
279 if(Surface == NULL)
280 Surface = G4LogicalSkinSurface::GetSurface
281 (pPostStepPoint->GetPhysicalVolume()->
282 GetLogicalVolume());
283 }
284 }
285
286 if (Surface) OpticalSurface =
287 dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
288
289 if (OpticalSurface) {
290
291 type = OpticalSurface->GetType();
292 theModel = OpticalSurface->GetModel();
293 theFinish = OpticalSurface->GetFinish();
294
295 aMaterialPropertiesTable = OpticalSurface->
296 GetMaterialPropertiesTable();
297
298 if (aMaterialPropertiesTable) {
299
300 if (theFinish == polishedbackpainted ||
301 theFinish == groundbackpainted ) {
302 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
303 if (Rindex) {
304 Rindex2 = Rindex->GetProperty(thePhotonMomentum);
305 }
306 else {
307 theStatus = NoRINDEX;
308 if ( verboseLevel > 0) BoundaryProcessVerbose();
309 aParticleChange.ProposeTrackStatus(fStopAndKill);
310 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
311 }
312 }
313
314 PropertyPointer =
315 aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
316 PropertyPointer1 =
317 aMaterialPropertiesTable->GetProperty("REALRINDEX");
318 PropertyPointer2 =
319 aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
320
321 iTE = 1;
322 iTM = 1;
323
324 if (PropertyPointer) {
325
326 theReflectivity =
327 PropertyPointer->GetProperty(thePhotonMomentum);
328
329 } else if (PropertyPointer1 && PropertyPointer2) {
330
331 CalculateReflectivity();
332
333 }
334
335 PropertyPointer =
336 aMaterialPropertiesTable->GetProperty("EFFICIENCY");
337 if (PropertyPointer) {
338 theEfficiency =
339 PropertyPointer->GetProperty(thePhotonMomentum);
340 }
341
342 if ( theModel == unified ) {
343 PropertyPointer =
344 aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
345 if (PropertyPointer) {
346 prob_sl =
347 PropertyPointer->GetProperty(thePhotonMomentum);
348 } else {
349 prob_sl = 0.0;
350 }
351
352 PropertyPointer =
353 aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
354 if (PropertyPointer) {
355 prob_ss =
356 PropertyPointer->GetProperty(thePhotonMomentum);
357 } else {
358 prob_ss = 0.0;
359 }
360
361 PropertyPointer =
362 aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
363 if (PropertyPointer) {
364 prob_bs =
365 PropertyPointer->GetProperty(thePhotonMomentum);
366 } else {
367 prob_bs = 0.0;
368 }
369 }
370 }
371 else if (theFinish == polishedbackpainted ||
372 theFinish == groundbackpainted ) {
373 aParticleChange.ProposeTrackStatus(fStopAndKill);
374 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
375 }
376 }
377
378 if (type == dielectric_dielectric ) {
379 if (theFinish == polished || theFinish == ground ) {
380
381 if (Material1 == Material2){
382 theStatus = SameMaterial;
383 if ( verboseLevel > 0) BoundaryProcessVerbose();
384 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
385 }
386 aMaterialPropertiesTable =
387 Material2->GetMaterialPropertiesTable();
388 if (aMaterialPropertiesTable)
389 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
390 if (Rindex) {
391 Rindex2 = Rindex->GetProperty(thePhotonMomentum);
392 }
393 else {
394 theStatus = NoRINDEX;
395 if ( verboseLevel > 0) BoundaryProcessVerbose();
396 aParticleChange.ProposeTrackStatus(fStopAndKill);
397 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
398 }
399 }
400 }
401
402 if (type == dielectric_metal) {
403
404 DielectricMetal();
405
406 // Uncomment the following lines if you wish to have
407 // Transmission instead of Absorption
408 // if (theStatus == Absorption) {
409 // return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
410 // }
411
412 }
413 else if (type == dielectric_LUT) {
414
415 DielectricLUT();
416
417 }
418 else if (type == dielectric_dielectric) {
419
420 if ( theFinish == polishedfrontpainted ||
421 theFinish == groundfrontpainted ) {
422 if( !G4BooleanRand(theReflectivity) ) {
423 DoAbsorption();
424 }
425 else {
426 if ( theFinish == groundfrontpainted )
427 theStatus = LambertianReflection;
428 DoReflection();
429 }
430 }
431 else {
432 if( !G4BooleanRand(theReflectivity) ) {
433 DoAbsorption();
434 }
435 else {
436 DielectricDielectric();
437 }
438 }
439 }
440 else {
441
442 G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
443 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
444
445 }
446
447 NewMomentum = NewMomentum.unit();
448 NewPolarization = NewPolarization.unit();
449
450 if ( verboseLevel > 0) {
451 G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
452 G4cout << " New Polarization: " << NewPolarization << G4endl;
453 BoundaryProcessVerbose();
454 }
455
456 aParticleChange.ProposeMomentumDirection(NewMomentum);
457 aParticleChange.ProposePolarization(NewPolarization);
458
459 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
460}
461
462void G4OpBoundaryProcess::BoundaryProcessVerbose() const
463{
464 if ( theStatus == Undefined )
465 G4cout << " *** Undefined *** " << G4endl;
466 if ( theStatus == FresnelRefraction )
467 G4cout << " *** FresnelRefraction *** " << G4endl;
468 if ( theStatus == FresnelReflection )
469 G4cout << " *** FresnelReflection *** " << G4endl;
470 if ( theStatus == TotalInternalReflection )
471 G4cout << " *** TotalInternalReflection *** " << G4endl;
472 if ( theStatus == LambertianReflection )
473 G4cout << " *** LambertianReflection *** " << G4endl;
474 if ( theStatus == LobeReflection )
475 G4cout << " *** LobeReflection *** " << G4endl;
476 if ( theStatus == SpikeReflection )
477 G4cout << " *** SpikeReflection *** " << G4endl;
478 if ( theStatus == BackScattering )
479 G4cout << " *** BackScattering *** " << G4endl;
480 if ( theStatus == PolishedLumirrorAirReflection )
481 G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
482 if ( theStatus == PolishedLumirrorGlueReflection )
483 G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
484 if ( theStatus == PolishedAirReflection )
485 G4cout << " *** PolishedAirReflection *** " << G4endl;
486 if ( theStatus == PolishedTeflonAirReflection )
487 G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
488 if ( theStatus == PolishedTiOAirReflection )
489 G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
490 if ( theStatus == PolishedTyvekAirReflection )
491 G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
492 if ( theStatus == PolishedVM2000AirReflection )
493 G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
494 if ( theStatus == PolishedVM2000GlueReflection )
495 G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
496 if ( theStatus == EtchedLumirrorAirReflection )
497 G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
498 if ( theStatus == EtchedLumirrorGlueReflection )
499 G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
500 if ( theStatus == EtchedAirReflection )
501 G4cout << " *** EtchedAirReflection *** " << G4endl;
502 if ( theStatus == EtchedTeflonAirReflection )
503 G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
504 if ( theStatus == EtchedTiOAirReflection )
505 G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
506 if ( theStatus == EtchedTyvekAirReflection )
507 G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
508 if ( theStatus == EtchedVM2000AirReflection )
509 G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
510 if ( theStatus == EtchedVM2000GlueReflection )
511 G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
512 if ( theStatus == GroundLumirrorAirReflection )
513 G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
514 if ( theStatus == GroundLumirrorGlueReflection )
515 G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
516 if ( theStatus == GroundAirReflection )
517 G4cout << " *** GroundAirReflection *** " << G4endl;
518 if ( theStatus == GroundTeflonAirReflection )
519 G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
520 if ( theStatus == GroundTiOAirReflection )
521 G4cout << " *** GroundTiOAirReflection *** " << G4endl;
522 if ( theStatus == GroundTyvekAirReflection )
523 G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
524 if ( theStatus == GroundVM2000AirReflection )
525 G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
526 if ( theStatus == GroundVM2000GlueReflection )
527 G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
528 if ( theStatus == Absorption )
529 G4cout << " *** Absorption *** " << G4endl;
530 if ( theStatus == Detection )
531 G4cout << " *** Detection *** " << G4endl;
532 if ( theStatus == NotAtBoundary )
533 G4cout << " *** NotAtBoundary *** " << G4endl;
534 if ( theStatus == SameMaterial )
535 G4cout << " *** SameMaterial *** " << G4endl;
536 if ( theStatus == StepTooSmall )
537 G4cout << " *** StepTooSmall *** " << G4endl;
538 if ( theStatus == NoRINDEX )
539 G4cout << " *** NoRINDEX *** " << G4endl;
540}
541
542G4ThreeVector
543G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
544 const G4ThreeVector& Normal ) const
545{
546 G4ThreeVector FacetNormal;
547
548 if (theModel == unified || theModel == LUT) {
549
550 /* This function code alpha to a random value taken from the
551 distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
552 for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
553 is a gaussian distribution with mean 0 and standard deviation
554 sigma_alpha. */
555
556 G4double alpha;
557
558 G4double sigma_alpha = 0.0;
559 if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
560
561 G4double f_max = std::min(1.0,4.*sigma_alpha);
562
563 do {
564 do {
565 alpha = G4RandGauss::shoot(0.0,sigma_alpha);
566 } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
567
568 G4double phi = G4UniformRand()*twopi;
569
570 G4double SinAlpha = std::sin(alpha);
571 G4double CosAlpha = std::cos(alpha);
572 G4double SinPhi = std::sin(phi);
573 G4double CosPhi = std::cos(phi);
574
575 G4double unit_x = SinAlpha * CosPhi;
576 G4double unit_y = SinAlpha * SinPhi;
577 G4double unit_z = CosAlpha;
578
579 FacetNormal.setX(unit_x);
580 FacetNormal.setY(unit_y);
581 FacetNormal.setZ(unit_z);
582
583 G4ThreeVector tmpNormal = Normal;
584
585 FacetNormal.rotateUz(tmpNormal);
586 } while (Momentum * FacetNormal >= 0.0);
587 }
588 else {
589
590 G4double polish = 1.0;
591 if (OpticalSurface) polish = OpticalSurface->GetPolish();
592
593 if (polish < 1.0) {
594 do {
595 G4ThreeVector smear;
596 do {
597 smear.setX(2.*G4UniformRand()-1.0);
598 smear.setY(2.*G4UniformRand()-1.0);
599 smear.setZ(2.*G4UniformRand()-1.0);
600 } while (smear.mag()>1.0);
601 smear = (1.-polish) * smear;
602 FacetNormal = Normal + smear;
603 } while (Momentum * FacetNormal >= 0.0);
604 FacetNormal = FacetNormal.unit();
605 }
606 else {
607 FacetNormal = Normal;
608 }
609 }
610 return FacetNormal;
611}
612
613void G4OpBoundaryProcess::DielectricMetal()
614{
615 G4int n = 0;
616
617 do {
618
619 n++;
620
621 if( !G4BooleanRand(theReflectivity) && n == 1 ) {
622
623 // Comment out DoAbsorption and uncomment theStatus = Absorption;
624 // if you wish to have Transmission instead of Absorption
625
626 DoAbsorption();
627 // theStatus = Absorption;
628 break;
629
630 }
631 else {
632
633 if (PropertyPointer1 && PropertyPointer2) {
634 if ( n > 1 ) {
635 CalculateReflectivity();
636 if ( !G4BooleanRand(theReflectivity) ) {
637 DoAbsorption();
638 break;
639 }
640 }
641 }
642
643 if ( theModel == glisur || theFinish == polished ) {
644
645 DoReflection();
646
647 } else {
648
649 if ( n == 1 ) ChooseReflection();
650
651 if ( theStatus == LambertianReflection ) {
652 DoReflection();
653 }
654 else if ( theStatus == BackScattering ) {
655 NewMomentum = -OldMomentum;
656 NewPolarization = -OldPolarization;
657 }
658 else {
659
660 if(theStatus==LobeReflection){
661 if ( PropertyPointer1 && PropertyPointer2 ){
662 } else {
663 theFacetNormal =
664 GetFacetNormal(OldMomentum,theGlobalNormal);
665 }
666 }
667
668 G4double PdotN = OldMomentum * theFacetNormal;
669 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
670 G4double EdotN = OldPolarization * theFacetNormal;
671
672 G4ThreeVector A_trans, A_paral;
673
674 if (sint1 > 0.0 ) {
675 A_trans = OldMomentum.cross(theFacetNormal);
676 A_trans = A_trans.unit();
677 } else {
678 A_trans = OldPolarization;
679 }
680 A_paral = NewMomentum.cross(A_trans);
681 A_paral = A_paral.unit();
682
683 if(iTE>0&&iTM>0) {
684 NewPolarization =
685 -OldPolarization + (2.*EdotN)*theFacetNormal;
686 } else if (iTE>0) {
687 NewPolarization = -A_trans;
688 } else if (iTM>0) {
689 NewPolarization = -A_paral;
690 }
691
692 }
693
694 }
695
696 OldMomentum = NewMomentum;
697 OldPolarization = NewPolarization;
698
699 }
700
701 } while (NewMomentum * theGlobalNormal < 0.0);
702}
703
704void G4OpBoundaryProcess::DielectricLUT()
705{
706 G4int thetaIndex, phiIndex;
707 G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
708 G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
709
710 theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) +
711 (G4int(NoRINDEX)-G4int(groundbackpainted)));
712
713 G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
714 G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
715
716 do {
717 if ( !G4BooleanRand(theReflectivity) ) // Not reflected, so Absorbed
718 DoAbsorption();
719 else {
720 // Calculate Angle between Normal and Photon Momentum
721 G4double anglePhotonToNormal =
722 OldMomentum.angle(-theGlobalNormal);
723 // Round it to closest integer
724 G4int angleIncident = G4int(floor(180/pi*anglePhotonToNormal+0.5));
725
726 // Take random angles THETA and PHI,
727 // and see if below Probability - if not - Redo
728 do {
729 thetaIndex = CLHEP::RandFlat::shootInt(thetaIndexMax-1);
730 phiIndex = CLHEP::RandFlat::shootInt(phiIndexMax-1);
731 // Find probability with the new indeces from LUT
732 AngularDistributionValue = OpticalSurface ->
733 GetAngularDistributionValue(angleIncident,
734 thetaIndex,
735 phiIndex);
736 } while ( !G4BooleanRand(AngularDistributionValue) );
737
738 thetaRad = (-90 + 4*thetaIndex)*pi/180;
739 phiRad = (-90 + 5*phiIndex)*pi/180;
740 // Rotate Photon Momentum in Theta, then in Phi
741 NewMomentum = -OldMomentum;
742 PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
743 if (PerpendicularVectorTheta.mag() > kCarTolerance ) {
744 PerpendicularVectorPhi =
745 PerpendicularVectorTheta.cross(NewMomentum);
746 }
747 else {
748 PerpendicularVectorTheta = NewMomentum.orthogonal();
749 PerpendicularVectorPhi =
750 PerpendicularVectorTheta.cross(NewMomentum);
751 }
752 NewMomentum =
753 NewMomentum.rotate(anglePhotonToNormal-thetaRad,
754 PerpendicularVectorTheta);
755 NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
756 // Rotate Polarization too:
757 theFacetNormal = (NewMomentum - OldMomentum).unit();
758 EdotN = OldPolarization * theFacetNormal;
759 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
760 }
761 } while (NewMomentum * theGlobalNormal <= 0.0);
762}
763
764void G4OpBoundaryProcess::DielectricDielectric()
765{
766 G4bool Inside = false;
767 G4bool Swap = false;
768
769 leap:
770
771 G4bool Through = false;
772 G4bool Done = false;
773
774 do {
775
776 if (Through) {
777 Swap = !Swap;
778 Through = false;
779 theGlobalNormal = -theGlobalNormal;
780 G4SwapPtr(Material1,Material2);
781 G4SwapObj(&Rindex1,&Rindex2);
782 }
783
784 if ( theFinish == ground || theFinish == groundbackpainted ) {
785 theFacetNormal =
786 GetFacetNormal(OldMomentum,theGlobalNormal);
787 }
788 else {
789 theFacetNormal = theGlobalNormal;
790 }
791
792 G4double PdotN = OldMomentum * theFacetNormal;
793 G4double EdotN = OldPolarization * theFacetNormal;
794
795 cost1 = - PdotN;
796 if (std::abs(cost1) < 1.0-kCarTolerance){
797 sint1 = std::sqrt(1.-cost1*cost1);
798 sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
799 }
800 else {
801 sint1 = 0.0;
802 sint2 = 0.0;
803 }
804
805 if (sint2 >= 1.0) {
806
807 // Simulate total internal reflection
808
809 if (Swap) Swap = !Swap;
810
811 theStatus = TotalInternalReflection;
812
813 if ( theModel == unified && theFinish != polished )
814 ChooseReflection();
815
816 if ( theStatus == LambertianReflection ) {
817 DoReflection();
818 }
819 else if ( theStatus == BackScattering ) {
820 NewMomentum = -OldMomentum;
821 NewPolarization = -OldPolarization;
822 }
823 else {
824
825 PdotN = OldMomentum * theFacetNormal;
826 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
827 EdotN = OldPolarization * theFacetNormal;
828 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
829
830 }
831 }
832 else if (sint2 < 1.0) {
833
834 // Calculate amplitude for transmission (Q = P x N)
835
836 if (cost1 > 0.0) {
837 cost2 = std::sqrt(1.-sint2*sint2);
838 }
839 else {
840 cost2 = -std::sqrt(1.-sint2*sint2);
841 }
842
843 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
844 G4double E1_perp, E1_parl;
845
846 if (sint1 > 0.0) {
847 A_trans = OldMomentum.cross(theFacetNormal);
848 A_trans = A_trans.unit();
849 E1_perp = OldPolarization * A_trans;
850 E1pp = E1_perp * A_trans;
851 E1pl = OldPolarization - E1pp;
852 E1_parl = E1pl.mag();
853 }
854 else {
855 A_trans = OldPolarization;
856 // Here we Follow Jackson's conventions and we set the
857 // parallel component = 1 in case of a ray perpendicular
858 // to the surface
859 E1_perp = 0.0;
860 E1_parl = 1.0;
861 }
862
863 G4double s1 = Rindex1*cost1;
864 G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
865 G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
866 G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
867 G4double s2 = Rindex2*cost2*E2_total;
868
869 G4double TransCoeff;
870
871 if (cost1 != 0.0) {
872 TransCoeff = s2/s1;
873 }
874 else {
875 TransCoeff = 0.0;
876 }
877
878 G4double E2_abs, C_parl, C_perp;
879
880 if ( !G4BooleanRand(TransCoeff) ) {
881
882 // Simulate reflection
883
884 if (Swap) Swap = !Swap;
885
886 theStatus = FresnelReflection;
887
888 if ( theModel == unified && theFinish != polished )
889 ChooseReflection();
890
891 if ( theStatus == LambertianReflection ) {
892 DoReflection();
893 }
894 else if ( theStatus == BackScattering ) {
895 NewMomentum = -OldMomentum;
896 NewPolarization = -OldPolarization;
897 }
898 else {
899
900 PdotN = OldMomentum * theFacetNormal;
901 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
902
903 if (sint1 > 0.0) { // incident ray oblique
904
905 E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
906 E2_perp = E2_perp - E1_perp;
907 E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
908 A_paral = NewMomentum.cross(A_trans);
909 A_paral = A_paral.unit();
910 E2_abs = std::sqrt(E2_total);
911 C_parl = E2_parl/E2_abs;
912 C_perp = E2_perp/E2_abs;
913
914 NewPolarization = C_parl*A_paral + C_perp*A_trans;
915
916 }
917
918 else { // incident ray perpendicular
919
920 if (Rindex2 > Rindex1) {
921 NewPolarization = - OldPolarization;
922 }
923 else {
924 NewPolarization = OldPolarization;
925 }
926
927 }
928 }
929 }
930 else { // photon gets transmitted
931
932 // Simulate transmission/refraction
933
934 Inside = !Inside;
935 Through = true;
936 theStatus = FresnelRefraction;
937
938 if (sint1 > 0.0) { // incident ray oblique
939
940 G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
941 NewMomentum = OldMomentum + alpha*theFacetNormal;
942 NewMomentum = NewMomentum.unit();
943 PdotN = -cost2;
944 A_paral = NewMomentum.cross(A_trans);
945 A_paral = A_paral.unit();
946 E2_abs = std::sqrt(E2_total);
947 C_parl = E2_parl/E2_abs;
948 C_perp = E2_perp/E2_abs;
949
950 NewPolarization = C_parl*A_paral + C_perp*A_trans;
951
952 }
953 else { // incident ray perpendicular
954
955 NewMomentum = OldMomentum;
956 NewPolarization = OldPolarization;
957
958 }
959 }
960 }
961
962 OldMomentum = NewMomentum.unit();
963 OldPolarization = NewPolarization.unit();
964
965 if (theStatus == FresnelRefraction) {
966 Done = (NewMomentum * theGlobalNormal <= 0.0);
967 }
968 else {
969 Done = (NewMomentum * theGlobalNormal >= 0.0);
970 }
971
972 } while (!Done);
973
974 if (Inside && !Swap) {
975 if( theFinish == polishedbackpainted ||
976 theFinish == groundbackpainted ) {
977
978 if( !G4BooleanRand(theReflectivity) ) {
979 DoAbsorption();
980 }
981 else {
982 if (theStatus != FresnelRefraction ) {
983 theGlobalNormal = -theGlobalNormal;
984 }
985 else {
986 Swap = !Swap;
987 G4SwapPtr(Material1,Material2);
988 G4SwapObj(&Rindex1,&Rindex2);
989 }
990 if ( theFinish == groundbackpainted )
991 theStatus = LambertianReflection;
992
993 DoReflection();
994
995 theGlobalNormal = -theGlobalNormal;
996 OldMomentum = NewMomentum;
997
998 goto leap;
999 }
1000 }
1001 }
1002}
1003
1004// GetMeanFreePath
1005// ---------------
1006//
1007G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
1008 G4double ,
1009 G4ForceCondition* condition)
1010{
1011 *condition = Forced;
1012
1013 return DBL_MAX;
1014}
1015
1016G4double G4OpBoundaryProcess::GetIncidentAngle()
1017{
1018 G4double PdotN = OldMomentum * theFacetNormal;
1019 G4double magP= OldMomentum.mag();
1020 G4double magN= theFacetNormal.mag();
1021 G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
1022
1023 return incidentangle;
1024}
1025
1026G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1027 G4double E1_parl,
1028 G4double incidentangle,
1029 G4double RealRindex,
1030 G4double ImaginaryRindex)
1031{
1032
1033 G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1034 G4complex N(RealRindex, ImaginaryRindex);
1035 G4complex CosPhi;
1036
1037 G4complex u(1,0); //unit number 1
1038
1039 G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1040 G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1041 G4complex denominatorTE, denominatorTM;
1042 G4complex rTM, rTE;
1043
1044 // Following two equations, rTM and rTE, are from: "Introduction To Modern
1045 // Optics" written by Fowles
1046
1047 CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));
1048
1049 numeratorTE = std::cos(incidentangle) - N*CosPhi;
1050 denominatorTE = std::cos(incidentangle) + N*CosPhi;
1051 rTE = numeratorTE/denominatorTE;
1052
1053 numeratorTM = N*std::cos(incidentangle) - CosPhi;
1054 denominatorTM = N*std::cos(incidentangle) + CosPhi;
1055 rTM = numeratorTM/denominatorTM;
1056
1057 // This is my calculaton for reflectivity on a metalic surface
1058 // depending on the fraction of TE and TM polarization
1059 // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1060 // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1061
1062 Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
1063 / (E1_perp*E1_perp + E1_parl*E1_parl);
1064 Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
1065 / (E1_perp*E1_perp + E1_parl*E1_parl);
1066 Reflectivity = Reflectivity_TE + Reflectivity_TM;
1067
1068 do {
1069 if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1070 {iTE = -1;}else{iTE = 1;}
1071 if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1072 {iTM = -1;}else{iTM = 1;}
1073 } while(iTE<0&&iTM<0);
1074
1075 return real(Reflectivity);
1076
1077}
1078
1079void G4OpBoundaryProcess::CalculateReflectivity()
1080{
1081 G4double RealRindex =
1082 PropertyPointer1->GetProperty(thePhotonMomentum);
1083 G4double ImaginaryRindex =
1084 PropertyPointer2->GetProperty(thePhotonMomentum);
1085
1086 // calculate FacetNormal
1087 if ( theFinish == ground ) {
1088 theFacetNormal =
1089 GetFacetNormal(OldMomentum, theGlobalNormal);
1090 } else {
1091 theFacetNormal = theGlobalNormal;
1092 }
1093
1094 G4double PdotN = OldMomentum * theFacetNormal;
1095 cost1 = -PdotN;
1096
1097 if (std::abs(cost1) < 1.0 - kCarTolerance) {
1098 sint1 = std::sqrt(1. - cost1*cost1);
1099 } else {
1100 sint1 = 0.0;
1101 }
1102
1103 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1104 G4double E1_perp, E1_parl;
1105
1106 if (sint1 > 0.0 ) {
1107 A_trans = OldMomentum.cross(theFacetNormal);
1108 A_trans = A_trans.unit();
1109 E1_perp = OldPolarization * A_trans;
1110 E1pp = E1_perp * A_trans;
1111 E1pl = OldPolarization - E1pp;
1112 E1_parl = E1pl.mag();
1113 }
1114 else {
1115 A_trans = OldPolarization;
1116 // Here we Follow Jackson's conventions and we set the
1117 // parallel component = 1 in case of a ray perpendicular
1118 // to the surface
1119 E1_perp = 0.0;
1120 E1_parl = 1.0;
1121 }
1122
1123 //calculate incident angle
1124 G4double incidentangle = GetIncidentAngle();
1125
1126 //calculate the reflectivity depending on incident angle,
1127 //polarization and complex refractive
1128
1129 theReflectivity =
1130 GetReflectivity(E1_perp, E1_parl, incidentangle,
1131 RealRindex, ImaginaryRindex);
1132}
Note: See TracBrowser for help on using the repository browser.