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

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

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