Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/optical/src/G4OpBoundaryProcess.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/optical/src/G4OpBoundaryProcess.cc (Version 11.3.0) and /processes/optical/src/G4OpBoundaryProcess.cc (Version 11.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //////////////////////////////////////////////     26 ////////////////////////////////////////////////////////////////////////
 27 // Optical Photon Boundary Process Class Imple     27 // Optical Photon Boundary Process Class Implementation
 28 //////////////////////////////////////////////     28 ////////////////////////////////////////////////////////////////////////
 29 //                                                 29 //
 30 // File:        G4OpBoundaryProcess.cc             30 // File:        G4OpBoundaryProcess.cc
 31 // Description: Discrete Process -- reflection     31 // Description: Discrete Process -- reflection/refraction at
 32 //                                  optical in     32 //                                  optical interfaces
 33 // Version:     1.1                                33 // Version:     1.1
 34 // Created:     1997-06-18                         34 // Created:     1997-06-18
 35 // Modified:    1998-05-25 - Correct parallel      35 // Modified:    1998-05-25 - Correct parallel component of polarization
 36 //                           (thanks to: Stefa     36 //                           (thanks to: Stefano Magni + Giovanni Pieri)
 37 //              1998-05-28 - NULL Rindex point     37 //              1998-05-28 - NULL Rindex pointer before reuse
 38 //                           (thanks to: Stefa     38 //                           (thanks to: Stefano Magni)
 39 //              1998-06-11 - delete *sint1 in      39 //              1998-06-11 - delete *sint1 in oblique reflection
 40 //                           (thanks to: Giova     40 //                           (thanks to: Giovanni Pieri)
 41 //              1998-06-19 - move from GetLoca     41 //              1998-06-19 - move from GetLocalExitNormal() to the new
 42 //                           method: GetLocalE     42 //                           method: GetLocalExitNormal(&valid) to get
 43 //                           the surface norma     43 //                           the surface normal in all cases
 44 //              1998-11-07 - NULL OpticalSurfa     44 //              1998-11-07 - NULL OpticalSurface pointer before use
 45 //                           comparison not sh     45 //                           comparison not sharp for: std::abs(cost1) < 1.0
 46 //                           remove sin1, sin2     46 //                           remove sin1, sin2 in lines 556,567
 47 //                           (thanks to Stefan     47 //                           (thanks to Stefano Magni)
 48 //              1999-10-10 - Accommodate chang     48 //              1999-10-10 - Accommodate changes done in DoAbsorption by
 49 //                           changing logic in     49 //                           changing logic in DielectricMetal
 50 //              2001-10-18 - avoid Linux (gcc-     50 //              2001-10-18 - avoid Linux (gcc-2.95.2) warning about variables
 51 //                           might be used uni     51 //                           might be used uninitialized in this function
 52 //                           moved E2_perp, E2     52 //                           moved E2_perp, E2_parl and E2_total out of 'if'
 53 //              2003-11-27 - Modified line 168     53 //              2003-11-27 - Modified line 168-9 to reflect changes made to
 54 //                           G4OpticalSurface      54 //                           G4OpticalSurface class ( by Fan Lei)
 55 //              2004-02-02 - Set theStatus = U     55 //              2004-02-02 - Set theStatus = Undefined at start of DoIt
 56 //              2005-07-28 - add G4ProcessType     56 //              2005-07-28 - add G4ProcessType to constructor
 57 //              2006-11-04 - add capability of     57 //              2006-11-04 - add capability of calculating the reflectivity
 58 //                           off a metal surfa     58 //                           off a metal surface by way of a complex index
 59 //                           of refraction - T     59 //                           of refraction - Thanks to Sehwook Lee and John
 60 //                           Hauptman (Dept. o     60 //                           Hauptman (Dept. of Physics - Iowa State Univ.)
 61 //              2009-11-10 - add capability of     61 //              2009-11-10 - add capability of simulating surface reflections
 62 //                           with Look-Up-Tabl     62 //                           with Look-Up-Tables (LUT) containing measured
 63 //                           optical reflectan     63 //                           optical reflectance for a variety of surface
 64 //                           treatments - Than     64 //                           treatments - Thanks to Martin Janecek and
 65 //                           William Moses (La     65 //                           William Moses (Lawrence Berkeley National Lab.)
 66 //              2013-06-01 - add the capabilit     66 //              2013-06-01 - add the capability of simulating the transmission
 67 //                           of a dichronic fi     67 //                           of a dichronic filter
 68 //              2017-02-24 - add capability of     68 //              2017-02-24 - add capability of simulating surface reflections
 69 //                           with Look-Up-Tabl     69 //                           with Look-Up-Tables (LUT) developed in DAVIS
 70 //                                                 70 //
 71 // Author:      Peter Gumplinger                   71 // Author:      Peter Gumplinger
 72 //    adopted from work by Werner Keil - April     72 //    adopted from work by Werner Keil - April 2/96
 73 //                                                 73 //
 74 //////////////////////////////////////////////     74 ////////////////////////////////////////////////////////////////////////
 75                                                    75 
 76 #include "G4OpBoundaryProcess.hh"                  76 #include "G4OpBoundaryProcess.hh"
 77                                                    77 
 78 #include "G4ios.hh"                                78 #include "G4ios.hh"
 79 #include "G4GeometryTolerance.hh"                  79 #include "G4GeometryTolerance.hh"
 80 #include "G4LogicalBorderSurface.hh"               80 #include "G4LogicalBorderSurface.hh"
 81 #include "G4LogicalSkinSurface.hh"                 81 #include "G4LogicalSkinSurface.hh"
 82 #include "G4OpProcessSubType.hh"                   82 #include "G4OpProcessSubType.hh"
 83 #include "G4OpticalParameters.hh"                  83 #include "G4OpticalParameters.hh"
 84 #include "G4ParallelWorldProcess.hh"               84 #include "G4ParallelWorldProcess.hh"
 85 #include "G4PhysicalConstants.hh"                  85 #include "G4PhysicalConstants.hh"
 86 #include "G4SystemOfUnits.hh"                      86 #include "G4SystemOfUnits.hh"
 87 #include "G4TransportationManager.hh"              87 #include "G4TransportationManager.hh"
 88 #include "G4VSensitiveDetector.hh"                 88 #include "G4VSensitiveDetector.hh"
 89                                                    89 
 90 //....oooOO0OOooo........oooOO0OOooo........oo     90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 91 G4OpBoundaryProcess::G4OpBoundaryProcess(const     91 G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
 92                                          G4Pro <<  92                                          G4ProcessType type)
 93   : G4VDiscreteProcess(processName, ptype)     <<  93   : G4VDiscreteProcess(processName, type)
 94 {                                                  94 {
 95   Initialise();                                    95   Initialise();
 96                                                    96 
 97   if(verboseLevel > 0)                             97   if(verboseLevel > 0)
 98   {                                                98   {
 99     G4cout << GetProcessName() << " is created     99     G4cout << GetProcessName() << " is created " << G4endl;
100   }                                               100   }
101   SetProcessSubType(fOpBoundary);                 101   SetProcessSubType(fOpBoundary);
102                                                   102 
103   fStatus           = Undefined;                  103   fStatus           = Undefined;
104   fModel            = glisur;                     104   fModel            = glisur;
105   fFinish           = polished;                   105   fFinish           = polished;
106   fReflectivity     = 1.;                         106   fReflectivity     = 1.;
107   fEfficiency       = 0.;                         107   fEfficiency       = 0.;
108   fTransmittance    = 0.;                         108   fTransmittance    = 0.;
109   fSurfaceRoughness = 0.;                         109   fSurfaceRoughness = 0.;
110   fProb_sl          = 0.;                         110   fProb_sl          = 0.;
111   fProb_ss          = 0.;                         111   fProb_ss          = 0.;
112   fProb_bs          = 0.;                         112   fProb_bs          = 0.;
113                                                   113 
114   fRealRIndexMPV  = nullptr;                      114   fRealRIndexMPV  = nullptr;
115   fImagRIndexMPV  = nullptr;                      115   fImagRIndexMPV  = nullptr;
116   fMaterial1      = nullptr;                      116   fMaterial1      = nullptr;
117   fMaterial2      = nullptr;                      117   fMaterial2      = nullptr;
118   fOpticalSurface = nullptr;                      118   fOpticalSurface = nullptr;
119   fCarTolerance   = G4GeometryTolerance::GetIn    119   fCarTolerance   = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
120                                                   120 
121   f_iTE = f_iTM   = 0;                            121   f_iTE = f_iTM   = 0;
122   fPhotonMomentum = 0.;                           122   fPhotonMomentum = 0.;
123   fRindex1 = fRindex2 = 1.;                       123   fRindex1 = fRindex2 = 1.;
124   fSint1              = 0.;                       124   fSint1              = 0.;
125   fDichroicVector     = nullptr;                  125   fDichroicVector     = nullptr;
126 }                                                 126 }
127                                                   127 
128 //....oooOO0OOooo........oooOO0OOooo........oo    128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 G4OpBoundaryProcess::~G4OpBoundaryProcess() =  << 129 G4OpBoundaryProcess::~G4OpBoundaryProcess() {}
130                                                   130 
131 //....oooOO0OOooo........oooOO0OOooo........oo    131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 void G4OpBoundaryProcess::PreparePhysicsTable(    132 void G4OpBoundaryProcess::PreparePhysicsTable(const G4ParticleDefinition&)
133 {                                                 133 {
134   Initialise();                                   134   Initialise();
135 }                                                 135 }
136                                                   136 
137 //....oooOO0OOooo........oooOO0OOooo........oo    137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
138 void G4OpBoundaryProcess::Initialise()            138 void G4OpBoundaryProcess::Initialise()
139 {                                                 139 {
140   G4OpticalParameters* params = G4OpticalParam    140   G4OpticalParameters* params = G4OpticalParameters::Instance();
141   SetInvokeSD(params->GetBoundaryInvokeSD());     141   SetInvokeSD(params->GetBoundaryInvokeSD());
142   SetVerboseLevel(params->GetBoundaryVerboseLe    142   SetVerboseLevel(params->GetBoundaryVerboseLevel());
143 }                                                 143 }
144                                                   144 
145 //....oooOO0OOooo........oooOO0OOooo........oo    145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146 G4VParticleChange* G4OpBoundaryProcess::PostSt    146 G4VParticleChange* G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack,
147                                                   147                                                      const G4Step& aStep)
148 {                                                 148 {
149   fStatus = Undefined;                            149   fStatus = Undefined;
150   aParticleChange.Initialize(aTrack);             150   aParticleChange.Initialize(aTrack);
151   aParticleChange.ProposeVelocity(aTrack.GetVe    151   aParticleChange.ProposeVelocity(aTrack.GetVelocity());
152                                                   152 
153   // Get hyperStep from  G4ParallelWorldProces    153   // Get hyperStep from  G4ParallelWorldProcess
154   //  NOTE: PostSetpDoIt of this process to be    154   //  NOTE: PostSetpDoIt of this process to be invoked after
155   //  G4ParallelWorldProcess!                     155   //  G4ParallelWorldProcess!
156   const G4Step* pStep = &aStep;                   156   const G4Step* pStep = &aStep;
157   const G4Step* hStep = G4ParallelWorldProcess    157   const G4Step* hStep = G4ParallelWorldProcess::GetHyperStep();
158   if(hStep != nullptr)                            158   if(hStep != nullptr)
159     pStep = hStep;                                159     pStep = hStep;
160                                                   160 
161   if(pStep->GetPostStepPoint()->GetStepStatus(    161   if(pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
162   {                                               162   {
163     fMaterial1 = pStep->GetPreStepPoint()->Get    163     fMaterial1 = pStep->GetPreStepPoint()->GetMaterial();
164     fMaterial2 = pStep->GetPostStepPoint()->Ge    164     fMaterial2 = pStep->GetPostStepPoint()->GetMaterial();
165   }                                               165   }
166   else                                            166   else
167   {                                               167   {
168     fStatus = NotAtBoundary;                      168     fStatus = NotAtBoundary;
169     if(verboseLevel > 1)                          169     if(verboseLevel > 1)
170       BoundaryProcessVerbose();                   170       BoundaryProcessVerbose();
171     return G4VDiscreteProcess::PostStepDoIt(aT    171     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
172   }                                               172   }
173                                                   173 
174   G4VPhysicalVolume* thePrePV  = pStep->GetPre    174   G4VPhysicalVolume* thePrePV  = pStep->GetPreStepPoint()->GetPhysicalVolume();
175   G4VPhysicalVolume* thePostPV = pStep->GetPos    175   G4VPhysicalVolume* thePostPV = pStep->GetPostStepPoint()->GetPhysicalVolume();
176                                                   176 
177   if(verboseLevel > 1)                            177   if(verboseLevel > 1)
178   {                                               178   {
179     G4cout << " Photon at Boundary! " << G4end    179     G4cout << " Photon at Boundary! " << G4endl;
180     if(thePrePV != nullptr)                       180     if(thePrePV != nullptr)
181       G4cout << " thePrePV:  " << thePrePV->Ge    181       G4cout << " thePrePV:  " << thePrePV->GetName() << G4endl;
182     if(thePostPV != nullptr)                      182     if(thePostPV != nullptr)
183       G4cout << " thePostPV: " << thePostPV->G    183       G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
184   }                                               184   }
185                                                   185 
186   G4double stepLength = aTrack.GetStepLength() << 186   if(aTrack.GetStepLength() <= fCarTolerance)
187   if(stepLength <= fCarTolerance)              << 
188   {                                               187   {
189     fStatus = StepTooSmall;                       188     fStatus = StepTooSmall;
190     if(verboseLevel > 1)                          189     if(verboseLevel > 1)
191       BoundaryProcessVerbose();                   190       BoundaryProcessVerbose();
192                                                   191 
193     G4MaterialPropertyVector* groupvel = nullp << 192     G4MaterialPropertyVector* groupvel =
194     G4MaterialPropertiesTable* aMPT = fMateria << 193       fMaterial2->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
195     if(aMPT != nullptr)                        << 
196     {                                          << 
197       groupvel = aMPT->GetProperty(kGROUPVEL); << 
198     }                                          << 
199                                                << 
200     if(groupvel != nullptr)                       194     if(groupvel != nullptr)
201     {                                             195     {
202       aParticleChange.ProposeVelocity(            196       aParticleChange.ProposeVelocity(
203         groupvel->Value(fPhotonMomentum, idx_g    197         groupvel->Value(fPhotonMomentum, idx_groupvel));
204     }                                             198     }
205     return G4VDiscreteProcess::PostStepDoIt(aT    199     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
206   }                                               200   }
207   else if (stepLength <= 10.*fCarTolerance &&  << 
208   {  // see bug 2510                           << 
209     ++fNumSmallStepWarnings;                   << 
210     if(verboseLevel > 0)                       << 
211     {                                          << 
212       G4ExceptionDescription ed;               << 
213       ed << "G4OpBoundaryProcess: "            << 
214          << "Opticalphoton step length: " << s << 
215          << "This is larger than the threshold << 
216             "to set status StepTooSmall." << G << 
217          << "Boundary scattering may be incorr << 
218       if(fNumSmallStepWarnings == 10)          << 
219       {                                        << 
220         ed << G4endl << "*** Step size warning << 
221       }                                        << 
222       G4Exception("G4OpBoundaryProcess", "OpBo << 
223     }                                          << 
224   }                                            << 
225                                                   201 
226   const G4DynamicParticle* aParticle = aTrack.    202   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
227                                                   203 
228   fPhotonMomentum  = aParticle->GetTotalMoment    204   fPhotonMomentum  = aParticle->GetTotalMomentum();
229   fOldMomentum     = aParticle->GetMomentumDir    205   fOldMomentum     = aParticle->GetMomentumDirection();
230   fOldPolarization = aParticle->GetPolarizatio    206   fOldPolarization = aParticle->GetPolarization();
231                                                   207 
232   if(verboseLevel > 1)                            208   if(verboseLevel > 1)
233   {                                               209   {
234     G4cout << " Old Momentum Direction: " << f    210     G4cout << " Old Momentum Direction: " << fOldMomentum << G4endl
235            << " Old Polarization:       " << f    211            << " Old Polarization:       " << fOldPolarization << G4endl;
236   }                                               212   }
237                                                   213 
238   G4ThreeVector theGlobalPoint = pStep->GetPos    214   G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
239   G4bool valid;                                   215   G4bool valid;
240                                                   216 
241   // ID of Navigator which limits step            217   // ID of Navigator which limits step
242   G4int hNavId = G4ParallelWorldProcess::GetHy    218   G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
243   auto iNav    = G4TransportationManager::GetT    219   auto iNav    = G4TransportationManager::GetTransportationManager()
244                 ->GetActiveNavigatorsIterator(    220                 ->GetActiveNavigatorsIterator();
245   fGlobalNormal = (iNav[hNavId])->GetGlobalExi    221   fGlobalNormal = (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint, &valid);
246                                                   222 
247   if(valid)                                       223   if(valid)
248   {                                               224   {
249     fGlobalNormal = -fGlobalNormal;               225     fGlobalNormal = -fGlobalNormal;
250   }                                               226   }
251   else                                            227   else
252   {                                               228   {
253     G4ExceptionDescription ed;                    229     G4ExceptionDescription ed;
254     ed << " G4OpBoundaryProcess/PostStepDoIt()    230     ed << " G4OpBoundaryProcess/PostStepDoIt(): "
255        << " The Navigator reports that it retu    231        << " The Navigator reports that it returned an invalid normal" << G4endl;
256     G4Exception(                                  232     G4Exception(
257       "G4OpBoundaryProcess::PostStepDoIt", "Op    233       "G4OpBoundaryProcess::PostStepDoIt", "OpBoun01", EventMustBeAborted, ed,
258       "Invalid Surface Normal - Geometry must     234       "Invalid Surface Normal - Geometry must return valid surface normal");
259   }                                               235   }
260                                                   236 
261   if(fOldMomentum * fGlobalNormal > 0.0)          237   if(fOldMomentum * fGlobalNormal > 0.0)
262   {                                               238   {
263 #ifdef G4OPTICAL_DEBUG                            239 #ifdef G4OPTICAL_DEBUG
264     G4ExceptionDescription ed;                    240     G4ExceptionDescription ed;
265     ed << " G4OpBoundaryProcess/PostStepDoIt()    241     ed << " G4OpBoundaryProcess/PostStepDoIt(): fGlobalNormal points in a "
266           "wrong direction. "                     242           "wrong direction. "
267        << G4endl                                  243        << G4endl
268        << "   The momentum of the photon arriv    244        << "   The momentum of the photon arriving at interface (oldMomentum)"
269        << "   must exit the volume cross in th    245        << "   must exit the volume cross in the step. " << G4endl
270        << "   So it MUST have dot < 0 with the    246        << "   So it MUST have dot < 0 with the normal that Exits the new "
271           "volume (globalNormal)."                247           "volume (globalNormal)."
272        << G4endl << "   >> The dot product of     248        << G4endl << "   >> The dot product of oldMomentum and global Normal is "
273        << fOldMomentum * fGlobalNormal << G4en    249        << fOldMomentum * fGlobalNormal << G4endl
274        << "     Old Momentum  (during step)       250        << "     Old Momentum  (during step)     = " << fOldMomentum << G4endl
275        << "     Global Normal (Exiting New Vol    251        << "     Global Normal (Exiting New Vol) = " << fGlobalNormal << G4endl
276        << G4endl;                                 252        << G4endl;
277     G4Exception("G4OpBoundaryProcess::PostStep    253     G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
278                 EventMustBeAborted,  // Or Jus    254                 EventMustBeAborted,  // Or JustWarning to see if it happens
279                                      // repeat    255                                      // repeatedly on one ray
280                 ed,                               256                 ed,
281                 "Invalid Surface Normal - Geom    257                 "Invalid Surface Normal - Geometry must return valid surface "
282                 "normal pointing in the right     258                 "normal pointing in the right direction");
283 #else                                             259 #else
284     fGlobalNormal = -fGlobalNormal;               260     fGlobalNormal = -fGlobalNormal;
285 #endif                                            261 #endif
286   }                                               262   }
287                                                   263 
288   G4MaterialPropertyVector* rIndexMPV = nullpt    264   G4MaterialPropertyVector* rIndexMPV = nullptr;
289   G4MaterialPropertiesTable* MPT = fMaterial1-    265   G4MaterialPropertiesTable* MPT = fMaterial1->GetMaterialPropertiesTable();
290   if(MPT != nullptr)                              266   if(MPT != nullptr)
291   {                                               267   {
292     rIndexMPV = MPT->GetProperty(kRINDEX);        268     rIndexMPV = MPT->GetProperty(kRINDEX);
293   }                                               269   }
294   if(rIndexMPV != nullptr)                        270   if(rIndexMPV != nullptr)
295   {                                               271   {
296     fRindex1 = rIndexMPV->Value(fPhotonMomentu    272     fRindex1 = rIndexMPV->Value(fPhotonMomentum, idx_rindex1);
297   }                                               273   }
298   else                                            274   else
299   {                                               275   {
300     fStatus = NoRINDEX;                           276     fStatus = NoRINDEX;
301     if(verboseLevel > 1)                          277     if(verboseLevel > 1)
302       BoundaryProcessVerbose();                   278       BoundaryProcessVerbose();
303     aParticleChange.ProposeLocalEnergyDeposit(    279     aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
304     aParticleChange.ProposeTrackStatus(fStopAn    280     aParticleChange.ProposeTrackStatus(fStopAndKill);
305     return G4VDiscreteProcess::PostStepDoIt(aT    281     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
306   }                                               282   }
307                                                   283 
308   fReflectivity      = 1.;                        284   fReflectivity      = 1.;
309   fEfficiency        = 0.;                        285   fEfficiency        = 0.;
310   fTransmittance     = 0.;                        286   fTransmittance     = 0.;
311   fSurfaceRoughness  = 0.;                        287   fSurfaceRoughness  = 0.;
312   fModel             = glisur;                    288   fModel             = glisur;
313   fFinish            = polished;                  289   fFinish            = polished;
314   G4SurfaceType type = dielectric_dielectric;     290   G4SurfaceType type = dielectric_dielectric;
315                                                   291 
316   rIndexMPV       = nullptr;                      292   rIndexMPV       = nullptr;
317   fOpticalSurface = nullptr;                      293   fOpticalSurface = nullptr;
318                                                   294 
319   G4LogicalSurface* surface =                     295   G4LogicalSurface* surface =
320     G4LogicalBorderSurface::GetSurface(thePreP    296     G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
321   if(surface == nullptr)                          297   if(surface == nullptr)
322   {                                               298   {
323     if(thePostPV->GetMotherLogical() == thePre    299     if(thePostPV->GetMotherLogical() == thePrePV->GetLogicalVolume())
324     {                                             300     {
325       surface = G4LogicalSkinSurface::GetSurfa    301       surface = G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
326       if(surface == nullptr)                      302       if(surface == nullptr)
327       {                                           303       {
328         surface =                                 304         surface =
329           G4LogicalSkinSurface::GetSurface(the    305           G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
330       }                                           306       }
331     }                                             307     }
332     else                                          308     else
333     {                                             309     {
334       surface = G4LogicalSkinSurface::GetSurfa    310       surface = G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
335       if(surface == nullptr)                      311       if(surface == nullptr)
336       {                                           312       {
337         surface =                                 313         surface =
338           G4LogicalSkinSurface::GetSurface(the    314           G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
339       }                                           315       }
340     }                                             316     }
341   }                                               317   }
342                                                   318 
343   if(surface != nullptr)                          319   if(surface != nullptr)
344   {                                               320   {
345     fOpticalSurface =                             321     fOpticalSurface =
346       dynamic_cast<G4OpticalSurface*>(surface-    322       dynamic_cast<G4OpticalSurface*>(surface->GetSurfaceProperty());
347   }                                               323   }
348   if(fOpticalSurface != nullptr)                  324   if(fOpticalSurface != nullptr)
349   {                                               325   {
350     type    = fOpticalSurface->GetType();         326     type    = fOpticalSurface->GetType();
351     fModel  = fOpticalSurface->GetModel();        327     fModel  = fOpticalSurface->GetModel();
352     fFinish = fOpticalSurface->GetFinish();       328     fFinish = fOpticalSurface->GetFinish();
353                                                   329 
354     G4MaterialPropertiesTable* sMPT =             330     G4MaterialPropertiesTable* sMPT =
355       fOpticalSurface->GetMaterialPropertiesTa    331       fOpticalSurface->GetMaterialPropertiesTable();
356     if(sMPT != nullptr)                           332     if(sMPT != nullptr)
357     {                                             333     {
358       if(fFinish == polishedbackpainted || fFi    334       if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
359       {                                           335       {
360         rIndexMPV = sMPT->GetProperty(kRINDEX)    336         rIndexMPV = sMPT->GetProperty(kRINDEX);
361         if(rIndexMPV != nullptr)                  337         if(rIndexMPV != nullptr)
362         {                                         338         {
363           fRindex2 = rIndexMPV->Value(fPhotonM    339           fRindex2 = rIndexMPV->Value(fPhotonMomentum, idx_rindex_surface);
364         }                                         340         }
365         else                                      341         else
366         {                                         342         {
367           fStatus = NoRINDEX;                     343           fStatus = NoRINDEX;
368           if(verboseLevel > 1)                    344           if(verboseLevel > 1)
369             BoundaryProcessVerbose();             345             BoundaryProcessVerbose();
370           aParticleChange.ProposeLocalEnergyDe    346           aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
371           aParticleChange.ProposeTrackStatus(f    347           aParticleChange.ProposeTrackStatus(fStopAndKill);
372           return G4VDiscreteProcess::PostStepD    348           return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
373         }                                         349         }
374       }                                           350       }
375                                                   351 
376       fRealRIndexMPV = sMPT->GetProperty(kREAL    352       fRealRIndexMPV = sMPT->GetProperty(kREALRINDEX);
377       fImagRIndexMPV = sMPT->GetProperty(kIMAG    353       fImagRIndexMPV = sMPT->GetProperty(kIMAGINARYRINDEX);
378       f_iTE = f_iTM = 1;                          354       f_iTE = f_iTM = 1;
379                                                   355 
380       G4MaterialPropertyVector* pp;               356       G4MaterialPropertyVector* pp;
381       if((pp = sMPT->GetProperty(kREFLECTIVITY    357       if((pp = sMPT->GetProperty(kREFLECTIVITY)))
382       {                                           358       {
383         fReflectivity = pp->Value(fPhotonMomen    359         fReflectivity = pp->Value(fPhotonMomentum, idx_reflect);
384       }                                           360       }
385       else if(fRealRIndexMPV && fImagRIndexMPV    361       else if(fRealRIndexMPV && fImagRIndexMPV)
386       {                                           362       {
387         CalculateReflectivity();                  363         CalculateReflectivity();
388       }                                           364       }
389                                                   365 
390       if((pp = sMPT->GetProperty(kEFFICIENCY))    366       if((pp = sMPT->GetProperty(kEFFICIENCY)))
391       {                                           367       {
392         fEfficiency = pp->Value(fPhotonMomentu    368         fEfficiency = pp->Value(fPhotonMomentum, idx_eff);
393       }                                           369       }
394       if((pp = sMPT->GetProperty(kTRANSMITTANC    370       if((pp = sMPT->GetProperty(kTRANSMITTANCE)))
395       {                                           371       {
396         fTransmittance = pp->Value(fPhotonMome    372         fTransmittance = pp->Value(fPhotonMomentum, idx_trans);
397       }                                           373       }
398       if(sMPT->ConstPropertyExists(kSURFACEROU    374       if(sMPT->ConstPropertyExists(kSURFACEROUGHNESS))
399       {                                           375       {
400         fSurfaceRoughness = sMPT->GetConstProp    376         fSurfaceRoughness = sMPT->GetConstProperty(kSURFACEROUGHNESS);
401       }                                           377       }
402                                                   378 
403       if(fModel == unified)                       379       if(fModel == unified)
404       {                                           380       {
405         fProb_sl = (pp = sMPT->GetProperty(kSP    381         fProb_sl = (pp = sMPT->GetProperty(kSPECULARLOBECONSTANT))
406                      ? pp->Value(fPhotonMoment    382                      ? pp->Value(fPhotonMomentum, idx_lobe)
407                      : 0.;                        383                      : 0.;
408         fProb_ss = (pp = sMPT->GetProperty(kSP    384         fProb_ss = (pp = sMPT->GetProperty(kSPECULARSPIKECONSTANT))
409                      ? pp->Value(fPhotonMoment    385                      ? pp->Value(fPhotonMomentum, idx_spike)
410                      : 0.;                        386                      : 0.;
411         fProb_bs = (pp = sMPT->GetProperty(kBA    387         fProb_bs = (pp = sMPT->GetProperty(kBACKSCATTERCONSTANT))
412                      ? pp->Value(fPhotonMoment    388                      ? pp->Value(fPhotonMomentum, idx_back)
413                      : 0.;                        389                      : 0.;
414       }                                           390       }
415     }  // end of if(sMPT)                         391     }  // end of if(sMPT)
416     else if(fFinish == polishedbackpainted ||     392     else if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
417     {                                             393     {
418       aParticleChange.ProposeLocalEnergyDeposi    394       aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
419       aParticleChange.ProposeTrackStatus(fStop    395       aParticleChange.ProposeTrackStatus(fStopAndKill);
420       return G4VDiscreteProcess::PostStepDoIt(    396       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
421     }                                             397     }
422   }  // end of if(fOpticalSurface)                398   }  // end of if(fOpticalSurface)
423                                                   399 
424   //  DIELECTRIC-DIELECTRIC                       400   //  DIELECTRIC-DIELECTRIC
425   if(type == dielectric_dielectric)               401   if(type == dielectric_dielectric)
426   {                                               402   {
427     if(fFinish == polished || fFinish == groun    403     if(fFinish == polished || fFinish == ground)
428     {                                             404     {
429       if(fMaterial1 == fMaterial2)                405       if(fMaterial1 == fMaterial2)
430       {                                           406       {
431         fStatus = SameMaterial;                   407         fStatus = SameMaterial;
432         if(verboseLevel > 1)                      408         if(verboseLevel > 1)
433           BoundaryProcessVerbose();               409           BoundaryProcessVerbose();
434         return G4VDiscreteProcess::PostStepDoI    410         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
435       }                                           411       }
436       MPT       = fMaterial2->GetMaterialPrope << 412       MPT = fMaterial2->GetMaterialPropertiesTable();
437       rIndexMPV = nullptr;                     << 
438       if(MPT != nullptr)                          413       if(MPT != nullptr)
439       {                                           414       {
440         rIndexMPV = MPT->GetProperty(kRINDEX);    415         rIndexMPV = MPT->GetProperty(kRINDEX);
441       }                                           416       }
442       if(rIndexMPV != nullptr)                    417       if(rIndexMPV != nullptr)
443       {                                           418       {
444         fRindex2 = rIndexMPV->Value(fPhotonMom    419         fRindex2 = rIndexMPV->Value(fPhotonMomentum, idx_rindex2);
445       }                                           420       }
446       else                                        421       else
447       {                                           422       {
448         fStatus = NoRINDEX;                       423         fStatus = NoRINDEX;
449         if(verboseLevel > 1)                      424         if(verboseLevel > 1)
450           BoundaryProcessVerbose();               425           BoundaryProcessVerbose();
451         aParticleChange.ProposeLocalEnergyDepo    426         aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
452         aParticleChange.ProposeTrackStatus(fSt    427         aParticleChange.ProposeTrackStatus(fStopAndKill);
453         return G4VDiscreteProcess::PostStepDoI    428         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
454       }                                           429       }
455     }                                             430     }
456     if(fFinish == polishedbackpainted || fFini    431     if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
457     {                                             432     {
458       DielectricDielectric();                     433       DielectricDielectric();
459     }                                             434     }
460     else                                          435     else
461     {                                             436     {
462       G4double rand = G4UniformRand();            437       G4double rand = G4UniformRand();
463       if(rand > fReflectivity + fTransmittance    438       if(rand > fReflectivity + fTransmittance)
464       {                                           439       {
465         DoAbsorption();                           440         DoAbsorption();
466       }                                           441       }
467       else if(rand > fReflectivity)               442       else if(rand > fReflectivity)
468       {                                           443       {
469         fStatus          = Transmission;          444         fStatus          = Transmission;
470         fNewMomentum     = fOldMomentum;          445         fNewMomentum     = fOldMomentum;
471         fNewPolarization = fOldPolarization;      446         fNewPolarization = fOldPolarization;
472       }                                           447       }
473       else                                        448       else
474       {                                           449       {
475         if(fFinish == polishedfrontpainted)       450         if(fFinish == polishedfrontpainted)
476         {                                         451         {
477           DoReflection();                         452           DoReflection();
478         }                                         453         }
479         else if(fFinish == groundfrontpainted)    454         else if(fFinish == groundfrontpainted)
480         {                                         455         {
481           fStatus = LambertianReflection;         456           fStatus = LambertianReflection;
482           DoReflection();                         457           DoReflection();
483         }                                         458         }
484         else                                      459         else
485         {                                         460         {
486           DielectricDielectric();                 461           DielectricDielectric();
487         }                                         462         }
488       }                                           463       }
489     }                                             464     }
490   }                                               465   }
491   else if(type == dielectric_metal)               466   else if(type == dielectric_metal)
492   {                                               467   {
493     DielectricMetal();                            468     DielectricMetal();
494   }                                               469   }
495   else if(type == dielectric_LUT)                 470   else if(type == dielectric_LUT)
496   {                                               471   {
497     DielectricLUT();                              472     DielectricLUT();
498   }                                               473   }
499   else if(type == dielectric_LUTDAVIS)            474   else if(type == dielectric_LUTDAVIS)
500   {                                               475   {
501     DielectricLUTDAVIS();                         476     DielectricLUTDAVIS();
502   }                                               477   }
503   else if(type == dielectric_dichroic)            478   else if(type == dielectric_dichroic)
504   {                                               479   {
505     DielectricDichroic();                         480     DielectricDichroic();
506   }                                               481   }
507   else if(type == coated)                      << 
508   {                                            << 
509     CoatedDielectricDielectric();              << 
510   }                                            << 
511   else                                            482   else
512   {                                               483   {
513     if(fNumBdryTypeWarnings <= 10)             << 484     G4ExceptionDescription ed;
514     {                                          << 485     ed << " PostStepDoIt(): Illegal boundary type." << G4endl;
515       ++fNumBdryTypeWarnings;                  << 486     G4Exception("G4OpBoundaryProcess", "OpBoun04", JustWarning, ed);
516       if(verboseLevel > 0)                     << 
517       {                                        << 
518         G4ExceptionDescription ed;             << 
519         ed << " PostStepDoIt(): Illegal bounda << 
520         if(fNumBdryTypeWarnings == 10)         << 
521         {                                      << 
522           ed << "** Boundary type warnings sto << 
523         }                                      << 
524         G4Exception("G4OpBoundaryProcess", "Op << 
525       }                                        << 
526     }                                          << 
527     return G4VDiscreteProcess::PostStepDoIt(aT    487     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
528   }                                               488   }
529                                                   489 
530   fNewMomentum     = fNewMomentum.unit();         490   fNewMomentum     = fNewMomentum.unit();
531   fNewPolarization = fNewPolarization.unit();     491   fNewPolarization = fNewPolarization.unit();
532                                                   492 
533   if(verboseLevel > 1)                            493   if(verboseLevel > 1)
534   {                                               494   {
535     G4cout << " New Momentum Direction: " << f    495     G4cout << " New Momentum Direction: " << fNewMomentum << G4endl
536            << " New Polarization:       " << f    496            << " New Polarization:       " << fNewPolarization << G4endl;
537     BoundaryProcessVerbose();                     497     BoundaryProcessVerbose();
538   }                                               498   }
539                                                   499 
540   aParticleChange.ProposeMomentumDirection(fNe    500   aParticleChange.ProposeMomentumDirection(fNewMomentum);
541   aParticleChange.ProposePolarization(fNewPola    501   aParticleChange.ProposePolarization(fNewPolarization);
542                                                   502 
543   if(fStatus == FresnelRefraction || fStatus =    503   if(fStatus == FresnelRefraction || fStatus == Transmission)
544   {                                               504   {
545     // not all surface types check that fMater << 505     G4MaterialPropertyVector* groupvel =
546     G4MaterialPropertiesTable* aMPT = fMateria << 506       fMaterial2->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
547     G4MaterialPropertyVector* groupvel = nullp << 
548     if(aMPT != nullptr)                        << 
549     {                                          << 
550       groupvel = aMPT->GetProperty(kGROUPVEL); << 
551     }                                          << 
552     if(groupvel != nullptr)                       507     if(groupvel != nullptr)
553     {                                             508     {
554       aParticleChange.ProposeVelocity(            509       aParticleChange.ProposeVelocity(
555         groupvel->Value(fPhotonMomentum, idx_g    510         groupvel->Value(fPhotonMomentum, idx_groupvel));
556     }                                             511     }
557   }                                               512   }
558                                                   513 
559   if(fStatus == Detection && fInvokeSD)           514   if(fStatus == Detection && fInvokeSD)
560     InvokeSD(pStep);                              515     InvokeSD(pStep);
561   return G4VDiscreteProcess::PostStepDoIt(aTra    516   return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
562 }                                                 517 }
563                                                   518 
564 //....oooOO0OOooo........oooOO0OOooo........oo    519 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
565 void G4OpBoundaryProcess::BoundaryProcessVerbo    520 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
566 {                                                 521 {
567   G4cout << " *** ";                              522   G4cout << " *** ";
568   if(fStatus == Undefined)                        523   if(fStatus == Undefined)
569     G4cout << "Undefined";                        524     G4cout << "Undefined";
570   else if(fStatus == Transmission)                525   else if(fStatus == Transmission)
571     G4cout << "Transmission";                     526     G4cout << "Transmission";
572   else if(fStatus == FresnelRefraction)           527   else if(fStatus == FresnelRefraction)
573     G4cout << "FresnelRefraction";                528     G4cout << "FresnelRefraction";
574   else if(fStatus == FresnelReflection)           529   else if(fStatus == FresnelReflection)
575     G4cout << "FresnelReflection";                530     G4cout << "FresnelReflection";
576   else if(fStatus == TotalInternalReflection)     531   else if(fStatus == TotalInternalReflection)
577     G4cout << "TotalInternalReflection";          532     G4cout << "TotalInternalReflection";
578   else if(fStatus == LambertianReflection)        533   else if(fStatus == LambertianReflection)
579     G4cout << "LambertianReflection";             534     G4cout << "LambertianReflection";
580   else if(fStatus == LobeReflection)              535   else if(fStatus == LobeReflection)
581     G4cout << "LobeReflection";                   536     G4cout << "LobeReflection";
582   else if(fStatus == SpikeReflection)             537   else if(fStatus == SpikeReflection)
583     G4cout << "SpikeReflection";                  538     G4cout << "SpikeReflection";
584   else if(fStatus == BackScattering)              539   else if(fStatus == BackScattering)
585     G4cout << "BackScattering";                   540     G4cout << "BackScattering";
586   else if(fStatus == PolishedLumirrorAirReflec    541   else if(fStatus == PolishedLumirrorAirReflection)
587     G4cout << "PolishedLumirrorAirReflection";    542     G4cout << "PolishedLumirrorAirReflection";
588   else if(fStatus == PolishedLumirrorGlueRefle    543   else if(fStatus == PolishedLumirrorGlueReflection)
589     G4cout << "PolishedLumirrorGlueReflection"    544     G4cout << "PolishedLumirrorGlueReflection";
590   else if(fStatus == PolishedAirReflection)       545   else if(fStatus == PolishedAirReflection)
591     G4cout << "PolishedAirReflection";            546     G4cout << "PolishedAirReflection";
592   else if(fStatus == PolishedTeflonAirReflecti    547   else if(fStatus == PolishedTeflonAirReflection)
593     G4cout << "PolishedTeflonAirReflection";      548     G4cout << "PolishedTeflonAirReflection";
594   else if(fStatus == PolishedTiOAirReflection)    549   else if(fStatus == PolishedTiOAirReflection)
595     G4cout << "PolishedTiOAirReflection";         550     G4cout << "PolishedTiOAirReflection";
596   else if(fStatus == PolishedTyvekAirReflectio    551   else if(fStatus == PolishedTyvekAirReflection)
597     G4cout << "PolishedTyvekAirReflection";       552     G4cout << "PolishedTyvekAirReflection";
598   else if(fStatus == PolishedVM2000AirReflecti    553   else if(fStatus == PolishedVM2000AirReflection)
599     G4cout << "PolishedVM2000AirReflection";      554     G4cout << "PolishedVM2000AirReflection";
600   else if(fStatus == PolishedVM2000GlueReflect    555   else if(fStatus == PolishedVM2000GlueReflection)
601     G4cout << "PolishedVM2000GlueReflection";     556     G4cout << "PolishedVM2000GlueReflection";
602   else if(fStatus == EtchedLumirrorAirReflecti    557   else if(fStatus == EtchedLumirrorAirReflection)
603     G4cout << "EtchedLumirrorAirReflection";      558     G4cout << "EtchedLumirrorAirReflection";
604   else if(fStatus == EtchedLumirrorGlueReflect    559   else if(fStatus == EtchedLumirrorGlueReflection)
605     G4cout << "EtchedLumirrorGlueReflection";     560     G4cout << "EtchedLumirrorGlueReflection";
606   else if(fStatus == EtchedAirReflection)         561   else if(fStatus == EtchedAirReflection)
607     G4cout << "EtchedAirReflection";              562     G4cout << "EtchedAirReflection";
608   else if(fStatus == EtchedTeflonAirReflection    563   else if(fStatus == EtchedTeflonAirReflection)
609     G4cout << "EtchedTeflonAirReflection";        564     G4cout << "EtchedTeflonAirReflection";
610   else if(fStatus == EtchedTiOAirReflection)      565   else if(fStatus == EtchedTiOAirReflection)
611     G4cout << "EtchedTiOAirReflection";           566     G4cout << "EtchedTiOAirReflection";
612   else if(fStatus == EtchedTyvekAirReflection)    567   else if(fStatus == EtchedTyvekAirReflection)
613     G4cout << "EtchedTyvekAirReflection";         568     G4cout << "EtchedTyvekAirReflection";
614   else if(fStatus == EtchedVM2000AirReflection    569   else if(fStatus == EtchedVM2000AirReflection)
615     G4cout << "EtchedVM2000AirReflection";        570     G4cout << "EtchedVM2000AirReflection";
616   else if(fStatus == EtchedVM2000GlueReflectio    571   else if(fStatus == EtchedVM2000GlueReflection)
617     G4cout << "EtchedVM2000GlueReflection";       572     G4cout << "EtchedVM2000GlueReflection";
618   else if(fStatus == GroundLumirrorAirReflecti    573   else if(fStatus == GroundLumirrorAirReflection)
619     G4cout << "GroundLumirrorAirReflection";      574     G4cout << "GroundLumirrorAirReflection";
620   else if(fStatus == GroundLumirrorGlueReflect    575   else if(fStatus == GroundLumirrorGlueReflection)
621     G4cout << "GroundLumirrorGlueReflection";     576     G4cout << "GroundLumirrorGlueReflection";
622   else if(fStatus == GroundAirReflection)         577   else if(fStatus == GroundAirReflection)
623     G4cout << "GroundAirReflection";              578     G4cout << "GroundAirReflection";
624   else if(fStatus == GroundTeflonAirReflection    579   else if(fStatus == GroundTeflonAirReflection)
625     G4cout << "GroundTeflonAirReflection";        580     G4cout << "GroundTeflonAirReflection";
626   else if(fStatus == GroundTiOAirReflection)      581   else if(fStatus == GroundTiOAirReflection)
627     G4cout << "GroundTiOAirReflection";           582     G4cout << "GroundTiOAirReflection";
628   else if(fStatus == GroundTyvekAirReflection)    583   else if(fStatus == GroundTyvekAirReflection)
629     G4cout << "GroundTyvekAirReflection";         584     G4cout << "GroundTyvekAirReflection";
630   else if(fStatus == GroundVM2000AirReflection    585   else if(fStatus == GroundVM2000AirReflection)
631     G4cout << "GroundVM2000AirReflection";        586     G4cout << "GroundVM2000AirReflection";
632   else if(fStatus == GroundVM2000GlueReflectio    587   else if(fStatus == GroundVM2000GlueReflection)
633     G4cout << "GroundVM2000GlueReflection";       588     G4cout << "GroundVM2000GlueReflection";
634   else if(fStatus == Absorption)                  589   else if(fStatus == Absorption)
635     G4cout << "Absorption";                       590     G4cout << "Absorption";
636   else if(fStatus == Detection)                   591   else if(fStatus == Detection)
637     G4cout << "Detection";                        592     G4cout << "Detection";
638   else if(fStatus == NotAtBoundary)               593   else if(fStatus == NotAtBoundary)
639     G4cout << "NotAtBoundary";                    594     G4cout << "NotAtBoundary";
640   else if(fStatus == SameMaterial)                595   else if(fStatus == SameMaterial)
641     G4cout << "SameMaterial";                     596     G4cout << "SameMaterial";
642   else if(fStatus == StepTooSmall)                597   else if(fStatus == StepTooSmall)
643     G4cout << "StepTooSmall";                     598     G4cout << "StepTooSmall";
644   else if(fStatus == NoRINDEX)                    599   else if(fStatus == NoRINDEX)
645     G4cout << "NoRINDEX";                         600     G4cout << "NoRINDEX";
646   else if(fStatus == Dichroic)                    601   else if(fStatus == Dichroic)
647     G4cout << "Dichroic Transmission";            602     G4cout << "Dichroic Transmission";
648   else if(fStatus == CoatedDielectricReflectio << 
649     G4cout << "Coated Dielectric Reflection";  << 
650   else if(fStatus == CoatedDielectricRefractio << 
651     G4cout << "Coated Dielectric Refraction";  << 
652   else if(fStatus == CoatedDielectricFrustrate << 
653     G4cout << "Coated Dielectric Frustrated Tr << 
654                                                << 
655   G4cout << " ***" << G4endl;                     603   G4cout << " ***" << G4endl;
656 }                                                 604 }
657                                                   605 
658 //....oooOO0OOooo........oooOO0OOooo........oo    606 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
659 G4ThreeVector G4OpBoundaryProcess::GetFacetNor    607 G4ThreeVector G4OpBoundaryProcess::GetFacetNormal(
660   const G4ThreeVector& momentum, const G4Three    608   const G4ThreeVector& momentum, const G4ThreeVector& normal) const
661 {                                                 609 {
662   G4ThreeVector facetNormal;                      610   G4ThreeVector facetNormal;
663   if(fModel == unified || fModel == LUT || fMo    611   if(fModel == unified || fModel == LUT || fModel == DAVIS)
664   {                                               612   {
665     /* This function codes alpha to a random v    613     /* This function codes alpha to a random value taken from the
666     distribution p(alpha) = g(alpha; 0, sigma_    614     distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
667     for alpha > 0 and alpha < 90, where g(alph    615     for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha) is a
668     gaussian distribution with mean 0 and stan    616     gaussian distribution with mean 0 and standard deviation sigma_alpha.  */
669                                                   617 
670     G4double sigma_alpha = 0.0;                   618     G4double sigma_alpha = 0.0;
671     if(fOpticalSurface)                           619     if(fOpticalSurface)
672       sigma_alpha = fOpticalSurface->GetSigmaA    620       sigma_alpha = fOpticalSurface->GetSigmaAlpha();
673     if(sigma_alpha == 0.0)                        621     if(sigma_alpha == 0.0)
674     {                                             622     {
675       return normal;                              623       return normal;
676     }                                             624     }
677                                                   625 
678     G4double f_max = std::min(1.0, 4. * sigma_    626     G4double f_max = std::min(1.0, 4. * sigma_alpha);
679     G4double alpha, phi, sinAlpha;                627     G4double alpha, phi, sinAlpha;
680                                                   628 
681     do                                            629     do
682     {  // Loop checking, 13-Aug-2015, Peter Gu    630     {  // Loop checking, 13-Aug-2015, Peter Gumplinger
683       do                                          631       do
684       {  // Loop checking, 13-Aug-2015, Peter     632       {  // Loop checking, 13-Aug-2015, Peter Gumplinger
685         alpha    = G4RandGauss::shoot(0.0, sig    633         alpha    = G4RandGauss::shoot(0.0, sigma_alpha);
686         sinAlpha = std::sin(alpha);               634         sinAlpha = std::sin(alpha);
687       } while(G4UniformRand() * f_max > sinAlp    635       } while(G4UniformRand() * f_max > sinAlpha || alpha >= halfpi);
688                                                   636 
689       phi = G4UniformRand() * twopi;              637       phi = G4UniformRand() * twopi;
690       facetNormal.set(sinAlpha * std::cos(phi)    638       facetNormal.set(sinAlpha * std::cos(phi), sinAlpha * std::sin(phi),
691                       std::cos(alpha));           639                       std::cos(alpha));
692       facetNormal.rotateUz(normal);               640       facetNormal.rotateUz(normal);
693     } while(momentum * facetNormal >= 0.0);       641     } while(momentum * facetNormal >= 0.0);
694   }                                               642   }
695   else                                            643   else
696   {                                               644   {
697     G4double polish = 1.0;                        645     G4double polish = 1.0;
698     if(fOpticalSurface)                           646     if(fOpticalSurface)
699       polish = fOpticalSurface->GetPolish();      647       polish = fOpticalSurface->GetPolish();
700     if(polish < 1.0)                              648     if(polish < 1.0)
701     {                                             649     {
702       do                                          650       do
703       {  // Loop checking, 13-Aug-2015, Peter     651       {  // Loop checking, 13-Aug-2015, Peter Gumplinger
704         G4ThreeVector smear;                      652         G4ThreeVector smear;
705         do                                        653         do
706         {  // Loop checking, 13-Aug-2015, Pete    654         {  // Loop checking, 13-Aug-2015, Peter Gumplinger
707           smear.setX(2. * G4UniformRand() - 1.    655           smear.setX(2. * G4UniformRand() - 1.);
708           smear.setY(2. * G4UniformRand() - 1.    656           smear.setY(2. * G4UniformRand() - 1.);
709           smear.setZ(2. * G4UniformRand() - 1.    657           smear.setZ(2. * G4UniformRand() - 1.);
710         } while(smear.mag2() > 1.0);              658         } while(smear.mag2() > 1.0);
711         facetNormal = normal + (1. - polish) *    659         facetNormal = normal + (1. - polish) * smear;
712       } while(momentum * facetNormal >= 0.0);     660       } while(momentum * facetNormal >= 0.0);
713       facetNormal = facetNormal.unit();           661       facetNormal = facetNormal.unit();
714     }                                             662     }
715     else                                          663     else
716     {                                             664     {
717       facetNormal = normal;                       665       facetNormal = normal;
718     }                                             666     }
719   }                                               667   }
720   return facetNormal;                             668   return facetNormal;
721 }                                                 669 }
722                                                   670 
723 //....oooOO0OOooo........oooOO0OOooo........oo    671 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
724 void G4OpBoundaryProcess::DielectricMetal()       672 void G4OpBoundaryProcess::DielectricMetal()
725 {                                                 673 {
726   G4int n = 0;                                    674   G4int n = 0;
727   G4double rand;                                  675   G4double rand;
728   G4ThreeVector A_trans;                          676   G4ThreeVector A_trans;
729                                                   677 
730   do                                              678   do
731   {                                               679   {
732     ++n;                                          680     ++n;
733     rand = G4UniformRand();                       681     rand = G4UniformRand();
734     if(rand > fReflectivity && n == 1)            682     if(rand > fReflectivity && n == 1)
735     {                                             683     {
736       if(rand > fReflectivity + fTransmittance    684       if(rand > fReflectivity + fTransmittance)
737       {                                           685       {
738         DoAbsorption();                           686         DoAbsorption();
739       }                                           687       }
740       else                                        688       else
741       {                                           689       {
742         fStatus          = Transmission;          690         fStatus          = Transmission;
743         fNewMomentum     = fOldMomentum;          691         fNewMomentum     = fOldMomentum;
744         fNewPolarization = fOldPolarization;      692         fNewPolarization = fOldPolarization;
745       }                                           693       }
746       break;                                      694       break;
747     }                                             695     }
748     else                                          696     else
749     {                                             697     {
750       if(fRealRIndexMPV && fImagRIndexMPV)        698       if(fRealRIndexMPV && fImagRIndexMPV)
751       {                                           699       {
752         if(n > 1)                                 700         if(n > 1)
753         {                                         701         {
754           CalculateReflectivity();                702           CalculateReflectivity();
755           if(!G4BooleanRand(fReflectivity))       703           if(!G4BooleanRand(fReflectivity))
756           {                                       704           {
757             DoAbsorption();                       705             DoAbsorption();
758             break;                                706             break;
759           }                                       707           }
760         }                                         708         }
761       }                                           709       }
762       if(fModel == glisur || fFinish == polish    710       if(fModel == glisur || fFinish == polished)
763       {                                           711       {
764         DoReflection();                           712         DoReflection();
765       }                                           713       }
766       else                                        714       else
767       {                                           715       {
768         if(n == 1)                                716         if(n == 1)
769           ChooseReflection();                     717           ChooseReflection();
770         if(fStatus == LambertianReflection)       718         if(fStatus == LambertianReflection)
771         {                                         719         {
772           DoReflection();                         720           DoReflection();
773         }                                         721         }
774         else if(fStatus == BackScattering)        722         else if(fStatus == BackScattering)
775         {                                         723         {
776           fNewMomentum     = -fOldMomentum;       724           fNewMomentum     = -fOldMomentum;
777           fNewPolarization = -fOldPolarization    725           fNewPolarization = -fOldPolarization;
778         }                                         726         }
779         else                                      727         else
780         {                                         728         {
781           if(fStatus == LobeReflection)           729           if(fStatus == LobeReflection)
782           {                                       730           {
783             if(!fRealRIndexMPV || !fImagRIndex    731             if(!fRealRIndexMPV || !fImagRIndexMPV)
784             {                                     732             {
785               fFacetNormal = GetFacetNormal(fO    733               fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
786             }                                     734             }
787             // else                            << 735             //else
788             //  case of complex rindex needs t    736             //  case of complex rindex needs to be implemented
789           }                                       737           }
790           fNewMomentum =                          738           fNewMomentum =
791             fOldMomentum - 2. * fOldMomentum *    739             fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
792                                                   740 
793           if(f_iTE > 0 && f_iTM > 0)              741           if(f_iTE > 0 && f_iTM > 0)
794           {                                       742           {
795             fNewPolarization =                    743             fNewPolarization =
796               -fOldPolarization +                 744               -fOldPolarization +
797               (2. * fOldPolarization * fFacetN    745               (2. * fOldPolarization * fFacetNormal * fFacetNormal);
798           }                                       746           }
799           else if(f_iTE > 0)                      747           else if(f_iTE > 0)
800           {                                       748           {
801             A_trans = (fSint1 > 0.0) ? fOldMom    749             A_trans = (fSint1 > 0.0) ? fOldMomentum.cross(fFacetNormal).unit()
802                                      : fOldPol    750                                      : fOldPolarization;
803             fNewPolarization = -A_trans;          751             fNewPolarization = -A_trans;
804           }                                       752           }
805           else if(f_iTM > 0)                      753           else if(f_iTM > 0)
806           {                                       754           {
807             fNewPolarization =                    755             fNewPolarization =
808               -fNewMomentum.cross(A_trans).uni    756               -fNewMomentum.cross(A_trans).unit();  // = -A_paral
809           }                                       757           }
810         }                                         758         }
811       }                                           759       }
812       fOldMomentum     = fNewMomentum;            760       fOldMomentum     = fNewMomentum;
813       fOldPolarization = fNewPolarization;        761       fOldPolarization = fNewPolarization;
814     }                                             762     }
815     // Loop checking, 13-Aug-2015, Peter Gumpl    763     // Loop checking, 13-Aug-2015, Peter Gumplinger
816   } while(fNewMomentum * fGlobalNormal < 0.0);    764   } while(fNewMomentum * fGlobalNormal < 0.0);
817 }                                                 765 }
818                                                   766 
819 //....oooOO0OOooo........oooOO0OOooo........oo    767 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
820 void G4OpBoundaryProcess::DielectricLUT()         768 void G4OpBoundaryProcess::DielectricLUT()
821 {                                                 769 {
822   G4int thetaIndex, phiIndex;                     770   G4int thetaIndex, phiIndex;
823   G4double angularDistVal, thetaRad, phiRad;      771   G4double angularDistVal, thetaRad, phiRad;
824   G4ThreeVector perpVectorTheta, perpVectorPhi    772   G4ThreeVector perpVectorTheta, perpVectorPhi;
825                                                   773 
826   fStatus = G4OpBoundaryProcessStatus(            774   fStatus = G4OpBoundaryProcessStatus(
827     G4int(fFinish) + (G4int(NoRINDEX) - G4int(    775     G4int(fFinish) + (G4int(NoRINDEX) - G4int(groundbackpainted)));
828                                                   776 
829   G4int thetaIndexMax = fOpticalSurface->GetTh    777   G4int thetaIndexMax = fOpticalSurface->GetThetaIndexMax();
830   G4int phiIndexMax   = fOpticalSurface->GetPh    778   G4int phiIndexMax   = fOpticalSurface->GetPhiIndexMax();
831                                                   779 
832   G4double rand;                                  780   G4double rand;
833                                                   781 
834   do                                              782   do
835   {                                               783   {
836     rand = G4UniformRand();                       784     rand = G4UniformRand();
837     if(rand > fReflectivity)                      785     if(rand > fReflectivity)
838     {                                             786     {
839       if(rand > fReflectivity + fTransmittance    787       if(rand > fReflectivity + fTransmittance)
840       {                                           788       {
841         DoAbsorption();                           789         DoAbsorption();
842       }                                           790       }
843       else                                        791       else
844       {                                           792       {
845         fStatus          = Transmission;          793         fStatus          = Transmission;
846         fNewMomentum     = fOldMomentum;          794         fNewMomentum     = fOldMomentum;
847         fNewPolarization = fOldPolarization;      795         fNewPolarization = fOldPolarization;
848       }                                           796       }
849       break;                                      797       break;
850     }                                             798     }
851     else                                          799     else
852     {                                             800     {
853       // Calculate Angle between Normal and Ph    801       // Calculate Angle between Normal and Photon Momentum
854       G4double anglePhotonToNormal = fOldMomen    802       G4double anglePhotonToNormal = fOldMomentum.angle(-fGlobalNormal);
855       // Round to closest integer: LBNL model     803       // Round to closest integer: LBNL model array has 91 values
856       G4int angleIncident = (G4int)std::lrint( << 804       G4int angleIncident = std::lrint(anglePhotonToNormal / CLHEP::deg);
857                                                   805 
858       // Take random angles THETA and PHI,        806       // Take random angles THETA and PHI,
859       // and see if below Probability - if not    807       // and see if below Probability - if not - Redo
860       do                                          808       do
861       {                                           809       {
862         thetaIndex = (G4int)G4RandFlat::shootI << 810         thetaIndex = G4RandFlat::shootInt(thetaIndexMax - 1);
863         phiIndex   = (G4int)G4RandFlat::shootI << 811         phiIndex   = G4RandFlat::shootInt(phiIndexMax - 1);
864         // Find probability with the new indec    812         // Find probability with the new indeces from LUT
865         angularDistVal = fOpticalSurface->GetA    813         angularDistVal = fOpticalSurface->GetAngularDistributionValue(
866           angleIncident, thetaIndex, phiIndex)    814           angleIncident, thetaIndex, phiIndex);
867         // Loop checking, 13-Aug-2015, Peter G    815         // Loop checking, 13-Aug-2015, Peter Gumplinger
868       } while(!G4BooleanRand(angularDistVal));    816       } while(!G4BooleanRand(angularDistVal));
869                                                   817 
870       thetaRad = G4double(-90 + 4 * thetaIndex    818       thetaRad = G4double(-90 + 4 * thetaIndex) * pi / 180.;
871       phiRad   = G4double(-90 + 5 * phiIndex)     819       phiRad   = G4double(-90 + 5 * phiIndex) * pi / 180.;
872       // Rotate Photon Momentum in Theta, then    820       // Rotate Photon Momentum in Theta, then in Phi
873       fNewMomentum = -fOldMomentum;               821       fNewMomentum = -fOldMomentum;
874                                                   822 
875       perpVectorTheta = fNewMomentum.cross(fGl    823       perpVectorTheta = fNewMomentum.cross(fGlobalNormal);
876       if(perpVectorTheta.mag() < fCarTolerance    824       if(perpVectorTheta.mag() < fCarTolerance)
877       {                                           825       {
878         perpVectorTheta = fNewMomentum.orthogo    826         perpVectorTheta = fNewMomentum.orthogonal();
879       }                                           827       }
880       fNewMomentum =                              828       fNewMomentum =
881         fNewMomentum.rotate(anglePhotonToNorma    829         fNewMomentum.rotate(anglePhotonToNormal - thetaRad, perpVectorTheta);
882       perpVectorPhi = perpVectorTheta.cross(fN    830       perpVectorPhi = perpVectorTheta.cross(fNewMomentum);
883       fNewMomentum  = fNewMomentum.rotate(-phi    831       fNewMomentum  = fNewMomentum.rotate(-phiRad, perpVectorPhi);
884                                                   832 
885       // Rotate Polarization too:                 833       // Rotate Polarization too:
886       fFacetNormal     = (fNewMomentum - fOldM    834       fFacetNormal     = (fNewMomentum - fOldMomentum).unit();
887       fNewPolarization = -fOldPolarization +      835       fNewPolarization = -fOldPolarization +
888                          (2. * fOldPolarizatio    836                          (2. * fOldPolarization * fFacetNormal * fFacetNormal);
889     }                                             837     }
890     // Loop checking, 13-Aug-2015, Peter Gumpl    838     // Loop checking, 13-Aug-2015, Peter Gumplinger
891   } while(fNewMomentum * fGlobalNormal <= 0.0)    839   } while(fNewMomentum * fGlobalNormal <= 0.0);
892 }                                                 840 }
893                                                   841 
894 //....oooOO0OOooo........oooOO0OOooo........oo    842 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
895 void G4OpBoundaryProcess::DielectricLUTDAVIS()    843 void G4OpBoundaryProcess::DielectricLUTDAVIS()
896 {                                                 844 {
897   G4int angindex, random, angleIncident;          845   G4int angindex, random, angleIncident;
898   G4double reflectivityValue, elevation, azimu    846   G4double reflectivityValue, elevation, azimuth;
899   G4double anglePhotonToNormal;                   847   G4double anglePhotonToNormal;
900                                                   848 
901   G4int lutbin  = fOpticalSurface->GetLUTbins(    849   G4int lutbin  = fOpticalSurface->GetLUTbins();
902   G4double rand = G4UniformRand();                850   G4double rand = G4UniformRand();
903                                                   851 
904   G4double sinEl;                                 852   G4double sinEl;
905   G4ThreeVector u, vNorm, w;                      853   G4ThreeVector u, vNorm, w;
906                                                   854 
907   do                                              855   do
908   {                                               856   {
909     anglePhotonToNormal = fOldMomentum.angle(-    857     anglePhotonToNormal = fOldMomentum.angle(-fGlobalNormal);
910                                                   858 
911     // Davis model has 90 reflection bins: rou    859     // Davis model has 90 reflection bins: round down
912     // don't allow angleIncident to be 90 for     860     // don't allow angleIncident to be 90 for anglePhotonToNormal close to 90
913     angleIncident = std::min(                  << 861     angleIncident = std::min(static_cast<G4int>(
914       static_cast<G4int>(std::floor(anglePhoto << 862       std::floor(anglePhotonToNormal / CLHEP::deg)), 89);
915     reflectivityValue = fOpticalSurface->GetRe    863     reflectivityValue = fOpticalSurface->GetReflectivityLUTValue(angleIncident);
916                                                   864 
917     if(rand > reflectivityValue)                  865     if(rand > reflectivityValue)
918     {                                             866     {
919       if(fEfficiency > 0.)                        867       if(fEfficiency > 0.)
920       {                                           868       {
921         DoAbsorption();                           869         DoAbsorption();
922         break;                                    870         break;
923       }                                           871       }
924       else                                        872       else
925       {                                           873       {
926         fStatus = Transmission;                   874         fStatus = Transmission;
927                                                   875 
928         if(angleIncident <= 0.01)                 876         if(angleIncident <= 0.01)
929         {                                         877         {
930           fNewMomentum = fOldMomentum;            878           fNewMomentum = fOldMomentum;
931           break;                                  879           break;
932         }                                         880         }
933                                                   881 
934         do                                        882         do
935         {                                         883         {
936           random = (G4int)G4RandFlat::shootInt << 884           random = G4RandFlat::shootInt(1, lutbin + 1);
937           angindex =                              885           angindex =
938             (((random * 2) - 1)) + angleIncide    886             (((random * 2) - 1)) + angleIncident * lutbin * 2 + 3640000;
939                                                   887 
940           azimuth =                               888           azimuth =
941             fOpticalSurface->GetAngularDistrib    889             fOpticalSurface->GetAngularDistributionValueLUT(angindex - 1);
942           elevation = fOpticalSurface->GetAngu    890           elevation = fOpticalSurface->GetAngularDistributionValueLUT(angindex);
943         } while(elevation == 0. && azimuth ==     891         } while(elevation == 0. && azimuth == 0.);
944                                                   892 
945         sinEl = std::sin(elevation);              893         sinEl = std::sin(elevation);
946         vNorm = (fGlobalNormal.cross(fOldMomen    894         vNorm = (fGlobalNormal.cross(fOldMomentum)).unit();
947         u     = vNorm.cross(fGlobalNormal) * (    895         u     = vNorm.cross(fGlobalNormal) * (sinEl * std::cos(azimuth));
948         vNorm *= (sinEl * std::sin(azimuth));     896         vNorm *= (sinEl * std::sin(azimuth));
949         // fGlobalNormal shouldn't be modified    897         // fGlobalNormal shouldn't be modified here
950         w            = (fGlobalNormal *= std::    898         w            = (fGlobalNormal *= std::cos(elevation));
951         fNewMomentum = u + vNorm + w;             899         fNewMomentum = u + vNorm + w;
952                                                   900 
953         // Rotate Polarization too:               901         // Rotate Polarization too:
954         fFacetNormal     = (fNewMomentum - fOl    902         fFacetNormal     = (fNewMomentum - fOldMomentum).unit();
955         fNewPolarization = -fOldPolarization +    903         fNewPolarization = -fOldPolarization + (2. * fOldPolarization *
956                                                   904                                                 fFacetNormal * fFacetNormal);
957       }                                           905       }
958     }                                             906     }
959     else                                          907     else
960     {                                             908     {
961       fStatus = LobeReflection;                   909       fStatus = LobeReflection;
962                                                   910 
963       if(angleIncident == 0)                      911       if(angleIncident == 0)
964       {                                           912       {
965         fNewMomentum = -fOldMomentum;             913         fNewMomentum = -fOldMomentum;
966         break;                                    914         break;
967       }                                           915       }
968                                                   916 
969       do                                          917       do
970       {                                           918       {
971         random   = (G4int)G4RandFlat::shootInt << 919         random   = G4RandFlat::shootInt(1, lutbin + 1);
972         angindex = (((random * 2) - 1)) + (ang    920         angindex = (((random * 2) - 1)) + (angleIncident - 1) * lutbin * 2;
973                                                   921 
974         azimuth = fOpticalSurface->GetAngularD    922         azimuth = fOpticalSurface->GetAngularDistributionValueLUT(angindex - 1);
975         elevation = fOpticalSurface->GetAngula    923         elevation = fOpticalSurface->GetAngularDistributionValueLUT(angindex);
976       } while(elevation == 0. && azimuth == 0.    924       } while(elevation == 0. && azimuth == 0.);
977                                                   925 
978       sinEl = std::sin(elevation);                926       sinEl = std::sin(elevation);
979       vNorm = (fGlobalNormal.cross(fOldMomentu    927       vNorm = (fGlobalNormal.cross(fOldMomentum)).unit();
980       u     = vNorm.cross(fGlobalNormal) * (si    928       u     = vNorm.cross(fGlobalNormal) * (sinEl * std::cos(azimuth));
981       vNorm *= (sinEl * std::sin(azimuth));       929       vNorm *= (sinEl * std::sin(azimuth));
982       // fGlobalNormal shouldn't be modified h    930       // fGlobalNormal shouldn't be modified here
983       w = (fGlobalNormal *= std::cos(elevation    931       w = (fGlobalNormal *= std::cos(elevation));
984                                                   932 
985       fNewMomentum = u + vNorm + w;               933       fNewMomentum = u + vNorm + w;
986                                                   934 
987       // Rotate Polarization too: (needs revis    935       // Rotate Polarization too: (needs revision)
988       fNewPolarization = fOldPolarization;        936       fNewPolarization = fOldPolarization;
989     }                                             937     }
990   } while(fNewMomentum * fGlobalNormal <= 0.0)    938   } while(fNewMomentum * fGlobalNormal <= 0.0);
991 }                                                 939 }
992                                                   940 
993 //....oooOO0OOooo........oooOO0OOooo........oo    941 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
994 void G4OpBoundaryProcess::DielectricDichroic()    942 void G4OpBoundaryProcess::DielectricDichroic()
995 {                                                 943 {
996   // Calculate Angle between Normal and Photon    944   // Calculate Angle between Normal and Photon Momentum
997   G4double anglePhotonToNormal = fOldMomentum.    945   G4double anglePhotonToNormal = fOldMomentum.angle(-fGlobalNormal);
998                                                   946 
999   // Round it to closest integer                  947   // Round it to closest integer
1000   G4double angleIncident = std::floor(180. /     948   G4double angleIncident = std::floor(180. / pi * anglePhotonToNormal + 0.5);
1001                                                  949 
1002   if(!fDichroicVector)                           950   if(!fDichroicVector)
1003   {                                              951   {
1004     if(fOpticalSurface)                          952     if(fOpticalSurface)
1005       fDichroicVector = fOpticalSurface->GetD    953       fDichroicVector = fOpticalSurface->GetDichroicVector();
1006   }                                              954   }
1007                                                  955 
1008   if(fDichroicVector)                            956   if(fDichroicVector)
1009   {                                              957   {
1010     G4double wavelength = h_Planck * c_light     958     G4double wavelength = h_Planck * c_light / fPhotonMomentum;
1011     fTransmittance      = fDichroicVector->Va    959     fTransmittance      = fDichroicVector->Value(wavelength / nm, angleIncident,
1012                                             i    960                                             idx_dichroicX, idx_dichroicY) *
1013                      perCent;                    961                      perCent;
1014     //   G4cout << "wavelength: " << std::flo    962     //   G4cout << "wavelength: " << std::floor(wavelength/nm)
1015     //                            << "nm" <<     963     //                            << "nm" << G4endl;
1016     //   G4cout << "Incident angle: " << angl    964     //   G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
1017     //   G4cout << "Transmittance: "             965     //   G4cout << "Transmittance: "
1018     //          << std::floor(fTransmittance/    966     //          << std::floor(fTransmittance/perCent) << "%" << G4endl;
1019   }                                              967   }
1020   else                                           968   else
1021   {                                              969   {
1022     G4ExceptionDescription ed;                   970     G4ExceptionDescription ed;
1023     ed << " G4OpBoundaryProcess/DielectricDic    971     ed << " G4OpBoundaryProcess/DielectricDichroic(): "
1024        << " The dichroic surface has no G4Phy    972        << " The dichroic surface has no G4Physics2DVector" << G4endl;
1025     G4Exception("G4OpBoundaryProcess::Dielect    973     G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
1026                 FatalException, ed,              974                 FatalException, ed,
1027                 "A dichroic surface must have    975                 "A dichroic surface must have an associated G4Physics2DVector");
1028   }                                              976   }
1029                                                  977 
1030   if(!G4BooleanRand(fTransmittance))             978   if(!G4BooleanRand(fTransmittance))
1031   {  // Not transmitted, so reflect              979   {  // Not transmitted, so reflect
1032     if(fModel == glisur || fFinish == polishe    980     if(fModel == glisur || fFinish == polished)
1033     {                                            981     {
1034       DoReflection();                            982       DoReflection();
1035     }                                            983     }
1036     else                                         984     else
1037     {                                            985     {
1038       ChooseReflection();                        986       ChooseReflection();
1039       if(fStatus == LambertianReflection)        987       if(fStatus == LambertianReflection)
1040       {                                          988       {
1041         DoReflection();                          989         DoReflection();
1042       }                                          990       }
1043       else if(fStatus == BackScattering)         991       else if(fStatus == BackScattering)
1044       {                                          992       {
1045         fNewMomentum     = -fOldMomentum;        993         fNewMomentum     = -fOldMomentum;
1046         fNewPolarization = -fOldPolarization;    994         fNewPolarization = -fOldPolarization;
1047       }                                          995       }
1048       else                                       996       else
1049       {                                          997       {
1050         G4double PdotN, EdotN;                   998         G4double PdotN, EdotN;
1051         do                                       999         do
1052         {                                        1000         {
1053           if(fStatus == LobeReflection)          1001           if(fStatus == LobeReflection)
1054           {                                      1002           {
1055             fFacetNormal = GetFacetNormal(fOl    1003             fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1056           }                                      1004           }
1057           PdotN        = fOldMomentum * fFace    1005           PdotN        = fOldMomentum * fFacetNormal;
1058           fNewMomentum = fOldMomentum - (2. *    1006           fNewMomentum = fOldMomentum - (2. * PdotN) * fFacetNormal;
1059           // Loop checking, 13-Aug-2015, Pete    1007           // Loop checking, 13-Aug-2015, Peter Gumplinger
1060         } while(fNewMomentum * fGlobalNormal     1008         } while(fNewMomentum * fGlobalNormal <= 0.0);
1061                                                  1009 
1062         EdotN            = fOldPolarization *    1010         EdotN            = fOldPolarization * fFacetNormal;
1063         fNewPolarization = -fOldPolarization     1011         fNewPolarization = -fOldPolarization + (2. * EdotN) * fFacetNormal;
1064       }                                          1012       }
1065     }                                            1013     }
1066   }                                              1014   }
1067   else                                           1015   else
1068   {                                              1016   {
1069     fStatus          = Dichroic;                 1017     fStatus          = Dichroic;
1070     fNewMomentum     = fOldMomentum;             1018     fNewMomentum     = fOldMomentum;
1071     fNewPolarization = fOldPolarization;         1019     fNewPolarization = fOldPolarization;
1072   }                                              1020   }
1073 }                                                1021 }
1074                                                  1022 
1075 //....oooOO0OOooo........oooOO0OOooo........o    1023 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1076 void G4OpBoundaryProcess::DielectricDielectri    1024 void G4OpBoundaryProcess::DielectricDielectric()
1077 {                                                1025 {
1078   G4bool inside = false;                         1026   G4bool inside = false;
1079   G4bool swap   = false;                         1027   G4bool swap   = false;
1080                                                  1028 
1081   if(fFinish == polished)                        1029   if(fFinish == polished)
1082   {                                              1030   {
1083     fFacetNormal = fGlobalNormal;                1031     fFacetNormal = fGlobalNormal;
1084   }                                              1032   }
1085   else                                           1033   else
1086   {                                              1034   {
1087     fFacetNormal = GetFacetNormal(fOldMomentu    1035     fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1088   }                                              1036   }
1089   G4double cost1 = -fOldMomentum * fFacetNorm    1037   G4double cost1 = -fOldMomentum * fFacetNormal;
1090   G4double cost2 = 0.;                           1038   G4double cost2 = 0.;
1091   G4double sint2 = 0.;                           1039   G4double sint2 = 0.;
1092                                                  1040 
1093   G4bool surfaceRoughnessCriterionPass = true    1041   G4bool surfaceRoughnessCriterionPass = true;
1094   if(fSurfaceRoughness != 0. && fRindex1 > fR    1042   if(fSurfaceRoughness != 0. && fRindex1 > fRindex2)
1095   {                                              1043   {
1096     G4double wavelength                = h_Pl    1044     G4double wavelength                = h_Planck * c_light / fPhotonMomentum;
1097     G4double surfaceRoughnessCriterion = std:    1045     G4double surfaceRoughnessCriterion = std::exp(-std::pow(
1098       (4. * pi * fSurfaceRoughness * fRindex1    1046       (4. * pi * fSurfaceRoughness * fRindex1 * cost1 / wavelength), 2));
1099     surfaceRoughnessCriterionPass = G4Boolean    1047     surfaceRoughnessCriterionPass = G4BooleanRand(surfaceRoughnessCriterion);
1100   }                                              1048   }
1101                                                  1049 
1102 leap:                                            1050 leap:
1103                                                  1051 
1104   G4bool through = false;                        1052   G4bool through = false;
1105   G4bool done    = false;                        1053   G4bool done    = false;
1106                                                  1054 
1107   G4ThreeVector A_trans, A_paral, E1pp, E1pl;    1055   G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1108   G4double E1_perp, E1_parl;                     1056   G4double E1_perp, E1_parl;
1109   G4double s1, s2, E2_perp, E2_parl, E2_total    1057   G4double s1, s2, E2_perp, E2_parl, E2_total, transCoeff;
1110   G4double E2_abs, C_parl, C_perp;               1058   G4double E2_abs, C_parl, C_perp;
1111   G4double alpha;                                1059   G4double alpha;
1112                                                  1060 
1113   do                                             1061   do
1114   {                                              1062   {
1115     if(through)                                  1063     if(through)
1116     {                                            1064     {
1117       swap          = !swap;                     1065       swap          = !swap;
1118       through       = false;                     1066       through       = false;
1119       fGlobalNormal = -fGlobalNormal;            1067       fGlobalNormal = -fGlobalNormal;
1120       G4SwapPtr(fMaterial1, fMaterial2);         1068       G4SwapPtr(fMaterial1, fMaterial2);
1121       G4SwapObj(&fRindex1, &fRindex2);           1069       G4SwapObj(&fRindex1, &fRindex2);
1122     }                                            1070     }
1123                                                  1071 
1124     if(fFinish == polished)                      1072     if(fFinish == polished)
1125     {                                            1073     {
1126       fFacetNormal = fGlobalNormal;              1074       fFacetNormal = fGlobalNormal;
1127     }                                            1075     }
1128     else                                         1076     else
1129     {                                            1077     {
1130       fFacetNormal = GetFacetNormal(fOldMomen    1078       fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1131     }                                            1079     }
1132                                                  1080 
1133     cost1 = -fOldMomentum * fFacetNormal;        1081     cost1 = -fOldMomentum * fFacetNormal;
1134     if(std::abs(cost1) < 1.0 - fCarTolerance)    1082     if(std::abs(cost1) < 1.0 - fCarTolerance)
1135     {                                            1083     {
1136       fSint1 = std::sqrt(1. - cost1 * cost1);    1084       fSint1 = std::sqrt(1. - cost1 * cost1);
1137       sint2  = fSint1 * fRindex1 / fRindex2;     1085       sint2  = fSint1 * fRindex1 / fRindex2;  // *** Snell's Law ***
1138       // this isn't a sine as we might expect    1086       // this isn't a sine as we might expect from the name; can be > 1
1139     }                                            1087     }
1140     else                                         1088     else
1141     {                                            1089     {
1142       fSint1 = 0.0;                              1090       fSint1 = 0.0;
1143       sint2  = 0.0;                              1091       sint2  = 0.0;
1144     }                                            1092     }
1145                                                  1093 
1146     // TOTAL INTERNAL REFLECTION                 1094     // TOTAL INTERNAL REFLECTION
1147     if(sint2 >= 1.0)                             1095     if(sint2 >= 1.0)
1148     {                                            1096     {
1149       swap = false;                              1097       swap = false;
1150                                                  1098 
1151       fStatus = TotalInternalReflection;         1099       fStatus = TotalInternalReflection;
1152       if(!surfaceRoughnessCriterionPass)         1100       if(!surfaceRoughnessCriterionPass)
1153         fStatus = LambertianReflection;          1101         fStatus = LambertianReflection;
1154       if(fModel == unified && fFinish != poli    1102       if(fModel == unified && fFinish != polished)
1155         ChooseReflection();                      1103         ChooseReflection();
1156       if(fStatus == LambertianReflection)        1104       if(fStatus == LambertianReflection)
1157       {                                          1105       {
1158         DoReflection();                          1106         DoReflection();
1159       }                                          1107       }
1160       else if(fStatus == BackScattering)         1108       else if(fStatus == BackScattering)
1161       {                                          1109       {
1162         fNewMomentum     = -fOldMomentum;        1110         fNewMomentum     = -fOldMomentum;
1163         fNewPolarization = -fOldPolarization;    1111         fNewPolarization = -fOldPolarization;
1164       }                                          1112       }
1165       else                                       1113       else
1166       {                                          1114       {
1167         fNewMomentum =                           1115         fNewMomentum =
1168           fOldMomentum - 2. * fOldMomentum *     1116           fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
1169         fNewPolarization = -fOldPolarization     1117         fNewPolarization = -fOldPolarization + (2. * fOldPolarization *
1170                                                  1118                                                 fFacetNormal * fFacetNormal);
1171       }                                          1119       }
1172     }                                            1120     }
1173     // NOT TIR                                   1121     // NOT TIR
1174     else if(sint2 < 1.0)                         1122     else if(sint2 < 1.0)
1175     {                                            1123     {
1176       // Calculate amplitude for transmission    1124       // Calculate amplitude for transmission (Q = P x N)
1177       if(cost1 > 0.0)                            1125       if(cost1 > 0.0)
1178       {                                          1126       {
1179         cost2 = std::sqrt(1. - sint2 * sint2)    1127         cost2 = std::sqrt(1. - sint2 * sint2);
1180       }                                          1128       }
1181       else                                       1129       else
1182       {                                          1130       {
1183         cost2 = -std::sqrt(1. - sint2 * sint2    1131         cost2 = -std::sqrt(1. - sint2 * sint2);
1184       }                                          1132       }
1185                                                  1133 
1186       if(fSint1 > 0.0)                           1134       if(fSint1 > 0.0)
1187       {                                          1135       {
1188         A_trans = (fOldMomentum.cross(fFacetN    1136         A_trans = (fOldMomentum.cross(fFacetNormal)).unit();
1189         E1_perp = fOldPolarization * A_trans;    1137         E1_perp = fOldPolarization * A_trans;
1190         E1pp    = E1_perp * A_trans;             1138         E1pp    = E1_perp * A_trans;
1191         E1pl    = fOldPolarization - E1pp;       1139         E1pl    = fOldPolarization - E1pp;
1192         E1_parl = E1pl.mag();                    1140         E1_parl = E1pl.mag();
1193       }                                          1141       }
1194       else                                       1142       else
1195       {                                          1143       {
1196         A_trans = fOldPolarization;              1144         A_trans = fOldPolarization;
1197         // Here we Follow Jackson's conventio    1145         // Here we Follow Jackson's conventions and set the parallel
1198         // component = 1 in case of a ray per    1146         // component = 1 in case of a ray perpendicular to the surface
1199         E1_perp = 0.0;                           1147         E1_perp = 0.0;
1200         E1_parl = 1.0;                           1148         E1_parl = 1.0;
1201       }                                          1149       }
1202                                                  1150 
1203       s1       = fRindex1 * cost1;               1151       s1       = fRindex1 * cost1;
1204       E2_perp  = 2. * s1 * E1_perp / (fRindex    1152       E2_perp  = 2. * s1 * E1_perp / (fRindex1 * cost1 + fRindex2 * cost2);
1205       E2_parl  = 2. * s1 * E1_parl / (fRindex    1153       E2_parl  = 2. * s1 * E1_parl / (fRindex2 * cost1 + fRindex1 * cost2);
1206       E2_total = E2_perp * E2_perp + E2_parl     1154       E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1207       s2       = fRindex2 * cost2 * E2_total;    1155       s2       = fRindex2 * cost2 * E2_total;
1208                                                  1156 
1209       // D.Sawkey, 24 May 24                  << 1157       if(fTransmittance > 0.)
1210       // Transmittance has already been taken << 1158         transCoeff = fTransmittance;
1211       // For e.g. specular surfaces, the rati << 1159       else if(cost1 != 0.0)
1212       // reflection should be given by the ma << 
1213       // TRANSMITTANCE                        << 
1214       //if(fTransmittance > 0.)               << 
1215       //  transCoeff = fTransmittance;        << 
1216       //else if(cost1 != 0.0)                 << 
1217       if(cost1 != 0.0)                        << 
1218         transCoeff = s2 / s1;                    1160         transCoeff = s2 / s1;
1219       else                                       1161       else
1220         transCoeff = 0.0;                        1162         transCoeff = 0.0;
1221                                                  1163 
1222       // NOT TIR: REFLECTION                     1164       // NOT TIR: REFLECTION
1223       if(!G4BooleanRand(transCoeff))             1165       if(!G4BooleanRand(transCoeff))
1224       {                                          1166       {
1225         swap    = false;                         1167         swap    = false;
1226         fStatus = FresnelReflection;             1168         fStatus = FresnelReflection;
1227                                                  1169 
1228         if(!surfaceRoughnessCriterionPass)       1170         if(!surfaceRoughnessCriterionPass)
1229           fStatus = LambertianReflection;        1171           fStatus = LambertianReflection;
1230         if(fModel == unified && fFinish != po    1172         if(fModel == unified && fFinish != polished)
1231           ChooseReflection();                    1173           ChooseReflection();
1232         if(fStatus == LambertianReflection)      1174         if(fStatus == LambertianReflection)
1233         {                                        1175         {
1234           DoReflection();                        1176           DoReflection();
1235         }                                        1177         }
1236         else if(fStatus == BackScattering)       1178         else if(fStatus == BackScattering)
1237         {                                        1179         {
1238           fNewMomentum     = -fOldMomentum;      1180           fNewMomentum     = -fOldMomentum;
1239           fNewPolarization = -fOldPolarizatio    1181           fNewPolarization = -fOldPolarization;
1240         }                                        1182         }
1241         else                                     1183         else
1242         {                                        1184         {
1243           fNewMomentum =                         1185           fNewMomentum =
1244             fOldMomentum - 2. * fOldMomentum     1186             fOldMomentum - 2. * fOldMomentum * fFacetNormal * fFacetNormal;
1245           if(fSint1 > 0.0)                       1187           if(fSint1 > 0.0)
1246           {  // incident ray oblique             1188           {  // incident ray oblique
1247             E2_parl  = fRindex2 * E2_parl / f    1189             E2_parl  = fRindex2 * E2_parl / fRindex1 - E1_parl;
1248             E2_perp  = E2_perp - E1_perp;        1190             E2_perp  = E2_perp - E1_perp;
1249             E2_total = E2_perp * E2_perp + E2    1191             E2_total = E2_perp * E2_perp + E2_parl * E2_parl;
1250             A_paral  = (fNewMomentum.cross(A_    1192             A_paral  = (fNewMomentum.cross(A_trans)).unit();
1251             E2_abs   = std::sqrt(E2_total);      1193             E2_abs   = std::sqrt(E2_total);
1252             C_parl   = E2_parl / E2_abs;         1194             C_parl   = E2_parl / E2_abs;
1253             C_perp   = E2_perp / E2_abs;         1195             C_perp   = E2_perp / E2_abs;
1254                                                  1196 
1255             fNewPolarization = C_parl * A_par    1197             fNewPolarization = C_parl * A_paral + C_perp * A_trans;
1256           }                                      1198           }
1257           else                                   1199           else
1258           {  // incident ray perpendicular       1200           {  // incident ray perpendicular
1259             if(fRindex2 > fRindex1)              1201             if(fRindex2 > fRindex1)
1260             {                                    1202             {
1261               fNewPolarization = -fOldPolariz    1203               fNewPolarization = -fOldPolarization;
1262             }                                    1204             }
1263             else                                 1205             else
1264             {                                    1206             {
1265               fNewPolarization = fOldPolariza    1207               fNewPolarization = fOldPolarization;
1266             }                                    1208             }
1267           }                                      1209           }
1268         }                                        1210         }
1269       }                                          1211       }
1270       // NOT TIR: TRANSMISSION                   1212       // NOT TIR: TRANSMISSION
1271       else                                       1213       else
1272       {                                          1214       {
1273         inside  = !inside;                       1215         inside  = !inside;
1274         through = true;                          1216         through = true;
1275         fStatus = FresnelRefraction;             1217         fStatus = FresnelRefraction;
1276                                                  1218 
1277         if(fSint1 > 0.0)                         1219         if(fSint1 > 0.0)
1278         {  // incident ray oblique               1220         {  // incident ray oblique
1279           alpha        = cost1 - cost2 * (fRi    1221           alpha        = cost1 - cost2 * (fRindex2 / fRindex1);
1280           fNewMomentum = (fOldMomentum + alph    1222           fNewMomentum = (fOldMomentum + alpha * fFacetNormal).unit();
1281           A_paral      = (fNewMomentum.cross(    1223           A_paral      = (fNewMomentum.cross(A_trans)).unit();
1282           E2_abs       = std::sqrt(E2_total);    1224           E2_abs       = std::sqrt(E2_total);
1283           C_parl       = E2_parl / E2_abs;       1225           C_parl       = E2_parl / E2_abs;
1284           C_perp       = E2_perp / E2_abs;       1226           C_perp       = E2_perp / E2_abs;
1285                                                  1227 
1286           fNewPolarization = C_parl * A_paral    1228           fNewPolarization = C_parl * A_paral + C_perp * A_trans;
1287         }                                        1229         }
1288         else                                     1230         else
1289         {  // incident ray perpendicular         1231         {  // incident ray perpendicular
1290           fNewMomentum     = fOldMomentum;       1232           fNewMomentum     = fOldMomentum;
1291           fNewPolarization = fOldPolarization    1233           fNewPolarization = fOldPolarization;
1292         }                                        1234         }
1293       }                                          1235       }
1294     }                                            1236     }
1295                                                  1237 
1296     fOldMomentum     = fNewMomentum.unit();      1238     fOldMomentum     = fNewMomentum.unit();
1297     fOldPolarization = fNewPolarization.unit(    1239     fOldPolarization = fNewPolarization.unit();
1298                                                  1240 
1299     if(fStatus == FresnelRefraction)             1241     if(fStatus == FresnelRefraction)
1300     {                                            1242     {
1301       done = (fNewMomentum * fGlobalNormal <=    1243       done = (fNewMomentum * fGlobalNormal <= 0.0);
1302     }                                            1244     }
1303     else                                         1245     else
1304     {                                            1246     {
1305       done = (fNewMomentum * fGlobalNormal >=    1247       done = (fNewMomentum * fGlobalNormal >= -fCarTolerance);
1306     }                                            1248     }
1307     // Loop checking, 13-Aug-2015, Peter Gump    1249     // Loop checking, 13-Aug-2015, Peter Gumplinger
1308   } while(!done);                                1250   } while(!done);
1309                                                  1251 
1310   if(inside && !swap)                            1252   if(inside && !swap)
1311   {                                              1253   {
1312     if(fFinish == polishedbackpainted || fFin    1254     if(fFinish == polishedbackpainted || fFinish == groundbackpainted)
1313     {                                            1255     {
1314       G4double rand = G4UniformRand();           1256       G4double rand = G4UniformRand();
1315       if(rand > fReflectivity + fTransmittanc    1257       if(rand > fReflectivity + fTransmittance)
1316       {                                          1258       {
1317         DoAbsorption();                          1259         DoAbsorption();
1318       }                                          1260       }
1319       else if(rand > fReflectivity)              1261       else if(rand > fReflectivity)
1320       {                                          1262       {
1321         fStatus          = Transmission;         1263         fStatus          = Transmission;
1322         fNewMomentum     = fOldMomentum;         1264         fNewMomentum     = fOldMomentum;
1323         fNewPolarization = fOldPolarization;     1265         fNewPolarization = fOldPolarization;
1324       }                                          1266       }
1325       else                                       1267       else
1326       {                                          1268       {
1327         if(fStatus != FresnelRefraction)         1269         if(fStatus != FresnelRefraction)
1328         {                                        1270         {
1329           fGlobalNormal = -fGlobalNormal;        1271           fGlobalNormal = -fGlobalNormal;
1330         }                                        1272         }
1331         else                                     1273         else
1332         {                                        1274         {
1333           swap = !swap;                          1275           swap = !swap;
1334           G4SwapPtr(fMaterial1, fMaterial2);     1276           G4SwapPtr(fMaterial1, fMaterial2);
1335           G4SwapObj(&fRindex1, &fRindex2);       1277           G4SwapObj(&fRindex1, &fRindex2);
1336         }                                        1278         }
1337         if(fFinish == groundbackpainted)         1279         if(fFinish == groundbackpainted)
1338           fStatus = LambertianReflection;        1280           fStatus = LambertianReflection;
1339                                                  1281 
1340         DoReflection();                          1282         DoReflection();
1341                                                  1283 
1342         fGlobalNormal = -fGlobalNormal;          1284         fGlobalNormal = -fGlobalNormal;
1343         fOldMomentum  = fNewMomentum;            1285         fOldMomentum  = fNewMomentum;
1344                                                  1286 
1345         goto leap;                               1287         goto leap;
1346       }                                          1288       }
1347     }                                            1289     }
1348   }                                              1290   }
1349 }                                                1291 }
1350                                                  1292 
1351 //....oooOO0OOooo........oooOO0OOooo........o    1293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1352 G4double G4OpBoundaryProcess::GetMeanFreePath    1294 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track&, G4double,
1353                                                  1295                                               G4ForceCondition* condition)
1354 {                                                1296 {
1355   *condition = Forced;                           1297   *condition = Forced;
1356   return DBL_MAX;                                1298   return DBL_MAX;
1357 }                                                1299 }
1358                                                  1300 
1359 //....oooOO0OOooo........oooOO0OOooo........o    1301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1360 G4double G4OpBoundaryProcess::GetIncidentAngl    1302 G4double G4OpBoundaryProcess::GetIncidentAngle()
1361 {                                                1303 {
1362   return pi - std::acos(fOldMomentum * fFacet    1304   return pi - std::acos(fOldMomentum * fFacetNormal /
1363                         (fOldMomentum.mag() *    1305                         (fOldMomentum.mag() * fFacetNormal.mag()));
1364 }                                                1306 }
1365                                                  1307 
1366 //....oooOO0OOooo........oooOO0OOooo........o    1308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1367 G4double G4OpBoundaryProcess::GetReflectivity    1309 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1368                                                  1310                                               G4double E1_parl,
1369                                                  1311                                               G4double incidentangle,
1370                                                  1312                                               G4double realRindex,
1371                                                  1313                                               G4double imaginaryRindex)
1372 {                                                1314 {
1373   G4complex reflectivity, reflectivity_TE, re    1315   G4complex reflectivity, reflectivity_TE, reflectivity_TM;
1374   G4complex N1(fRindex1, 0.), N2(realRindex,     1316   G4complex N1(fRindex1, 0.), N2(realRindex, imaginaryRindex);
1375   G4complex cosPhi;                              1317   G4complex cosPhi;
1376                                                  1318 
1377   G4complex u(1., 0.);  // unit number 1         1319   G4complex u(1., 0.);  // unit number 1
1378                                                  1320 
1379   G4complex numeratorTE;  // E1_perp=1 E1_par    1321   G4complex numeratorTE;  // E1_perp=1 E1_parl=0 -> TE polarization
1380   G4complex numeratorTM;  // E1_parl=1 E1_per    1322   G4complex numeratorTM;  // E1_parl=1 E1_perp=0 -> TM polarization
1381   G4complex denominatorTE, denominatorTM;        1323   G4complex denominatorTE, denominatorTM;
1382   G4complex rTM, rTE;                            1324   G4complex rTM, rTE;
1383                                                  1325 
1384   G4MaterialPropertiesTable* MPT = fMaterial1    1326   G4MaterialPropertiesTable* MPT = fMaterial1->GetMaterialPropertiesTable();
1385   G4MaterialPropertyVector* ppR  = MPT->GetPr    1327   G4MaterialPropertyVector* ppR  = MPT->GetProperty(kREALRINDEX);
1386   G4MaterialPropertyVector* ppI  = MPT->GetPr    1328   G4MaterialPropertyVector* ppI  = MPT->GetProperty(kIMAGINARYRINDEX);
1387   if(ppR && ppI)                                 1329   if(ppR && ppI)
1388   {                                              1330   {
1389     G4double rRindex = ppR->Value(fPhotonMome    1331     G4double rRindex = ppR->Value(fPhotonMomentum, idx_rrindex);
1390     G4double iRindex = ppI->Value(fPhotonMome    1332     G4double iRindex = ppI->Value(fPhotonMomentum, idx_irindex);
1391     N1               = G4complex(rRindex, iRi    1333     N1               = G4complex(rRindex, iRindex);
1392   }                                              1334   }
1393                                                  1335 
1394   // Following two equations, rTM and rTE, ar    1336   // Following two equations, rTM and rTE, are from: "Introduction To Modern
1395   // Optics" written by Fowles                   1337   // Optics" written by Fowles
1396   cosPhi = std::sqrt(u - ((std::sin(incidenta    1338   cosPhi = std::sqrt(u - ((std::sin(incidentangle) * std::sin(incidentangle)) *
1397                           (N1 * N1) / (N2 * N    1339                           (N1 * N1) / (N2 * N2)));
1398                                                  1340 
1399   numeratorTE   = N1 * std::cos(incidentangle    1341   numeratorTE   = N1 * std::cos(incidentangle) - N2 * cosPhi;
1400   denominatorTE = N1 * std::cos(incidentangle    1342   denominatorTE = N1 * std::cos(incidentangle) + N2 * cosPhi;
1401   rTE           = numeratorTE / denominatorTE    1343   rTE           = numeratorTE / denominatorTE;
1402                                                  1344 
1403   numeratorTM   = N2 * std::cos(incidentangle    1345   numeratorTM   = N2 * std::cos(incidentangle) - N1 * cosPhi;
1404   denominatorTM = N2 * std::cos(incidentangle    1346   denominatorTM = N2 * std::cos(incidentangle) + N1 * cosPhi;
1405   rTM           = numeratorTM / denominatorTM    1347   rTM           = numeratorTM / denominatorTM;
1406                                                  1348 
1407   // This is my (PG) calculaton for reflectiv    1349   // This is my (PG) calculaton for reflectivity on a metallic surface
1408   // depending on the fraction of TE and TM p    1350   // depending on the fraction of TE and TM polarization
1409   // when TE polarization, E1_parl=0 and E1_p    1351   // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1410   // when TM polarization, E1_parl=1 and E1_p    1352   // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1411                                                  1353 
1412   reflectivity_TE = (rTE * conj(rTE)) * (E1_p    1354   reflectivity_TE = (rTE * conj(rTE)) * (E1_perp * E1_perp) /
1413                     (E1_perp * E1_perp + E1_p    1355                     (E1_perp * E1_perp + E1_parl * E1_parl);
1414   reflectivity_TM = (rTM * conj(rTM)) * (E1_p    1356   reflectivity_TM = (rTM * conj(rTM)) * (E1_parl * E1_parl) /
1415                     (E1_perp * E1_perp + E1_p    1357                     (E1_perp * E1_perp + E1_parl * E1_parl);
1416   reflectivity = reflectivity_TE + reflectivi    1358   reflectivity = reflectivity_TE + reflectivity_TM;
1417                                                  1359 
1418   do                                             1360   do
1419   {                                              1361   {
1420     if(G4UniformRand() * real(reflectivity) >    1362     if(G4UniformRand() * real(reflectivity) > real(reflectivity_TE))
1421     {                                            1363     {
1422       f_iTE = -1;                                1364       f_iTE = -1;
1423     }                                            1365     }
1424     else                                         1366     else
1425     {                                            1367     {
1426       f_iTE = 1;                                 1368       f_iTE = 1;
1427     }                                            1369     }
1428     if(G4UniformRand() * real(reflectivity) >    1370     if(G4UniformRand() * real(reflectivity) > real(reflectivity_TM))
1429     {                                            1371     {
1430       f_iTM = -1;                                1372       f_iTM = -1;
1431     }                                            1373     }
1432     else                                         1374     else
1433     {                                            1375     {
1434       f_iTM = 1;                                 1376       f_iTM = 1;
1435     }                                            1377     }
1436     // Loop checking, 13-Aug-2015, Peter Gump    1378     // Loop checking, 13-Aug-2015, Peter Gumplinger
1437   } while(f_iTE < 0 && f_iTM < 0);               1379   } while(f_iTE < 0 && f_iTM < 0);
1438                                                  1380 
1439   return real(reflectivity);                     1381   return real(reflectivity);
1440 }                                                1382 }
1441                                                  1383 
1442 //....oooOO0OOooo........oooOO0OOooo........o    1384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 1385 
1443 void G4OpBoundaryProcess::CalculateReflectivi    1386 void G4OpBoundaryProcess::CalculateReflectivity()
1444 {                                                1387 {
1445   G4double realRindex = fRealRIndexMPV->Value    1388   G4double realRindex = fRealRIndexMPV->Value(fPhotonMomentum, idx_rrindex);
1446   G4double imaginaryRindex =                     1389   G4double imaginaryRindex =
1447     fImagRIndexMPV->Value(fPhotonMomentum, id    1390     fImagRIndexMPV->Value(fPhotonMomentum, idx_irindex);
1448                                                  1391 
1449   // calculate FacetNormal                       1392   // calculate FacetNormal
1450   if(fFinish == ground)                          1393   if(fFinish == ground)
1451   {                                              1394   {
1452     fFacetNormal = GetFacetNormal(fOldMomentu    1395     fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
1453   }                                              1396   }
1454   else                                           1397   else
1455   {                                              1398   {
1456     fFacetNormal = fGlobalNormal;                1399     fFacetNormal = fGlobalNormal;
1457   }                                              1400   }
1458                                                  1401 
1459   G4double cost1 = -fOldMomentum * fFacetNorm    1402   G4double cost1 = -fOldMomentum * fFacetNormal;
1460   if(std::abs(cost1) < 1.0 - fCarTolerance)      1403   if(std::abs(cost1) < 1.0 - fCarTolerance)
1461   {                                              1404   {
1462     fSint1 = std::sqrt(1. - cost1 * cost1);      1405     fSint1 = std::sqrt(1. - cost1 * cost1);
1463   }                                              1406   }
1464   else                                           1407   else
1465   {                                              1408   {
1466     fSint1 = 0.0;                                1409     fSint1 = 0.0;
1467   }                                              1410   }
1468                                                  1411 
1469   G4ThreeVector A_trans, A_paral, E1pp, E1pl;    1412   G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1470   G4double E1_perp, E1_parl;                     1413   G4double E1_perp, E1_parl;
1471                                                  1414 
1472   if(fSint1 > 0.0)                               1415   if(fSint1 > 0.0)
1473   {                                              1416   {
1474     A_trans = (fOldMomentum.cross(fFacetNorma    1417     A_trans = (fOldMomentum.cross(fFacetNormal)).unit();
1475     E1_perp = fOldPolarization * A_trans;        1418     E1_perp = fOldPolarization * A_trans;
1476     E1pp    = E1_perp * A_trans;                 1419     E1pp    = E1_perp * A_trans;
1477     E1pl    = fOldPolarization - E1pp;           1420     E1pl    = fOldPolarization - E1pp;
1478     E1_parl = E1pl.mag();                        1421     E1_parl = E1pl.mag();
1479   }                                              1422   }
1480   else                                           1423   else
1481   {                                              1424   {
1482     A_trans = fOldPolarization;                  1425     A_trans = fOldPolarization;
1483     // Here we Follow Jackson's conventions a    1426     // Here we Follow Jackson's conventions and we set the parallel
1484     // component = 1 in case of a ray perpend    1427     // component = 1 in case of a ray perpendicular to the surface
1485     E1_perp = 0.0;                               1428     E1_perp = 0.0;
1486     E1_parl = 1.0;                               1429     E1_parl = 1.0;
1487   }                                              1430   }
1488                                                  1431 
1489   G4double incidentangle = GetIncidentAngle()    1432   G4double incidentangle = GetIncidentAngle();
1490                                                  1433 
1491   // calculate the reflectivity depending on     1434   // calculate the reflectivity depending on incident angle,
1492   // polarization and complex refractive         1435   // polarization and complex refractive
1493   fReflectivity = GetReflectivity(E1_perp, E1    1436   fReflectivity = GetReflectivity(E1_perp, E1_parl, incidentangle, realRindex,
1494                                   imaginaryRi    1437                                   imaginaryRindex);
1495 }                                                1438 }
1496                                                  1439 
1497 //....oooOO0OOooo........oooOO0OOooo........o    1440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1498 G4bool G4OpBoundaryProcess::InvokeSD(const G4    1441 G4bool G4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
1499 {                                                1442 {
1500   G4Step aStep = *pStep;                         1443   G4Step aStep = *pStep;
1501   aStep.AddTotalEnergyDeposit(fPhotonMomentum    1444   aStep.AddTotalEnergyDeposit(fPhotonMomentum);
1502                                                  1445 
1503   G4VSensitiveDetector* sd = aStep.GetPostSte    1446   G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
1504   if(sd != nullptr)                              1447   if(sd != nullptr)
1505     return sd->Hit(&aStep);                      1448     return sd->Hit(&aStep);
1506   else                                           1449   else
1507     return false;                                1450     return false;
1508 }                                                1451 }
1509                                                  1452 
1510 //....oooOO0OOooo........oooOO0OOooo........o    1453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1511 inline void G4OpBoundaryProcess::SetInvokeSD(    1454 inline void G4OpBoundaryProcess::SetInvokeSD(G4bool flag)
1512 {                                                1455 {
1513   fInvokeSD = flag;                              1456   fInvokeSD = flag;
1514   G4OpticalParameters::Instance()->SetBoundar    1457   G4OpticalParameters::Instance()->SetBoundaryInvokeSD(fInvokeSD);
1515 }                                                1458 }
1516                                                  1459 
1517 //....oooOO0OOooo........oooOO0OOooo........o    1460 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1518 void G4OpBoundaryProcess::SetVerboseLevel(G4i    1461 void G4OpBoundaryProcess::SetVerboseLevel(G4int verbose)
1519 {                                                1462 {
1520   verboseLevel = verbose;                        1463   verboseLevel = verbose;
1521   G4OpticalParameters::Instance()->SetBoundar    1464   G4OpticalParameters::Instance()->SetBoundaryVerboseLevel(verboseLevel);
1522 }                                             << 
1523                                               << 
1524 //....oooOO0OOooo........oooOO0OOooo........o << 
1525 void G4OpBoundaryProcess::CoatedDielectricDie << 
1526 {                                             << 
1527   G4MaterialPropertyVector* pp = nullptr;     << 
1528                                               << 
1529   G4MaterialPropertiesTable* MPT = fMaterial2 << 
1530   if((pp = MPT->GetProperty(kRINDEX)))        << 
1531   {                                           << 
1532     fRindex2 = pp->Value(fPhotonMomentum, idx << 
1533   }                                           << 
1534                                               << 
1535   MPT = fOpticalSurface->GetMaterialPropertie << 
1536   if((pp = MPT->GetProperty(kCOATEDRINDEX)))  << 
1537   {                                           << 
1538     fCoatedRindex = pp->Value(fPhotonMomentum << 
1539   }                                           << 
1540   if(MPT->ConstPropertyExists(kCOATEDTHICKNES << 
1541   {                                           << 
1542     fCoatedThickness = MPT->GetConstProperty( << 
1543   }                                           << 
1544   if(MPT->ConstPropertyExists(kCOATEDFRUSTRAT << 
1545   {                                           << 
1546     fCoatedFrustratedTransmission =           << 
1547       (G4bool)MPT->GetConstProperty(kCOATEDFR << 
1548   }                                           << 
1549                                               << 
1550   G4double sintTL;                            << 
1551   G4double wavelength = h_Planck * c_light /  << 
1552   G4double PdotN;                             << 
1553   G4double E1_perp, E1_parl;                  << 
1554   G4double s1, E2_perp, E2_parl, E2_total, tr << 
1555   G4double E2_abs, C_parl, C_perp;            << 
1556   G4double alpha;                             << 
1557   G4ThreeVector A_trans, A_paral, E1pp, E1pl; << 
1558   //G4bool Inside  = false;                   << 
1559   //G4bool Swap    = false;                   << 
1560   G4bool through = false;                     << 
1561   G4bool done    = false;                     << 
1562                                               << 
1563   do {                                        << 
1564     if (through)                              << 
1565     {                                         << 
1566       //Swap = !Swap;                         << 
1567       through = false;                        << 
1568       fGlobalNormal = -fGlobalNormal;         << 
1569       G4SwapPtr(fMaterial1, fMaterial2);      << 
1570       G4SwapObj(&fRindex1, &fRindex2);        << 
1571     }                                         << 
1572                                               << 
1573     if(fFinish == polished)                   << 
1574     {                                         << 
1575       fFacetNormal = fGlobalNormal;           << 
1576     }                                         << 
1577     else                                      << 
1578     {                                         << 
1579       fFacetNormal = GetFacetNormal(fOldMomen << 
1580     }                                         << 
1581                                               << 
1582     PdotN = fOldMomentum * fFacetNormal;      << 
1583     G4double cost1 = -PdotN;                  << 
1584     G4double sint2, cost2 = 0.;               << 
1585                                               << 
1586     if (std::abs(cost1) < 1.0 - fCarTolerance << 
1587     {                                         << 
1588       fSint1 = std::sqrt(1. - cost1 * cost1); << 
1589       sint2 = fSint1 * fRindex1 / fRindex2;   << 
1590       sintTL = fSint1 * fRindex1 / fCoatedRin << 
1591     } else                                    << 
1592     {                                         << 
1593       fSint1 = 0.0;                           << 
1594       sint2 = 0.0;                            << 
1595       sintTL = 0.0;                           << 
1596     }                                         << 
1597                                               << 
1598     if (fSint1 > 0.0)                         << 
1599     {                                         << 
1600       A_trans = fOldMomentum.cross(fFacetNorm << 
1601       A_trans = A_trans.unit();               << 
1602       E1_perp = fOldPolarization * A_trans;   << 
1603       E1pp = E1_perp * A_trans;               << 
1604       E1pl = fOldPolarization - E1pp;         << 
1605       E1_parl = E1pl.mag();                   << 
1606     }                                         << 
1607     else                                      << 
1608     {                                         << 
1609       A_trans = fOldPolarization;             << 
1610       E1_perp = 0.0;                          << 
1611       E1_parl = 1.0;                          << 
1612     }                                         << 
1613                                               << 
1614     s1 = fRindex1 * cost1;                    << 
1615                                               << 
1616     if (cost1 > 0.0)                          << 
1617     {                                         << 
1618       cost2 = std::sqrt(1. - sint2 * sint2);  << 
1619     }                                         << 
1620     else                                      << 
1621     {                                         << 
1622       cost2 = -std::sqrt(1. - sint2 * sint2); << 
1623     }                                         << 
1624                                               << 
1625     transCoeff = 0.0;                         << 
1626                                               << 
1627     if (sintTL >= 1.0)                        << 
1628     { // --> Angle > Angle Limit              << 
1629       //Swap = false;                         << 
1630     }                                         << 
1631     E2_perp = 2. * s1 * E1_perp / (fRindex1 * << 
1632     E2_parl = 2. * s1 * E1_parl / (fRindex2 * << 
1633     E2_total = E2_perp * E2_perp + E2_parl *  << 
1634                                               << 
1635     transCoeff = 1. - GetReflectivityThroughT << 
1636                         sintTL, E1_perp, E1_p << 
1637     if (!G4BooleanRand(transCoeff))           << 
1638     {                                         << 
1639       if(verboseLevel > 2)                    << 
1640         G4cout << "Reflection from " << fMate << 
1641                << fMaterial2->GetName() << G4 << 
1642                                               << 
1643       //Swap = false;                         << 
1644                                               << 
1645       if (sintTL >= 1.0)                      << 
1646       {                                       << 
1647         fStatus = TotalInternalReflection;    << 
1648       }                                       << 
1649       else                                    << 
1650       {                                       << 
1651         fStatus = CoatedDielectricReflection; << 
1652       }                                       << 
1653                                               << 
1654       PdotN = fOldMomentum * fFacetNormal;    << 
1655       fNewMomentum = fOldMomentum - (2. * Pdo << 
1656                                               << 
1657       if (fSint1 > 0.0) {   // incident ray o << 
1658                                               << 
1659         E2_parl = fRindex2 * E2_parl / fRinde << 
1660         E2_perp = E2_perp - E1_perp;          << 
1661         E2_total = E2_perp * E2_perp + E2_par << 
1662         A_paral = fNewMomentum.cross(A_trans) << 
1663         A_paral = A_paral.unit();             << 
1664         E2_abs = std::sqrt(E2_total);         << 
1665         C_parl = E2_parl / E2_abs;            << 
1666         C_perp = E2_perp / E2_abs;            << 
1667                                               << 
1668         fNewPolarization = C_parl * A_paral + << 
1669                                               << 
1670       }                                       << 
1671       else                                    << 
1672       {               // incident ray perpend << 
1673         if (fRindex2 > fRindex1)              << 
1674         {                                     << 
1675           fNewPolarization = -fOldPolarizatio << 
1676         }                                     << 
1677         else                                  << 
1678         {                                     << 
1679           fNewPolarization = fOldPolarization << 
1680         }                                     << 
1681       }                                       << 
1682                                               << 
1683     } else { // photon gets transmitted       << 
1684       if (verboseLevel > 2)                   << 
1685         G4cout << "Transmission from " << fMa << 
1686                << fMaterial2->GetName() << G4 << 
1687                                               << 
1688       //Inside = !Inside;                     << 
1689       through = true;                         << 
1690                                               << 
1691       if (fEfficiency > 0.)                   << 
1692       {                                       << 
1693         DoAbsorption();                       << 
1694         return;                               << 
1695       }                                       << 
1696       else                                    << 
1697       {                                       << 
1698         if (sintTL >= 1.0)                    << 
1699         {                                     << 
1700           fStatus = CoatedDielectricFrustrate << 
1701         }                                     << 
1702         else                                  << 
1703         {                                     << 
1704           fStatus = CoatedDielectricRefractio << 
1705         }                                     << 
1706                                               << 
1707         if (fSint1 > 0.0) {      // incident  << 
1708                                               << 
1709           alpha = cost1 - cost2 * (fRindex2 / << 
1710           fNewMomentum = fOldMomentum + alpha << 
1711           fNewMomentum = fNewMomentum.unit(); << 
1712           A_paral = fNewMomentum.cross(A_tran << 
1713           A_paral = A_paral.unit();           << 
1714           E2_abs = std::sqrt(E2_total);       << 
1715           C_parl = E2_parl / E2_abs;          << 
1716           C_perp = E2_perp / E2_abs;          << 
1717                                               << 
1718           fNewPolarization = C_parl * A_paral << 
1719                                               << 
1720         }                                     << 
1721         else                                  << 
1722         {                  // incident ray pe << 
1723           fNewMomentum = fOldMomentum;        << 
1724           fNewPolarization = fOldPolarization << 
1725         }                                     << 
1726       }                                       << 
1727     }                                         << 
1728                                               << 
1729     fOldMomentum = fNewMomentum.unit();       << 
1730     fOldPolarization = fNewPolarization.unit( << 
1731     if ((fStatus == CoatedDielectricFrustrate << 
1732         (fStatus == CoatedDielectricRefractio << 
1733     {                                         << 
1734       done = (fNewMomentum * fGlobalNormal <= << 
1735     }                                         << 
1736     else                                      << 
1737     {                                         << 
1738       done = (fNewMomentum * fGlobalNormal >= << 
1739     }                                         << 
1740                                               << 
1741   } while (!done);                            << 
1742 }                                             << 
1743                                               << 
1744 //....oooOO0OOooo........oooOO0OOooo........o << 
1745 G4double G4OpBoundaryProcess::GetReflectivity << 
1746                    G4double E1_perp,          << 
1747                    G4double E1_parl,          << 
1748                    G4double wavelength, G4dou << 
1749   G4complex Reflectivity, Reflectivity_TE, Re << 
1750   G4double gammaTL, costTL;                   << 
1751                                               << 
1752   G4complex i(0, 1);                          << 
1753   G4complex rTM, rTE;                         << 
1754   G4complex r1toTL, rTLto2;                   << 
1755   G4double k0 = 2 * pi / wavelength;          << 
1756                                               << 
1757   // Angle > Angle limit                      << 
1758   if (sinTL >= 1.0) {                         << 
1759     if (fCoatedFrustratedTransmission) { //Fr << 
1760                                               << 
1761       if (cost1 > 0.0)                        << 
1762       {                                       << 
1763         gammaTL = std::sqrt(fRindex1 * fRinde << 
1764                    fCoatedRindex * fCoatedRin << 
1765       }                                       << 
1766       else                                    << 
1767       {                                       << 
1768         gammaTL = -std::sqrt(fRindex1 * fRind << 
1769                    fCoatedRindex * fCoatedRin << 
1770       }                                       << 
1771                                               << 
1772       // TE                                   << 
1773       r1toTL = (fRindex1 * cost1 - i * gammaT << 
1774       rTLto2 = (i * gammaTL - fRindex2 * cost << 
1775       if (cost1 != 0.0)                       << 
1776       {                                       << 
1777         rTE = (r1toTL + rTLto2 * std::exp(-2  << 
1778                  (1.0 + r1toTL * rTLto2 * std << 
1779       }                                       << 
1780       // TM                                   << 
1781       r1toTL = (fRindex1 * i * gammaTL - fCoa << 
1782                   (fRindex1 * i * gammaTL + f << 
1783       rTLto2 = (fCoatedRindex * fCoatedRindex << 
1784                   (fCoatedRindex * fCoatedRin << 
1785       if (cost1 != 0.0)                       << 
1786       {                                       << 
1787         rTM = (r1toTL + rTLto2 * std::exp(-2  << 
1788                  (1.0 + r1toTL * rTLto2 * std << 
1789       }                                       << 
1790     }                                         << 
1791     else                                      << 
1792     { //Total reflection                      << 
1793       return(1.);                             << 
1794     }                                         << 
1795   }                                           << 
1796                                               << 
1797   // Angle <= Angle limit                     << 
1798   else //if (sinTL < 1.0)                     << 
1799   {                                           << 
1800     if (cost1 > 0.0)                          << 
1801     {                                         << 
1802       costTL = std::sqrt(1. - sinTL * sinTL); << 
1803     }                                         << 
1804     else                                      << 
1805     {                                         << 
1806       costTL = -std::sqrt(1. - sinTL * sinTL) << 
1807     }                                         << 
1808     // TE                                     << 
1809     r1toTL = (fRindex1 * cost1 - fCoatedRinde << 
1810     rTLto2 = (fCoatedRindex * costTL - fRinde << 
1811     if (cost1 != 0.0)                         << 
1812     {                                         << 
1813       rTE = (r1toTL + rTLto2 * std::exp(2.0 * << 
1814             (1.0 + r1toTL * rTLto2 * std::exp << 
1815     }                                         << 
1816     // TM                                     << 
1817     r1toTL = (fRindex1 * costTL - fCoatedRind << 
1818     rTLto2 = (fCoatedRindex * cost2 - fRindex << 
1819     if (cost1 != 0.0)                         << 
1820     {                                         << 
1821       rTM = (r1toTL + rTLto2 * std::exp(2.0 * << 
1822             (1.0 + r1toTL * rTLto2 * std::exp << 
1823     }                                         << 
1824   }                                           << 
1825                                               << 
1826   Reflectivity_TE = (rTE * conj(rTE)) * (E1_p << 
1827   Reflectivity_TM = (rTM * conj(rTM)) * (E1_p << 
1828   Reflectivity = Reflectivity_TE + Reflectivi << 
1829                                               << 
1830   return real(Reflectivity);                  << 
1831 }                                                1465 }
1832                                                  1466