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

Last change on this file since 1337 was 1337, checked in by garnier, 14 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.