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

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

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