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

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago

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