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 9.5.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 << 
 67 //                           of a dichronic fi << 
 68 //              2017-02-24 - add capability of << 
 69 //                           with Look-Up-Tabl << 
 70 //                                                 66 //
 71 // Author:      Peter Gumplinger                   67 // Author:      Peter Gumplinger
 72 //    adopted from work by Werner Keil - April     68 //    adopted from work by Werner Keil - April 2/96
                                                   >>  69 // mail:        gum@triumf.ca
 73 //                                                 70 //
 74 //////////////////////////////////////////////     71 ////////////////////////////////////////////////////////////////////////
 75                                                    72 
 76 #include "G4OpBoundaryProcess.hh"              << 
 77                                                << 
 78 #include "G4ios.hh"                                73 #include "G4ios.hh"
 79 #include "G4GeometryTolerance.hh"              << 
 80 #include "G4LogicalBorderSurface.hh"           << 
 81 #include "G4LogicalSkinSurface.hh"             << 
 82 #include "G4OpProcessSubType.hh"                   74 #include "G4OpProcessSubType.hh"
 83 #include "G4OpticalParameters.hh"              << 
 84 #include "G4ParallelWorldProcess.hh"           << 
 85 #include "G4PhysicalConstants.hh"              << 
 86 #include "G4SystemOfUnits.hh"                  << 
 87 #include "G4TransportationManager.hh"          << 
 88 #include "G4VSensitiveDetector.hh"             << 
 89                                                    75 
 90 //....oooOO0OOooo........oooOO0OOooo........oo <<  76 #include "G4OpBoundaryProcess.hh"
                                                   >>  77 #include "G4GeometryTolerance.hh"
                                                   >>  78 
                                                   >>  79 /////////////////////////
                                                   >>  80 // Class Implementation
                                                   >>  81 /////////////////////////
                                                   >>  82 
                                                   >>  83         //////////////
                                                   >>  84         // Operators
                                                   >>  85         //////////////
                                                   >>  86 
                                                   >>  87 // G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
                                                   >>  88 // {
                                                   >>  89 // }
                                                   >>  90 
                                                   >>  91         /////////////////
                                                   >>  92         // Constructors
                                                   >>  93         /////////////////
                                                   >>  94 
 91 G4OpBoundaryProcess::G4OpBoundaryProcess(const     95 G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
 92                                          G4Pro <<  96                                                G4ProcessType type)
 93   : G4VDiscreteProcess(processName, ptype)     <<  97              : G4VDiscreteProcess(processName, type)
 94 {                                                  98 {
 95   Initialise();                                <<  99         if ( verboseLevel > 0) {
                                                   >> 100            G4cout << GetProcessName() << " is created " << G4endl;
                                                   >> 101         }
 96                                                   102 
 97   if(verboseLevel > 0)                         << 103         SetProcessSubType(fOpBoundary);
 98   {                                            << 
 99     G4cout << GetProcessName() << " is created << 
100   }                                            << 
101   SetProcessSubType(fOpBoundary);              << 
102                                                   104 
103   fStatus           = Undefined;               << 105   theStatus = Undefined;
104   fModel            = glisur;                  << 106   theModel = glisur;
105   fFinish           = polished;                << 107   theFinish = polished;
106   fReflectivity     = 1.;                      << 108         theReflectivity =  1.;
107   fEfficiency       = 0.;                      << 109         theEfficiency   =  0.;
108   fTransmittance    = 0.;                      << 110         theTransmittance = 0.;
109   fSurfaceRoughness = 0.;                      << 
110   fProb_sl          = 0.;                      << 
111   fProb_ss          = 0.;                      << 
112   fProb_bs          = 0.;                      << 
113                                                << 
114   fRealRIndexMPV  = nullptr;                   << 
115   fImagRIndexMPV  = nullptr;                   << 
116   fMaterial1      = nullptr;                   << 
117   fMaterial2      = nullptr;                   << 
118   fOpticalSurface = nullptr;                   << 
119   fCarTolerance   = G4GeometryTolerance::GetIn << 
120                                                << 
121   f_iTE = f_iTM   = 0;                         << 
122   fPhotonMomentum = 0.;                        << 
123   fRindex1 = fRindex2 = 1.;                    << 
124   fSint1              = 0.;                    << 
125   fDichroicVector     = nullptr;               << 
126 }                                              << 
127                                                   111 
128 //....oooOO0OOooo........oooOO0OOooo........oo << 112         prob_sl = 0.;
129 G4OpBoundaryProcess::~G4OpBoundaryProcess() =  << 113         prob_ss = 0.;
                                                   >> 114         prob_bs = 0.;
130                                                   115 
131 //....oooOO0OOooo........oooOO0OOooo........oo << 116         PropertyPointer  = NULL;
132 void G4OpBoundaryProcess::PreparePhysicsTable( << 117         PropertyPointer1 = NULL;
133 {                                              << 118         PropertyPointer2 = NULL;
134   Initialise();                                << 
135 }                                              << 
136                                                   119 
137 //....oooOO0OOooo........oooOO0OOooo........oo << 120         kCarTolerance = G4GeometryTolerance::GetInstance()
138 void G4OpBoundaryProcess::Initialise()         << 121                         ->GetSurfaceTolerance();
139 {                                              << 122 
140   G4OpticalParameters* params = G4OpticalParam << 123         iTE = iTM = 0;
141   SetInvokeSD(params->GetBoundaryInvokeSD());  << 124         thePhotonMomentum = 0.;
142   SetVerboseLevel(params->GetBoundaryVerboseLe << 125         Rindex1 = Rindex2 = cost1 = cost2 = sint1 = sint2 = 0.;
143 }                                                 126 }
144                                                   127 
145 //....oooOO0OOooo........oooOO0OOooo........oo << 128 // G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
146 G4VParticleChange* G4OpBoundaryProcess::PostSt << 129 // {
147                                                << 130 // }
                                                   >> 131 
                                                   >> 132         ////////////////
                                                   >> 133         // Destructors
                                                   >> 134         ////////////////
                                                   >> 135 
                                                   >> 136 G4OpBoundaryProcess::~G4OpBoundaryProcess(){}
                                                   >> 137 
                                                   >> 138         ////////////
                                                   >> 139         // Methods
                                                   >> 140         ////////////
                                                   >> 141 
                                                   >> 142 // PostStepDoIt
                                                   >> 143 // ------------
                                                   >> 144 //
                                                   >> 145 G4VParticleChange*
                                                   >> 146 G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
148 {                                                 147 {
149   fStatus = Undefined;                         << 148         theStatus = Undefined;
150   aParticleChange.Initialize(aTrack);          << 
151   aParticleChange.ProposeVelocity(aTrack.GetVe << 
152                                                << 
153   // Get hyperStep from  G4ParallelWorldProces << 
154   //  NOTE: PostSetpDoIt of this process to be << 
155   //  G4ParallelWorldProcess!                  << 
156   const G4Step* pStep = &aStep;                << 
157   const G4Step* hStep = G4ParallelWorldProcess << 
158   if(hStep != nullptr)                         << 
159     pStep = hStep;                             << 
160                                                << 
161   if(pStep->GetPostStepPoint()->GetStepStatus( << 
162   {                                            << 
163     fMaterial1 = pStep->GetPreStepPoint()->Get << 
164     fMaterial2 = pStep->GetPostStepPoint()->Ge << 
165   }                                            << 
166   else                                         << 
167   {                                            << 
168     fStatus = NotAtBoundary;                   << 
169     if(verboseLevel > 1)                       << 
170       BoundaryProcessVerbose();                << 
171     return G4VDiscreteProcess::PostStepDoIt(aT << 
172   }                                            << 
173                                                   149 
174   G4VPhysicalVolume* thePrePV  = pStep->GetPre << 150         aParticleChange.Initialize(aTrack);
175   G4VPhysicalVolume* thePostPV = pStep->GetPos << 151         aParticleChange.ProposeVelocity(aTrack.GetVelocity());
176                                                   152 
177   if(verboseLevel > 1)                         << 153         G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
178   {                                            << 154         G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
179     G4cout << " Photon at Boundary! " << G4end << 
180     if(thePrePV != nullptr)                    << 
181       G4cout << " thePrePV:  " << thePrePV->Ge << 
182     if(thePostPV != nullptr)                   << 
183       G4cout << " thePostPV: " << thePostPV->G << 
184   }                                            << 
185                                                   155 
186   G4double stepLength = aTrack.GetStepLength() << 156         if ( verboseLevel > 0 ) {
187   if(stepLength <= fCarTolerance)              << 157            G4cout << " Photon at Boundary! " << G4endl;
188   {                                            << 158            G4VPhysicalVolume* thePrePV = pPreStepPoint->GetPhysicalVolume();
189     fStatus = StepTooSmall;                    << 159            G4VPhysicalVolume* thePostPV = pPostStepPoint->GetPhysicalVolume();
190     if(verboseLevel > 1)                       << 160            if (thePrePV)  G4cout << " thePrePV:  " << thePrePV->GetName()  << G4endl;
191       BoundaryProcessVerbose();                << 161            if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
192                                                << 162         }
193     G4MaterialPropertyVector* groupvel = nullp << 
194     G4MaterialPropertiesTable* aMPT = fMateria << 
195     if(aMPT != nullptr)                        << 
196     {                                          << 
197       groupvel = aMPT->GetProperty(kGROUPVEL); << 
198     }                                          << 
199                                                << 
200     if(groupvel != nullptr)                    << 
201     {                                          << 
202       aParticleChange.ProposeVelocity(         << 
203         groupvel->Value(fPhotonMomentum, idx_g << 
204     }                                          << 
205     return G4VDiscreteProcess::PostStepDoIt(aT << 
206   }                                            << 
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                                                   163 
226   const G4DynamicParticle* aParticle = aTrack. << 164         if (pPostStepPoint->GetStepStatus() != fGeomBoundary){
                                                   >> 165           theStatus = NotAtBoundary;
                                                   >> 166                 if ( verboseLevel > 0) BoundaryProcessVerbose();
                                                   >> 167           return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 168   }
                                                   >> 169   if (aTrack.GetStepLength()<=kCarTolerance/2){
                                                   >> 170           theStatus = StepTooSmall;
                                                   >> 171                 if ( verboseLevel > 0) BoundaryProcessVerbose();
                                                   >> 172           return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 173   }
227                                                   174 
228   fPhotonMomentum  = aParticle->GetTotalMoment << 175   Material1 = pPreStepPoint  -> GetMaterial();
229   fOldMomentum     = aParticle->GetMomentumDir << 176   Material2 = pPostStepPoint -> GetMaterial();
230   fOldPolarization = aParticle->GetPolarizatio << 
231                                                << 
232   if(verboseLevel > 1)                         << 
233   {                                            << 
234     G4cout << " Old Momentum Direction: " << f << 
235            << " Old Polarization:       " << f << 
236   }                                            << 
237                                                   177 
238   G4ThreeVector theGlobalPoint = pStep->GetPos << 178         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
239   G4bool valid;                                << 
240                                                   179 
241   // ID of Navigator which limits step         << 180   thePhotonMomentum = aParticle->GetTotalMomentum();
242   G4int hNavId = G4ParallelWorldProcess::GetHy << 181         OldMomentum       = aParticle->GetMomentumDirection();
243   auto iNav    = G4TransportationManager::GetT << 182   OldPolarization   = aParticle->GetPolarization();
244                 ->GetActiveNavigatorsIterator( << 
245   fGlobalNormal = (iNav[hNavId])->GetGlobalExi << 
246                                                << 
247   if(valid)                                    << 
248   {                                            << 
249     fGlobalNormal = -fGlobalNormal;            << 
250   }                                            << 
251   else                                         << 
252   {                                            << 
253     G4ExceptionDescription ed;                 << 
254     ed << " G4OpBoundaryProcess/PostStepDoIt() << 
255        << " The Navigator reports that it retu << 
256     G4Exception(                               << 
257       "G4OpBoundaryProcess::PostStepDoIt", "Op << 
258       "Invalid Surface Normal - Geometry must  << 
259   }                                            << 
260                                                   183 
261   if(fOldMomentum * fGlobalNormal > 0.0)       << 184         if ( verboseLevel > 0 ) {
262   {                                            << 185            G4cout << " Old Momentum Direction: " << OldMomentum     << G4endl;
263 #ifdef G4OPTICAL_DEBUG                         << 186            G4cout << " Old Polarization:       " << OldPolarization << G4endl;
264     G4ExceptionDescription ed;                 << 187         }
265     ed << " G4OpBoundaryProcess/PostStepDoIt() << 
266           "wrong direction. "                  << 
267        << G4endl                               << 
268        << "   The momentum of the photon arriv << 
269        << "   must exit the volume cross in th << 
270        << "   So it MUST have dot < 0 with the << 
271           "volume (globalNormal)."             << 
272        << G4endl << "   >> The dot product of  << 
273        << fOldMomentum * fGlobalNormal << G4en << 
274        << "     Old Momentum  (during step)    << 
275        << "     Global Normal (Exiting New Vol << 
276        << G4endl;                              << 
277     G4Exception("G4OpBoundaryProcess::PostStep << 
278                 EventMustBeAborted,  // Or Jus << 
279                                      // repeat << 
280                 ed,                            << 
281                 "Invalid Surface Normal - Geom << 
282                 "normal pointing in the right  << 
283 #else                                          << 
284     fGlobalNormal = -fGlobalNormal;            << 
285 #endif                                         << 
286   }                                            << 
287                                                   188 
288   G4MaterialPropertyVector* rIndexMPV = nullpt << 189         G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
289   G4MaterialPropertiesTable* MPT = fMaterial1- << 
290   if(MPT != nullptr)                           << 
291   {                                            << 
292     rIndexMPV = MPT->GetProperty(kRINDEX);     << 
293   }                                            << 
294   if(rIndexMPV != nullptr)                     << 
295   {                                            << 
296     fRindex1 = rIndexMPV->Value(fPhotonMomentu << 
297   }                                            << 
298   else                                         << 
299   {                                            << 
300     fStatus = NoRINDEX;                        << 
301     if(verboseLevel > 1)                       << 
302       BoundaryProcessVerbose();                << 
303     aParticleChange.ProposeLocalEnergyDeposit( << 
304     aParticleChange.ProposeTrackStatus(fStopAn << 
305     return G4VDiscreteProcess::PostStepDoIt(aT << 
306   }                                            << 
307                                                   190 
308   fReflectivity      = 1.;                     << 191         G4Navigator* theNavigator =
309   fEfficiency        = 0.;                     << 192                      G4TransportationManager::GetTransportationManager()->
310   fTransmittance     = 0.;                     << 193                                               GetNavigatorForTracking();
311   fSurfaceRoughness  = 0.;                     << 
312   fModel             = glisur;                 << 
313   fFinish            = polished;               << 
314   G4SurfaceType type = dielectric_dielectric;  << 
315                                                << 
316   rIndexMPV       = nullptr;                   << 
317   fOpticalSurface = nullptr;                   << 
318                                                << 
319   G4LogicalSurface* surface =                  << 
320     G4LogicalBorderSurface::GetSurface(thePreP << 
321   if(surface == nullptr)                       << 
322   {                                            << 
323     if(thePostPV->GetMotherLogical() == thePre << 
324     {                                          << 
325       surface = G4LogicalSkinSurface::GetSurfa << 
326       if(surface == nullptr)                   << 
327       {                                        << 
328         surface =                              << 
329           G4LogicalSkinSurface::GetSurface(the << 
330       }                                        << 
331     }                                          << 
332     else                                       << 
333     {                                          << 
334       surface = G4LogicalSkinSurface::GetSurfa << 
335       if(surface == nullptr)                   << 
336       {                                        << 
337         surface =                              << 
338           G4LogicalSkinSurface::GetSurface(the << 
339       }                                        << 
340     }                                          << 
341   }                                            << 
342                                                   194 
343   if(surface != nullptr)                       << 195         G4ThreeVector theLocalPoint = theNavigator->
344   {                                            << 196                                       GetGlobalToLocalTransform().
345     fOpticalSurface =                          << 197                                       TransformPoint(theGlobalPoint);
346       dynamic_cast<G4OpticalSurface*>(surface- << 
347   }                                            << 
348   if(fOpticalSurface != nullptr)               << 
349   {                                            << 
350     type    = fOpticalSurface->GetType();      << 
351     fModel  = fOpticalSurface->GetModel();     << 
352     fFinish = fOpticalSurface->GetFinish();    << 
353                                                << 
354     G4MaterialPropertiesTable* sMPT =          << 
355       fOpticalSurface->GetMaterialPropertiesTa << 
356     if(sMPT != nullptr)                        << 
357     {                                          << 
358       if(fFinish == polishedbackpainted || fFi << 
359       {                                        << 
360         rIndexMPV = sMPT->GetProperty(kRINDEX) << 
361         if(rIndexMPV != nullptr)               << 
362         {                                      << 
363           fRindex2 = rIndexMPV->Value(fPhotonM << 
364         }                                      << 
365         else                                   << 
366         {                                      << 
367           fStatus = NoRINDEX;                  << 
368           if(verboseLevel > 1)                 << 
369             BoundaryProcessVerbose();          << 
370           aParticleChange.ProposeLocalEnergyDe << 
371           aParticleChange.ProposeTrackStatus(f << 
372           return G4VDiscreteProcess::PostStepD << 
373         }                                      << 
374       }                                        << 
375                                                   198 
376       fRealRIndexMPV = sMPT->GetProperty(kREAL << 199         G4ThreeVector theLocalNormal;   // Normal points back into volume
377       fImagRIndexMPV = sMPT->GetProperty(kIMAG << 200 
378       f_iTE = f_iTM = 1;                       << 201         G4bool valid;
379                                                << 202         theLocalNormal = theNavigator->GetLocalExitNormal(&valid);
380       G4MaterialPropertyVector* pp;            << 203 
381       if((pp = sMPT->GetProperty(kREFLECTIVITY << 204         if (valid) {
382       {                                        << 205           theLocalNormal = -theLocalNormal;
383         fReflectivity = pp->Value(fPhotonMomen << 
384       }                                        << 
385       else if(fRealRIndexMPV && fImagRIndexMPV << 
386       {                                        << 
387         CalculateReflectivity();               << 
388       }                                        << 
389                                                << 
390       if((pp = sMPT->GetProperty(kEFFICIENCY)) << 
391       {                                        << 
392         fEfficiency = pp->Value(fPhotonMomentu << 
393       }                                        << 
394       if((pp = sMPT->GetProperty(kTRANSMITTANC << 
395       {                                        << 
396         fTransmittance = pp->Value(fPhotonMome << 
397       }                                        << 
398       if(sMPT->ConstPropertyExists(kSURFACEROU << 
399       {                                        << 
400         fSurfaceRoughness = sMPT->GetConstProp << 
401       }                                        << 
402                                                << 
403       if(fModel == unified)                    << 
404       {                                        << 
405         fProb_sl = (pp = sMPT->GetProperty(kSP << 
406                      ? pp->Value(fPhotonMoment << 
407                      : 0.;                     << 
408         fProb_ss = (pp = sMPT->GetProperty(kSP << 
409                      ? pp->Value(fPhotonMoment << 
410                      : 0.;                     << 
411         fProb_bs = (pp = sMPT->GetProperty(kBA << 
412                      ? pp->Value(fPhotonMoment << 
413                      : 0.;                     << 
414       }                                        << 
415     }  // end of if(sMPT)                      << 
416     else if(fFinish == polishedbackpainted ||  << 
417     {                                          << 
418       aParticleChange.ProposeLocalEnergyDeposi << 
419       aParticleChange.ProposeTrackStatus(fStop << 
420       return G4VDiscreteProcess::PostStepDoIt( << 
421     }                                          << 
422   }  // end of if(fOpticalSurface)             << 
423                                                << 
424   //  DIELECTRIC-DIELECTRIC                    << 
425   if(type == dielectric_dielectric)            << 
426   {                                            << 
427     if(fFinish == polished || fFinish == groun << 
428     {                                          << 
429       if(fMaterial1 == fMaterial2)             << 
430       {                                        << 
431         fStatus = SameMaterial;                << 
432         if(verboseLevel > 1)                   << 
433           BoundaryProcessVerbose();            << 
434         return G4VDiscreteProcess::PostStepDoI << 
435       }                                        << 
436       MPT       = fMaterial2->GetMaterialPrope << 
437       rIndexMPV = nullptr;                     << 
438       if(MPT != nullptr)                       << 
439       {                                        << 
440         rIndexMPV = MPT->GetProperty(kRINDEX); << 
441       }                                        << 
442       if(rIndexMPV != nullptr)                 << 
443       {                                        << 
444         fRindex2 = rIndexMPV->Value(fPhotonMom << 
445       }                                        << 
446       else                                     << 
447       {                                        << 
448         fStatus = NoRINDEX;                    << 
449         if(verboseLevel > 1)                   << 
450           BoundaryProcessVerbose();            << 
451         aParticleChange.ProposeLocalEnergyDepo << 
452         aParticleChange.ProposeTrackStatus(fSt << 
453         return G4VDiscreteProcess::PostStepDoI << 
454       }                                        << 
455     }                                          << 
456     if(fFinish == polishedbackpainted || fFini << 
457     {                                          << 
458       DielectricDielectric();                  << 
459     }                                          << 
460     else                                       << 
461     {                                          << 
462       G4double rand = G4UniformRand();         << 
463       if(rand > fReflectivity + fTransmittance << 
464       {                                        << 
465         DoAbsorption();                        << 
466       }                                        << 
467       else if(rand > fReflectivity)            << 
468       {                                        << 
469         fStatus          = Transmission;       << 
470         fNewMomentum     = fOldMomentum;       << 
471         fNewPolarization = fOldPolarization;   << 
472       }                                        << 
473       else                                     << 
474       {                                        << 
475         if(fFinish == polishedfrontpainted)    << 
476         {                                      << 
477           DoReflection();                      << 
478         }                                      << 
479         else if(fFinish == groundfrontpainted) << 
480         {                                      << 
481           fStatus = LambertianReflection;      << 
482           DoReflection();                      << 
483         }                                         206         }
484         else                                   << 207         else {
485         {                                      << 208           G4ExceptionDescription ed;
486           DielectricDielectric();              << 209           ed << " G4OpBoundaryProcess/PostStepDoIt(): "
                                                   >> 210                  << " The Navigator reports that it returned an invalid normal"
                                                   >> 211                  << G4endl;
                                                   >> 212           G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
                                                   >> 213                       EventMustBeAborted,ed,
                                                   >> 214                       "Invalid Surface Normal - Geometry must return valid surface normal");
487         }                                         215         }
488       }                                        << 
489     }                                          << 
490   }                                            << 
491   else if(type == dielectric_metal)            << 
492   {                                            << 
493     DielectricMetal();                         << 
494   }                                            << 
495   else if(type == dielectric_LUT)              << 
496   {                                            << 
497     DielectricLUT();                           << 
498   }                                            << 
499   else if(type == dielectric_LUTDAVIS)         << 
500   {                                            << 
501     DielectricLUTDAVIS();                      << 
502   }                                            << 
503   else if(type == dielectric_dichroic)         << 
504   {                                            << 
505     DielectricDichroic();                      << 
506   }                                            << 
507   else if(type == coated)                      << 
508   {                                            << 
509     CoatedDielectricDielectric();              << 
510   }                                            << 
511   else                                         << 
512   {                                            << 
513     if(fNumBdryTypeWarnings <= 10)             << 
514     {                                          << 
515       ++fNumBdryTypeWarnings;                  << 
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 << 
528   }                                            << 
529                                                   216 
530   fNewMomentum     = fNewMomentum.unit();      << 217         theGlobalNormal = theNavigator->GetLocalToGlobalTransform().
531   fNewPolarization = fNewPolarization.unit();  << 218                                         TransformAxis(theLocalNormal);
532                                                   219 
533   if(verboseLevel > 1)                         << 220         if (OldMomentum * theGlobalNormal > 0.0) {
534   {                                            << 221 #ifdef G4DEBUG_OPTICAL
535     G4cout << " New Momentum Direction: " << f << 222            G4cerr << " G4OpBoundaryProcess/PostStepDoIt(): "
536            << " New Polarization:       " << f << 223                   << " theGlobalNormal points the wrong direction "
537     BoundaryProcessVerbose();                  << 224                   << G4endl;
538   }                                            << 225 #endif
539                                                << 226            theGlobalNormal = -theGlobalNormal;
540   aParticleChange.ProposeMomentumDirection(fNe << 227         }
541   aParticleChange.ProposePolarization(fNewPola << 
542                                                   228 
543   if(fStatus == FresnelRefraction || fStatus = << 229   G4MaterialPropertiesTable* aMaterialPropertiesTable;
544   {                                            << 230         G4MaterialPropertyVector* Rindex;
545     // not all surface types check that fMater << 
546     G4MaterialPropertiesTable* aMPT = fMateria << 
547     G4MaterialPropertyVector* groupvel = nullp << 
548     if(aMPT != nullptr)                        << 
549     {                                          << 
550       groupvel = aMPT->GetProperty(kGROUPVEL); << 
551     }                                          << 
552     if(groupvel != nullptr)                    << 
553     {                                          << 
554       aParticleChange.ProposeVelocity(         << 
555         groupvel->Value(fPhotonMomentum, idx_g << 
556     }                                          << 
557   }                                            << 
558                                                   231 
559   if(fStatus == Detection && fInvokeSD)        << 232   aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
560     InvokeSD(pStep);                           << 233         if (aMaterialPropertiesTable) {
561   return G4VDiscreteProcess::PostStepDoIt(aTra << 234     Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
562 }                                              << 235   }
                                                   >> 236   else {
                                                   >> 237           theStatus = NoRINDEX;
                                                   >> 238                 if ( verboseLevel > 0) BoundaryProcessVerbose();
                                                   >> 239                 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
                                                   >> 240     aParticleChange.ProposeTrackStatus(fStopAndKill);
                                                   >> 241     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 242   }
                                                   >> 243 
                                                   >> 244         if (Rindex) {
                                                   >> 245     Rindex1 = Rindex->Value(thePhotonMomentum);
                                                   >> 246   }
                                                   >> 247   else {
                                                   >> 248           theStatus = NoRINDEX;
                                                   >> 249                 if ( verboseLevel > 0) BoundaryProcessVerbose();
                                                   >> 250                 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
                                                   >> 251     aParticleChange.ProposeTrackStatus(fStopAndKill);
                                                   >> 252     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 253   }
                                                   >> 254 
                                                   >> 255         theReflectivity =  1.;
                                                   >> 256         theEfficiency   =  0.;
                                                   >> 257         theTransmittance = 0.;
                                                   >> 258 
                                                   >> 259         theModel = glisur;
                                                   >> 260         theFinish = polished;
                                                   >> 261 
                                                   >> 262         G4SurfaceType type = dielectric_dielectric;
                                                   >> 263 
                                                   >> 264         Rindex = NULL;
                                                   >> 265         OpticalSurface = NULL;
                                                   >> 266 
                                                   >> 267         G4LogicalSurface* Surface = NULL;
                                                   >> 268 
                                                   >> 269         Surface = G4LogicalBorderSurface::GetSurface
                                                   >> 270             (pPreStepPoint ->GetPhysicalVolume(),
                                                   >> 271              pPostStepPoint->GetPhysicalVolume());
                                                   >> 272 
                                                   >> 273         if (Surface == NULL){
                                                   >> 274     G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
                                                   >> 275           ->GetMotherLogical() ==
                                                   >> 276           pPreStepPoint->GetPhysicalVolume()
                                                   >> 277           ->GetLogicalVolume());
                                                   >> 278     if(enteredDaughter){
                                                   >> 279       Surface = G4LogicalSkinSurface::GetSurface
                                                   >> 280         (pPostStepPoint->GetPhysicalVolume()->
                                                   >> 281          GetLogicalVolume());
                                                   >> 282       if(Surface == NULL)
                                                   >> 283         Surface = G4LogicalSkinSurface::GetSurface
                                                   >> 284         (pPreStepPoint->GetPhysicalVolume()->
                                                   >> 285          GetLogicalVolume());
                                                   >> 286     }
                                                   >> 287     else {
                                                   >> 288       Surface = G4LogicalSkinSurface::GetSurface
                                                   >> 289         (pPreStepPoint->GetPhysicalVolume()->
                                                   >> 290          GetLogicalVolume());
                                                   >> 291       if(Surface == NULL)
                                                   >> 292         Surface = G4LogicalSkinSurface::GetSurface
                                                   >> 293         (pPostStepPoint->GetPhysicalVolume()->
                                                   >> 294          GetLogicalVolume());
                                                   >> 295     }
                                                   >> 296   }
                                                   >> 297 
                                                   >> 298   if (Surface) OpticalSurface = 
                                                   >> 299            dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
                                                   >> 300 
                                                   >> 301   if (OpticalSurface) {
                                                   >> 302 
                                                   >> 303            type      = OpticalSurface->GetType();
                                                   >> 304      theModel  = OpticalSurface->GetModel();
                                                   >> 305      theFinish = OpticalSurface->GetFinish();
                                                   >> 306 
                                                   >> 307      aMaterialPropertiesTable = OpticalSurface->
                                                   >> 308           GetMaterialPropertiesTable();
                                                   >> 309 
                                                   >> 310            if (aMaterialPropertiesTable) {
                                                   >> 311 
                                                   >> 312               if (theFinish == polishedbackpainted ||
                                                   >> 313                   theFinish == groundbackpainted ) {
                                                   >> 314                   Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
                                                   >> 315             if (Rindex) {
                                                   >> 316                      Rindex2 = Rindex->Value(thePhotonMomentum);
                                                   >> 317                   }
                                                   >> 318                   else {
                                                   >> 319          theStatus = NoRINDEX;
                                                   >> 320                      if ( verboseLevel > 0) BoundaryProcessVerbose();
                                                   >> 321                      aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
                                                   >> 322                      aParticleChange.ProposeTrackStatus(fStopAndKill);
                                                   >> 323                      return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 324                   }
                                                   >> 325               }
                                                   >> 326 
                                                   >> 327               PropertyPointer =
                                                   >> 328                       aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
                                                   >> 329               PropertyPointer1 =
                                                   >> 330                       aMaterialPropertiesTable->GetProperty("REALRINDEX");
                                                   >> 331               PropertyPointer2 =
                                                   >> 332                       aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
                                                   >> 333 
                                                   >> 334               iTE = 1;
                                                   >> 335               iTM = 1;
                                                   >> 336 
                                                   >> 337               if (PropertyPointer) {
                                                   >> 338 
                                                   >> 339                  theReflectivity =
                                                   >> 340                           PropertyPointer->Value(thePhotonMomentum);
                                                   >> 341 
                                                   >> 342               } else if (PropertyPointer1 && PropertyPointer2) {
                                                   >> 343 
                                                   >> 344                  CalculateReflectivity();
                                                   >> 345 
                                                   >> 346               }
                                                   >> 347 
                                                   >> 348               PropertyPointer =
                                                   >> 349               aMaterialPropertiesTable->GetProperty("EFFICIENCY");
                                                   >> 350               if (PropertyPointer) {
                                                   >> 351                       theEfficiency =
                                                   >> 352                       PropertyPointer->Value(thePhotonMomentum);
                                                   >> 353               }
                                                   >> 354 
                                                   >> 355               PropertyPointer =
                                                   >> 356               aMaterialPropertiesTable->GetProperty("TRANSMITTANCE");
                                                   >> 357               if (PropertyPointer) {
                                                   >> 358                       theTransmittance =
                                                   >> 359                       PropertyPointer->Value(thePhotonMomentum);
                                                   >> 360               }
                                                   >> 361 
                                                   >> 362         if ( theModel == unified ) {
                                                   >> 363           PropertyPointer =
                                                   >> 364     aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
                                                   >> 365           if (PropertyPointer) {
                                                   >> 366                          prob_sl =
                                                   >> 367        PropertyPointer->Value(thePhotonMomentum);
                                                   >> 368                 } else {
                                                   >> 369                          prob_sl = 0.0;
                                                   >> 370                 }
                                                   >> 371 
                                                   >> 372           PropertyPointer =
                                                   >> 373     aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
                                                   >> 374           if (PropertyPointer) {
                                                   >> 375                          prob_ss =
                                                   >> 376        PropertyPointer->Value(thePhotonMomentum);
                                                   >> 377                 } else {
                                                   >> 378                          prob_ss = 0.0;
                                                   >> 379                 }
                                                   >> 380 
                                                   >> 381           PropertyPointer =
                                                   >> 382     aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
                                                   >> 383           if (PropertyPointer) {
                                                   >> 384                          prob_bs =
                                                   >> 385        PropertyPointer->Value(thePhotonMomentum);
                                                   >> 386                 } else {
                                                   >> 387                          prob_bs = 0.0;
                                                   >> 388                 }
                                                   >> 389         }
                                                   >> 390      }
                                                   >> 391            else if (theFinish == polishedbackpainted ||
                                                   >> 392                     theFinish == groundbackpainted ) {
                                                   >> 393                       aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
                                                   >> 394                       aParticleChange.ProposeTrackStatus(fStopAndKill);
                                                   >> 395                       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 396            }
                                                   >> 397         }
                                                   >> 398 
                                                   >> 399         if (type == dielectric_dielectric ) {
                                                   >> 400            if (theFinish == polished || theFinish == ground ) {
                                                   >> 401 
                                                   >> 402         if (Material1 == Material2){
                                                   >> 403      theStatus = SameMaterial;
                                                   >> 404                  if ( verboseLevel > 0) BoundaryProcessVerbose();
                                                   >> 405      return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 406         }
                                                   >> 407               aMaterialPropertiesTable =
                                                   >> 408                      Material2->GetMaterialPropertiesTable();
                                                   >> 409               if (aMaterialPropertiesTable)
                                                   >> 410                  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
                                                   >> 411               if (Rindex) {
                                                   >> 412                  Rindex2 = Rindex->Value(thePhotonMomentum);
                                                   >> 413               }
                                                   >> 414               else {
                                                   >> 415      theStatus = NoRINDEX;
                                                   >> 416                  if ( verboseLevel > 0) BoundaryProcessVerbose();
                                                   >> 417                  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
                                                   >> 418                  aParticleChange.ProposeTrackStatus(fStopAndKill);
                                                   >> 419                  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 420         }
                                                   >> 421            }
                                                   >> 422         }
                                                   >> 423 
                                                   >> 424   if (type == dielectric_metal) {
                                                   >> 425 
                                                   >> 426     DielectricMetal();
                                                   >> 427 
                                                   >> 428           // Uncomment the following lines if you wish to have 
                                                   >> 429           //         Transmission instead of Absorption
                                                   >> 430           // if (theStatus == Absorption) {
                                                   >> 431           //    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 432           // }
                                                   >> 433 
                                                   >> 434   }
                                                   >> 435         else if (type == dielectric_LUT) {
                                                   >> 436 
                                                   >> 437           DielectricLUT();
                                                   >> 438 
                                                   >> 439         }
                                                   >> 440   else if (type == dielectric_dielectric) {
                                                   >> 441 
                                                   >> 442           if ( theFinish == polishedbackpainted ||
                                                   >> 443                theFinish == groundbackpainted ) {
                                                   >> 444              DielectricDielectric();
                                                   >> 445           }
                                                   >> 446     else {
                                                   >> 447              if ( !G4BooleanRand(theReflectivity) ) {
                                                   >> 448                 DoAbsorption();
                                                   >> 449              }
                                                   >> 450              else {
                                                   >> 451                 if ( theFinish == polishedfrontpainted ) {
                                                   >> 452                    DoReflection();
                                                   >> 453                 }
                                                   >> 454                 else if ( theFinish == groundfrontpainted ) {
                                                   >> 455                    theStatus = LambertianReflection;
                                                   >> 456                    DoReflection();
                                                   >> 457                 }
                                                   >> 458                 else {
                                                   >> 459                    DielectricDielectric();
                                                   >> 460                 }
                                                   >> 461              }
                                                   >> 462           }
                                                   >> 463         }
                                                   >> 464   else {
563                                                   465 
564 //....oooOO0OOooo........oooOO0OOooo........oo << 466     G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
565 void G4OpBoundaryProcess::BoundaryProcessVerbo << 467     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
566 {                                              << 
567   G4cout << " *** ";                           << 
568   if(fStatus == Undefined)                     << 
569     G4cout << "Undefined";                     << 
570   else if(fStatus == Transmission)             << 
571     G4cout << "Transmission";                  << 
572   else if(fStatus == FresnelRefraction)        << 
573     G4cout << "FresnelRefraction";             << 
574   else if(fStatus == FresnelReflection)        << 
575     G4cout << "FresnelReflection";             << 
576   else if(fStatus == TotalInternalReflection)  << 
577     G4cout << "TotalInternalReflection";       << 
578   else if(fStatus == LambertianReflection)     << 
579     G4cout << "LambertianReflection";          << 
580   else if(fStatus == LobeReflection)           << 
581     G4cout << "LobeReflection";                << 
582   else if(fStatus == SpikeReflection)          << 
583     G4cout << "SpikeReflection";               << 
584   else if(fStatus == BackScattering)           << 
585     G4cout << "BackScattering";                << 
586   else if(fStatus == PolishedLumirrorAirReflec << 
587     G4cout << "PolishedLumirrorAirReflection"; << 
588   else if(fStatus == PolishedLumirrorGlueRefle << 
589     G4cout << "PolishedLumirrorGlueReflection" << 
590   else if(fStatus == PolishedAirReflection)    << 
591     G4cout << "PolishedAirReflection";         << 
592   else if(fStatus == PolishedTeflonAirReflecti << 
593     G4cout << "PolishedTeflonAirReflection";   << 
594   else if(fStatus == PolishedTiOAirReflection) << 
595     G4cout << "PolishedTiOAirReflection";      << 
596   else if(fStatus == PolishedTyvekAirReflectio << 
597     G4cout << "PolishedTyvekAirReflection";    << 
598   else if(fStatus == PolishedVM2000AirReflecti << 
599     G4cout << "PolishedVM2000AirReflection";   << 
600   else if(fStatus == PolishedVM2000GlueReflect << 
601     G4cout << "PolishedVM2000GlueReflection";  << 
602   else if(fStatus == EtchedLumirrorAirReflecti << 
603     G4cout << "EtchedLumirrorAirReflection";   << 
604   else if(fStatus == EtchedLumirrorGlueReflect << 
605     G4cout << "EtchedLumirrorGlueReflection";  << 
606   else if(fStatus == EtchedAirReflection)      << 
607     G4cout << "EtchedAirReflection";           << 
608   else if(fStatus == EtchedTeflonAirReflection << 
609     G4cout << "EtchedTeflonAirReflection";     << 
610   else if(fStatus == EtchedTiOAirReflection)   << 
611     G4cout << "EtchedTiOAirReflection";        << 
612   else if(fStatus == EtchedTyvekAirReflection) << 
613     G4cout << "EtchedTyvekAirReflection";      << 
614   else if(fStatus == EtchedVM2000AirReflection << 
615     G4cout << "EtchedVM2000AirReflection";     << 
616   else if(fStatus == EtchedVM2000GlueReflectio << 
617     G4cout << "EtchedVM2000GlueReflection";    << 
618   else if(fStatus == GroundLumirrorAirReflecti << 
619     G4cout << "GroundLumirrorAirReflection";   << 
620   else if(fStatus == GroundLumirrorGlueReflect << 
621     G4cout << "GroundLumirrorGlueReflection";  << 
622   else if(fStatus == GroundAirReflection)      << 
623     G4cout << "GroundAirReflection";           << 
624   else if(fStatus == GroundTeflonAirReflection << 
625     G4cout << "GroundTeflonAirReflection";     << 
626   else if(fStatus == GroundTiOAirReflection)   << 
627     G4cout << "GroundTiOAirReflection";        << 
628   else if(fStatus == GroundTyvekAirReflection) << 
629     G4cout << "GroundTyvekAirReflection";      << 
630   else if(fStatus == GroundVM2000AirReflection << 
631     G4cout << "GroundVM2000AirReflection";     << 
632   else if(fStatus == GroundVM2000GlueReflectio << 
633     G4cout << "GroundVM2000GlueReflection";    << 
634   else if(fStatus == Absorption)               << 
635     G4cout << "Absorption";                    << 
636   else if(fStatus == Detection)                << 
637     G4cout << "Detection";                     << 
638   else if(fStatus == NotAtBoundary)            << 
639     G4cout << "NotAtBoundary";                 << 
640   else if(fStatus == SameMaterial)             << 
641     G4cout << "SameMaterial";                  << 
642   else if(fStatus == StepTooSmall)             << 
643     G4cout << "StepTooSmall";                  << 
644   else if(fStatus == NoRINDEX)                 << 
645     G4cout << "NoRINDEX";                      << 
646   else if(fStatus == Dichroic)                 << 
647     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                                                   468 
655   G4cout << " ***" << G4endl;                  << 469   }
656 }                                              << 
657                                                   470 
658 //....oooOO0OOooo........oooOO0OOooo........oo << 471         NewMomentum = NewMomentum.unit();
659 G4ThreeVector G4OpBoundaryProcess::GetFacetNor << 472         NewPolarization = NewPolarization.unit();
660   const G4ThreeVector& momentum, const G4Three << 
661 {                                              << 
662   G4ThreeVector facetNormal;                   << 
663   if(fModel == unified || fModel == LUT || fMo << 
664   {                                            << 
665     /* This function codes alpha to a random v << 
666     distribution p(alpha) = g(alpha; 0, sigma_ << 
667     for alpha > 0 and alpha < 90, where g(alph << 
668     gaussian distribution with mean 0 and stan << 
669                                                << 
670     G4double sigma_alpha = 0.0;                << 
671     if(fOpticalSurface)                        << 
672       sigma_alpha = fOpticalSurface->GetSigmaA << 
673     if(sigma_alpha == 0.0)                     << 
674     {                                          << 
675       return normal;                           << 
676     }                                          << 
677                                                << 
678     G4double f_max = std::min(1.0, 4. * sigma_ << 
679     G4double alpha, phi, sinAlpha;             << 
680                                                << 
681     do                                         << 
682     {  // Loop checking, 13-Aug-2015, Peter Gu << 
683       do                                       << 
684       {  // Loop checking, 13-Aug-2015, Peter  << 
685         alpha    = G4RandGauss::shoot(0.0, sig << 
686         sinAlpha = std::sin(alpha);            << 
687       } while(G4UniformRand() * f_max > sinAlp << 
688                                                << 
689       phi = G4UniformRand() * twopi;           << 
690       facetNormal.set(sinAlpha * std::cos(phi) << 
691                       std::cos(alpha));        << 
692       facetNormal.rotateUz(normal);            << 
693     } while(momentum * facetNormal >= 0.0);    << 
694   }                                            << 
695   else                                         << 
696   {                                            << 
697     G4double polish = 1.0;                     << 
698     if(fOpticalSurface)                        << 
699       polish = fOpticalSurface->GetPolish();   << 
700     if(polish < 1.0)                           << 
701     {                                          << 
702       do                                       << 
703       {  // Loop checking, 13-Aug-2015, Peter  << 
704         G4ThreeVector smear;                   << 
705         do                                     << 
706         {  // Loop checking, 13-Aug-2015, Pete << 
707           smear.setX(2. * G4UniformRand() - 1. << 
708           smear.setY(2. * G4UniformRand() - 1. << 
709           smear.setZ(2. * G4UniformRand() - 1. << 
710         } while(smear.mag2() > 1.0);           << 
711         facetNormal = normal + (1. - polish) * << 
712       } while(momentum * facetNormal >= 0.0);  << 
713       facetNormal = facetNormal.unit();        << 
714     }                                          << 
715     else                                       << 
716     {                                          << 
717       facetNormal = normal;                    << 
718     }                                          << 
719   }                                            << 
720   return facetNormal;                          << 
721 }                                              << 
722                                                   473 
723 //....oooOO0OOooo........oooOO0OOooo........oo << 474         if ( verboseLevel > 0) {
724 void G4OpBoundaryProcess::DielectricMetal()    << 475      G4cout << " New Momentum Direction: " << NewMomentum     << G4endl;
725 {                                              << 476      G4cout << " New Polarization:       " << NewPolarization << G4endl;
726   G4int n = 0;                                 << 477            BoundaryProcessVerbose();
727   G4double rand;                               << 
728   G4ThreeVector A_trans;                       << 
729                                                << 
730   do                                           << 
731   {                                            << 
732     ++n;                                       << 
733     rand = G4UniformRand();                    << 
734     if(rand > fReflectivity && n == 1)         << 
735     {                                          << 
736       if(rand > fReflectivity + fTransmittance << 
737       {                                        << 
738         DoAbsorption();                        << 
739       }                                        << 
740       else                                     << 
741       {                                        << 
742         fStatus          = Transmission;       << 
743         fNewMomentum     = fOldMomentum;       << 
744         fNewPolarization = fOldPolarization;   << 
745       }                                        << 
746       break;                                   << 
747     }                                          << 
748     else                                       << 
749     {                                          << 
750       if(fRealRIndexMPV && fImagRIndexMPV)     << 
751       {                                        << 
752         if(n > 1)                              << 
753         {                                      << 
754           CalculateReflectivity();             << 
755           if(!G4BooleanRand(fReflectivity))    << 
756           {                                    << 
757             DoAbsorption();                    << 
758             break;                             << 
759           }                                    << 
760         }                                         478         }
761       }                                        << 
762       if(fModel == glisur || fFinish == polish << 
763       {                                        << 
764         DoReflection();                        << 
765       }                                        << 
766       else                                     << 
767       {                                        << 
768         if(n == 1)                             << 
769           ChooseReflection();                  << 
770         if(fStatus == LambertianReflection)    << 
771         {                                      << 
772           DoReflection();                      << 
773         }                                      << 
774         else if(fStatus == BackScattering)     << 
775         {                                      << 
776           fNewMomentum     = -fOldMomentum;    << 
777           fNewPolarization = -fOldPolarization << 
778         }                                      << 
779         else                                   << 
780         {                                      << 
781           if(fStatus == LobeReflection)        << 
782           {                                    << 
783             if(!fRealRIndexMPV || !fImagRIndex << 
784             {                                  << 
785               fFacetNormal = GetFacetNormal(fO << 
786             }                                  << 
787             // else                            << 
788             //  case of complex rindex needs t << 
789           }                                    << 
790           fNewMomentum =                       << 
791             fOldMomentum - 2. * fOldMomentum * << 
792                                                   479 
793           if(f_iTE > 0 && f_iTM > 0)           << 480   aParticleChange.ProposeMomentumDirection(NewMomentum);
794           {                                    << 481   aParticleChange.ProposePolarization(NewPolarization);
795             fNewPolarization =                 << 482 
796               -fOldPolarization +              << 483         if ( theStatus == FresnelRefraction ) {
797               (2. * fOldPolarization * fFacetN << 484            G4MaterialPropertyVector* groupvel =
798           }                                    << 485            Material2->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
799           else if(f_iTE > 0)                   << 486            G4double finalVelocity = groupvel->Value(thePhotonMomentum);
800           {                                    << 487            aParticleChange.ProposeVelocity(finalVelocity);
801             A_trans = (fSint1 > 0.0) ? fOldMom << 
802                                      : fOldPol << 
803             fNewPolarization = -A_trans;       << 
804           }                                    << 
805           else if(f_iTM > 0)                   << 
806           {                                    << 
807             fNewPolarization =                 << 
808               -fNewMomentum.cross(A_trans).uni << 
809           }                                    << 
810         }                                         488         }
811       }                                        << 
812       fOldMomentum     = fNewMomentum;         << 
813       fOldPolarization = fNewPolarization;     << 
814     }                                          << 
815     // Loop checking, 13-Aug-2015, Peter Gumpl << 
816   } while(fNewMomentum * fGlobalNormal < 0.0); << 
817 }                                              << 
818                                                   489 
819 //....oooOO0OOooo........oooOO0OOooo........oo << 490         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
820 void G4OpBoundaryProcess::DielectricLUT()      << 
821 {                                              << 
822   G4int thetaIndex, phiIndex;                  << 
823   G4double angularDistVal, thetaRad, phiRad;   << 
824   G4ThreeVector perpVectorTheta, perpVectorPhi << 
825                                                << 
826   fStatus = G4OpBoundaryProcessStatus(         << 
827     G4int(fFinish) + (G4int(NoRINDEX) - G4int( << 
828                                                << 
829   G4int thetaIndexMax = fOpticalSurface->GetTh << 
830   G4int phiIndexMax   = fOpticalSurface->GetPh << 
831                                                << 
832   G4double rand;                               << 
833                                                << 
834   do                                           << 
835   {                                            << 
836     rand = G4UniformRand();                    << 
837     if(rand > fReflectivity)                   << 
838     {                                          << 
839       if(rand > fReflectivity + fTransmittance << 
840       {                                        << 
841         DoAbsorption();                        << 
842       }                                        << 
843       else                                     << 
844       {                                        << 
845         fStatus          = Transmission;       << 
846         fNewMomentum     = fOldMomentum;       << 
847         fNewPolarization = fOldPolarization;   << 
848       }                                        << 
849       break;                                   << 
850     }                                          << 
851     else                                       << 
852     {                                          << 
853       // Calculate Angle between Normal and Ph << 
854       G4double anglePhotonToNormal = fOldMomen << 
855       // Round to closest integer: LBNL model  << 
856       G4int angleIncident = (G4int)std::lrint( << 
857                                                << 
858       // Take random angles THETA and PHI,     << 
859       // and see if below Probability - if not << 
860       do                                       << 
861       {                                        << 
862         thetaIndex = (G4int)G4RandFlat::shootI << 
863         phiIndex   = (G4int)G4RandFlat::shootI << 
864         // Find probability with the new indec << 
865         angularDistVal = fOpticalSurface->GetA << 
866           angleIncident, thetaIndex, phiIndex) << 
867         // Loop checking, 13-Aug-2015, Peter G << 
868       } while(!G4BooleanRand(angularDistVal)); << 
869                                                << 
870       thetaRad = G4double(-90 + 4 * thetaIndex << 
871       phiRad   = G4double(-90 + 5 * phiIndex)  << 
872       // Rotate Photon Momentum in Theta, then << 
873       fNewMomentum = -fOldMomentum;            << 
874                                                << 
875       perpVectorTheta = fNewMomentum.cross(fGl << 
876       if(perpVectorTheta.mag() < fCarTolerance << 
877       {                                        << 
878         perpVectorTheta = fNewMomentum.orthogo << 
879       }                                        << 
880       fNewMomentum =                           << 
881         fNewMomentum.rotate(anglePhotonToNorma << 
882       perpVectorPhi = perpVectorTheta.cross(fN << 
883       fNewMomentum  = fNewMomentum.rotate(-phi << 
884                                                << 
885       // Rotate Polarization too:              << 
886       fFacetNormal     = (fNewMomentum - fOldM << 
887       fNewPolarization = -fOldPolarization +   << 
888                          (2. * fOldPolarizatio << 
889     }                                          << 
890     // Loop checking, 13-Aug-2015, Peter Gumpl << 
891   } while(fNewMomentum * fGlobalNormal <= 0.0) << 
892 }                                                 491 }
893                                                   492 
894 //....oooOO0OOooo........oooOO0OOooo........oo << 493 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
895 void G4OpBoundaryProcess::DielectricLUTDAVIS() << 
896 {                                                 494 {
897   G4int angindex, random, angleIncident;       << 495         if ( theStatus == Undefined )
898   G4double reflectivityValue, elevation, azimu << 496                 G4cout << " *** Undefined *** " << G4endl;
899   G4double anglePhotonToNormal;                << 497         if ( theStatus == FresnelRefraction )
900                                                << 498                 G4cout << " *** FresnelRefraction *** " << G4endl;
901   G4int lutbin  = fOpticalSurface->GetLUTbins( << 499         if ( theStatus == FresnelReflection )
902   G4double rand = G4UniformRand();             << 500                 G4cout << " *** FresnelReflection *** " << G4endl;
903                                                << 501         if ( theStatus == TotalInternalReflection )
904   G4double sinEl;                              << 502                 G4cout << " *** TotalInternalReflection *** " << G4endl;
905   G4ThreeVector u, vNorm, w;                   << 503         if ( theStatus == LambertianReflection )
906                                                << 504                 G4cout << " *** LambertianReflection *** " << G4endl;
907   do                                           << 505         if ( theStatus == LobeReflection )
908   {                                            << 506                 G4cout << " *** LobeReflection *** " << G4endl;
909     anglePhotonToNormal = fOldMomentum.angle(- << 507         if ( theStatus == SpikeReflection )
910                                                << 508                 G4cout << " *** SpikeReflection *** " << G4endl;
911     // Davis model has 90 reflection bins: rou << 509         if ( theStatus == BackScattering )
912     // don't allow angleIncident to be 90 for  << 510                 G4cout << " *** BackScattering *** " << G4endl;
913     angleIncident = std::min(                  << 511         if ( theStatus == PolishedLumirrorAirReflection )
914       static_cast<G4int>(std::floor(anglePhoto << 512                 G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
915     reflectivityValue = fOpticalSurface->GetRe << 513         if ( theStatus == PolishedLumirrorGlueReflection )
916                                                << 514                 G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
917     if(rand > reflectivityValue)               << 515         if ( theStatus == PolishedAirReflection )
918     {                                          << 516                 G4cout << " *** PolishedAirReflection *** " << G4endl;
919       if(fEfficiency > 0.)                     << 517         if ( theStatus == PolishedTeflonAirReflection )
920       {                                        << 518                 G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
921         DoAbsorption();                        << 519         if ( theStatus == PolishedTiOAirReflection )
922         break;                                 << 520                 G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
923       }                                        << 521         if ( theStatus == PolishedTyvekAirReflection )
924       else                                     << 522                 G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
925       {                                        << 523         if ( theStatus == PolishedVM2000AirReflection )
926         fStatus = Transmission;                << 524                 G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
927                                                << 525         if ( theStatus == PolishedVM2000GlueReflection )
928         if(angleIncident <= 0.01)              << 526                 G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
929         {                                      << 527         if ( theStatus == EtchedLumirrorAirReflection )
930           fNewMomentum = fOldMomentum;         << 528                 G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
931           break;                               << 529         if ( theStatus == EtchedLumirrorGlueReflection )
932         }                                      << 530                 G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
933                                                << 531         if ( theStatus == EtchedAirReflection )
934         do                                     << 532                 G4cout << " *** EtchedAirReflection *** " << G4endl;
935         {                                      << 533         if ( theStatus == EtchedTeflonAirReflection )
936           random = (G4int)G4RandFlat::shootInt << 534                 G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
937           angindex =                           << 535         if ( theStatus == EtchedTiOAirReflection )
938             (((random * 2) - 1)) + angleIncide << 536                 G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
939                                                << 537         if ( theStatus == EtchedTyvekAirReflection )
940           azimuth =                            << 538                 G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
941             fOpticalSurface->GetAngularDistrib << 539         if ( theStatus == EtchedVM2000AirReflection )
942           elevation = fOpticalSurface->GetAngu << 540                 G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
943         } while(elevation == 0. && azimuth ==  << 541         if ( theStatus == EtchedVM2000GlueReflection )
944                                                << 542                 G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
945         sinEl = std::sin(elevation);           << 543         if ( theStatus == GroundLumirrorAirReflection )
946         vNorm = (fGlobalNormal.cross(fOldMomen << 544                 G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
947         u     = vNorm.cross(fGlobalNormal) * ( << 545         if ( theStatus == GroundLumirrorGlueReflection )
948         vNorm *= (sinEl * std::sin(azimuth));  << 546                 G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
949         // fGlobalNormal shouldn't be modified << 547         if ( theStatus == GroundAirReflection )
950         w            = (fGlobalNormal *= std:: << 548                 G4cout << " *** GroundAirReflection *** " << G4endl;
951         fNewMomentum = u + vNorm + w;          << 549         if ( theStatus == GroundTeflonAirReflection )
952                                                << 550                 G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
953         // Rotate Polarization too:            << 551         if ( theStatus == GroundTiOAirReflection )
954         fFacetNormal     = (fNewMomentum - fOl << 552                 G4cout << " *** GroundTiOAirReflection *** " << G4endl;
955         fNewPolarization = -fOldPolarization + << 553         if ( theStatus == GroundTyvekAirReflection )
956                                                << 554                 G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
957       }                                        << 555         if ( theStatus == GroundVM2000AirReflection )
958     }                                          << 556                 G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
959     else                                       << 557         if ( theStatus == GroundVM2000GlueReflection )
960     {                                          << 558                 G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
961       fStatus = LobeReflection;                << 559         if ( theStatus == Absorption )
962                                                << 560                 G4cout << " *** Absorption *** " << G4endl;
963       if(angleIncident == 0)                   << 561         if ( theStatus == Detection )
964       {                                        << 562                 G4cout << " *** Detection *** " << G4endl;
965         fNewMomentum = -fOldMomentum;          << 563         if ( theStatus == NotAtBoundary )
966         break;                                 << 564                 G4cout << " *** NotAtBoundary *** " << G4endl;
967       }                                        << 565         if ( theStatus == SameMaterial )
968                                                << 566                 G4cout << " *** SameMaterial *** " << G4endl;
969       do                                       << 567         if ( theStatus == StepTooSmall )
970       {                                        << 568                 G4cout << " *** StepTooSmall *** " << G4endl;
971         random   = (G4int)G4RandFlat::shootInt << 569         if ( theStatus == NoRINDEX )
972         angindex = (((random * 2) - 1)) + (ang << 570                 G4cout << " *** NoRINDEX *** " << G4endl;
973                                                << 
974         azimuth = fOpticalSurface->GetAngularD << 
975         elevation = fOpticalSurface->GetAngula << 
976       } while(elevation == 0. && azimuth == 0. << 
977                                                << 
978       sinEl = std::sin(elevation);             << 
979       vNorm = (fGlobalNormal.cross(fOldMomentu << 
980       u     = vNorm.cross(fGlobalNormal) * (si << 
981       vNorm *= (sinEl * std::sin(azimuth));    << 
982       // fGlobalNormal shouldn't be modified h << 
983       w = (fGlobalNormal *= std::cos(elevation << 
984                                                << 
985       fNewMomentum = u + vNorm + w;            << 
986                                                << 
987       // Rotate Polarization too: (needs revis << 
988       fNewPolarization = fOldPolarization;     << 
989     }                                          << 
990   } while(fNewMomentum * fGlobalNormal <= 0.0) << 
991 }                                                 571 }
992                                                   572 
993 //....oooOO0OOooo........oooOO0OOooo........oo << 573 G4ThreeVector
994 void G4OpBoundaryProcess::DielectricDichroic() << 574 G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
                                                   >> 575                   const G4ThreeVector&  Normal ) const
995 {                                                 576 {
996   // Calculate Angle between Normal and Photon << 577         G4ThreeVector FacetNormal;
997   G4double anglePhotonToNormal = fOldMomentum. << 
998                                                   578 
999   // Round it to closest integer               << 579   if (theModel == unified || theModel == LUT) {
1000   G4double angleIncident = std::floor(180. /  << 
1001                                                  580 
1002   if(!fDichroicVector)                        << 581   /* This function code alpha to a random value taken from the
1003   {                                           << 582            distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
1004     if(fOpticalSurface)                       << 583            for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
1005       fDichroicVector = fOpticalSurface->GetD << 584            is a gaussian distribution with mean 0 and standard deviation
1006   }                                           << 585            sigma_alpha.  */
1007                                               << 586 
1008   if(fDichroicVector)                         << 587      G4double alpha;
1009   {                                           << 588 
1010     G4double wavelength = h_Planck * c_light  << 589      G4double sigma_alpha = 0.0;
1011     fTransmittance      = fDichroicVector->Va << 590      if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
1012                                             i << 591 
1013                      perCent;                 << 592      G4double f_max = std::min(1.0,4.*sigma_alpha);
1014     //   G4cout << "wavelength: " << std::flo << 593 
1015     //                            << "nm" <<  << 594      do {
1016     //   G4cout << "Incident angle: " << angl << 595         do {
1017     //   G4cout << "Transmittance: "          << 596            alpha = G4RandGauss::shoot(0.0,sigma_alpha);
1018     //          << std::floor(fTransmittance/ << 597         } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
1019   }                                           << 598 
1020   else                                        << 599         G4double phi = G4UniformRand()*twopi;
1021   {                                           << 600 
1022     G4ExceptionDescription ed;                << 601         G4double SinAlpha = std::sin(alpha);
1023     ed << " G4OpBoundaryProcess/DielectricDic << 602         G4double CosAlpha = std::cos(alpha);
1024        << " The dichroic surface has no G4Phy << 603               G4double SinPhi = std::sin(phi);
1025     G4Exception("G4OpBoundaryProcess::Dielect << 604               G4double CosPhi = std::cos(phi);
1026                 FatalException, ed,           << 605 
1027                 "A dichroic surface must have << 606               G4double unit_x = SinAlpha * CosPhi;
1028   }                                           << 607               G4double unit_y = SinAlpha * SinPhi;
1029                                               << 608               G4double unit_z = CosAlpha;
1030   if(!G4BooleanRand(fTransmittance))          << 609 
1031   {  // Not transmitted, so reflect           << 610         FacetNormal.setX(unit_x);
1032     if(fModel == glisur || fFinish == polishe << 611         FacetNormal.setY(unit_y);
1033     {                                         << 612         FacetNormal.setZ(unit_z);
1034       DoReflection();                         << 613 
1035     }                                         << 614         G4ThreeVector tmpNormal = Normal;
1036     else                                      << 615 
1037     {                                         << 616         FacetNormal.rotateUz(tmpNormal);
1038       ChooseReflection();                     << 617      } while (Momentum * FacetNormal >= 0.0);
1039       if(fStatus == LambertianReflection)     << 618   }
1040       {                                       << 619   else {
1041         DoReflection();                       << 620 
1042       }                                       << 621      G4double  polish = 1.0;
1043       else if(fStatus == BackScattering)      << 622      if (OpticalSurface) polish = OpticalSurface->GetPolish();
1044       {                                       << 623 
1045         fNewMomentum     = -fOldMomentum;     << 624            if (polish < 1.0) {
1046         fNewPolarization = -fOldPolarization; << 625               do {
1047       }                                       << 626                  G4ThreeVector smear;
1048       else                                    << 627                  do {
1049       {                                       << 628                     smear.setX(2.*G4UniformRand()-1.0);
1050         G4double PdotN, EdotN;                << 629                     smear.setY(2.*G4UniformRand()-1.0);
1051         do                                    << 630                     smear.setZ(2.*G4UniformRand()-1.0);
1052         {                                     << 631                  } while (smear.mag()>1.0);
1053           if(fStatus == LobeReflection)       << 632                  smear = (1.-polish) * smear;
1054           {                                   << 633                  FacetNormal = Normal + smear;
1055             fFacetNormal = GetFacetNormal(fOl << 634               } while (Momentum * FacetNormal >= 0.0);
1056           }                                   << 635               FacetNormal = FacetNormal.unit();
1057           PdotN        = fOldMomentum * fFace << 636      }
1058           fNewMomentum = fOldMomentum - (2. * << 637            else {
1059           // Loop checking, 13-Aug-2015, Pete << 638               FacetNormal = Normal;
1060         } while(fNewMomentum * fGlobalNormal  << 639            }
1061                                               << 640   }
1062         EdotN            = fOldPolarization * << 641   return FacetNormal;
1063         fNewPolarization = -fOldPolarization  << 
1064       }                                       << 
1065     }                                         << 
1066   }                                           << 
1067   else                                        << 
1068   {                                           << 
1069     fStatus          = Dichroic;              << 
1070     fNewMomentum     = fOldMomentum;          << 
1071     fNewPolarization = fOldPolarization;      << 
1072   }                                           << 
1073 }                                                642 }
1074                                                  643 
1075 //....oooOO0OOooo........oooOO0OOooo........o << 644 void G4OpBoundaryProcess::DielectricMetal()
1076 void G4OpBoundaryProcess::DielectricDielectri << 
1077 {                                                645 {
1078   G4bool inside = false;                      << 646         G4int n = 0;
1079   G4bool swap   = false;                      << 
1080                                                  647 
1081   if(fFinish == polished)                     << 648   do {
1082   {                                           << 
1083     fFacetNormal = fGlobalNormal;             << 
1084   }                                           << 
1085   else                                        << 
1086   {                                           << 
1087     fFacetNormal = GetFacetNormal(fOldMomentu << 
1088   }                                           << 
1089   G4double cost1 = -fOldMomentum * fFacetNorm << 
1090   G4double cost2 = 0.;                        << 
1091   G4double sint2 = 0.;                        << 
1092                                               << 
1093   G4bool surfaceRoughnessCriterionPass = true << 
1094   if(fSurfaceRoughness != 0. && fRindex1 > fR << 
1095   {                                           << 
1096     G4double wavelength                = h_Pl << 
1097     G4double surfaceRoughnessCriterion = std: << 
1098       (4. * pi * fSurfaceRoughness * fRindex1 << 
1099     surfaceRoughnessCriterionPass = G4Boolean << 
1100   }                                           << 
1101                                                  649 
1102 leap:                                         << 650            n++;
1103                                                  651 
1104   G4bool through = false;                     << 652            if( !G4BooleanRand(theReflectivity) && n == 1 ) {
1105   G4bool done    = false;                     << 
1106                                                  653 
1107   G4ThreeVector A_trans, A_paral, E1pp, E1pl; << 654              // Comment out DoAbsorption and uncomment theStatus = Absorption;
1108   G4double E1_perp, E1_parl;                  << 655              // if you wish to have Transmission instead of Absorption
1109   G4double s1, s2, E2_perp, E2_parl, E2_total << 656 
1110   G4double E2_abs, C_parl, C_perp;            << 657              DoAbsorption();
1111   G4double alpha;                             << 658              // theStatus = Absorption;
1112                                               << 659              break;
1113   do                                          << 660 
1114   {                                           << 661            }
1115     if(through)                               << 662            else {
1116     {                                         << 663 
1117       swap          = !swap;                  << 664              if (PropertyPointer1 && PropertyPointer2) {
1118       through       = false;                  << 665                 if ( n > 1 ) {
1119       fGlobalNormal = -fGlobalNormal;         << 666                    CalculateReflectivity();
1120       G4SwapPtr(fMaterial1, fMaterial2);      << 667                    if ( !G4BooleanRand(theReflectivity) ) {
1121       G4SwapObj(&fRindex1, &fRindex2);        << 668                       DoAbsorption();
1122     }                                         << 669                       break;
1123                                               << 670                    }
1124     if(fFinish == polished)                   << 671                 }
1125     {                                         << 672              }
1126       fFacetNormal = fGlobalNormal;           << 673 
1127     }                                         << 674              if ( theModel == glisur || theFinish == polished ) {
1128     else                                      << 675 
1129     {                                         << 676                 DoReflection();
1130       fFacetNormal = GetFacetNormal(fOldMomen << 677 
1131     }                                         << 678              } else {
1132                                               << 679 
1133     cost1 = -fOldMomentum * fFacetNormal;     << 680                 if ( n == 1 ) ChooseReflection();
1134     if(std::abs(cost1) < 1.0 - fCarTolerance) << 681                                                                                 
1135     {                                         << 682                 if ( theStatus == LambertianReflection ) {
1136       fSint1 = std::sqrt(1. - cost1 * cost1); << 683                    DoReflection();
1137       sint2  = fSint1 * fRindex1 / fRindex2;  << 684                 }
1138       // this isn't a sine as we might expect << 685                 else if ( theStatus == BackScattering ) {
1139     }                                         << 686                    NewMomentum = -OldMomentum;
1140     else                                      << 687                    NewPolarization = -OldPolarization;
1141     {                                         << 688                 }
1142       fSint1 = 0.0;                           << 689                 else {
1143       sint2  = 0.0;                           << 690 
1144     }                                         << 691                    if(theStatus==LobeReflection){
1145                                               << 692                      if ( PropertyPointer1 && PropertyPointer2 ){
1146     // TOTAL INTERNAL REFLECTION              << 693                      } else {
1147     if(sint2 >= 1.0)                          << 694                         theFacetNormal =
1148     {                                         << 695                             GetFacetNormal(OldMomentum,theGlobalNormal);
1149       swap = false;                           << 696                      }
1150                                               << 697                    }
1151       fStatus = TotalInternalReflection;      << 698 
1152       if(!surfaceRoughnessCriterionPass)      << 699                    G4double PdotN = OldMomentum * theFacetNormal;
1153         fStatus = LambertianReflection;       << 700                    NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1154       if(fModel == unified && fFinish != poli << 701                    G4double EdotN = OldPolarization * theFacetNormal;
1155         ChooseReflection();                   << 702 
1156       if(fStatus == LambertianReflection)     << 703                    G4ThreeVector A_trans, A_paral;
1157       {                                       << 704 
1158         DoReflection();                       << 705                    if (sint1 > 0.0 ) {
1159       }                                       << 706                       A_trans = OldMomentum.cross(theFacetNormal);
1160       else if(fStatus == BackScattering)      << 707                       A_trans = A_trans.unit();
1161       {                                       << 708                    } else {
1162         fNewMomentum     = -fOldMomentum;     << 709                       A_trans  = OldPolarization;
1163         fNewPolarization = -fOldPolarization; << 710                    }
1164       }                                       << 711                    A_paral   = NewMomentum.cross(A_trans);
1165       else                                    << 712                    A_paral   = A_paral.unit();
1166       {                                       << 713 
1167         fNewMomentum =                        << 714                    if(iTE>0&&iTM>0) {
1168           fOldMomentum - 2. * fOldMomentum *  << 715                      NewPolarization = 
1169         fNewPolarization = -fOldPolarization  << 716                            -OldPolarization + (2.*EdotN)*theFacetNormal;
1170                                               << 717                    } else if (iTE>0) {
1171       }                                       << 718                      NewPolarization = -A_trans;
1172     }                                         << 719                    } else if (iTM>0) {
1173     // NOT TIR                                << 720                      NewPolarization = -A_paral;
1174     else if(sint2 < 1.0)                      << 721                    }
1175     {                                         << 722 
1176       // Calculate amplitude for transmission << 723                 }
1177       if(cost1 > 0.0)                         << 
1178       {                                       << 
1179         cost2 = std::sqrt(1. - sint2 * sint2) << 
1180       }                                       << 
1181       else                                    << 
1182       {                                       << 
1183         cost2 = -std::sqrt(1. - sint2 * sint2 << 
1184       }                                       << 
1185                                               << 
1186       if(fSint1 > 0.0)                        << 
1187       {                                       << 
1188         A_trans = (fOldMomentum.cross(fFacetN << 
1189         E1_perp = fOldPolarization * A_trans; << 
1190         E1pp    = E1_perp * A_trans;          << 
1191         E1pl    = fOldPolarization - E1pp;    << 
1192         E1_parl = E1pl.mag();                 << 
1193       }                                       << 
1194       else                                    << 
1195       {                                       << 
1196         A_trans = fOldPolarization;           << 
1197         // Here we Follow Jackson's conventio << 
1198         // component = 1 in case of a ray per << 
1199         E1_perp = 0.0;                        << 
1200         E1_parl = 1.0;                        << 
1201       }                                       << 
1202                                               << 
1203       s1       = fRindex1 * cost1;            << 
1204       E2_perp  = 2. * s1 * E1_perp / (fRindex << 
1205       E2_parl  = 2. * s1 * E1_parl / (fRindex << 
1206       E2_total = E2_perp * E2_perp + E2_parl  << 
1207       s2       = fRindex2 * cost2 * E2_total; << 
1208                                               << 
1209       // D.Sawkey, 24 May 24                  << 
1210       // Transmittance has already been taken << 
1211       // For e.g. specular surfaces, the rati << 
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;                 << 
1219       else                                    << 
1220         transCoeff = 0.0;                     << 
1221                                               << 
1222       // NOT TIR: REFLECTION                  << 
1223       if(!G4BooleanRand(transCoeff))          << 
1224       {                                       << 
1225         swap    = false;                      << 
1226         fStatus = FresnelReflection;          << 
1227                                               << 
1228         if(!surfaceRoughnessCriterionPass)    << 
1229           fStatus = LambertianReflection;     << 
1230         if(fModel == unified && fFinish != po << 
1231           ChooseReflection();                 << 
1232         if(fStatus == LambertianReflection)   << 
1233         {                                     << 
1234           DoReflection();                     << 
1235         }                                     << 
1236         else if(fStatus == BackScattering)    << 
1237         {                                     << 
1238           fNewMomentum     = -fOldMomentum;   << 
1239           fNewPolarization = -fOldPolarizatio << 
1240         }                                     << 
1241         else                                  << 
1242         {                                     << 
1243           fNewMomentum =                      << 
1244             fOldMomentum - 2. * fOldMomentum  << 
1245           if(fSint1 > 0.0)                    << 
1246           {  // incident ray oblique          << 
1247             E2_parl  = fRindex2 * E2_parl / f << 
1248             E2_perp  = E2_perp - E1_perp;     << 
1249             E2_total = E2_perp * E2_perp + E2 << 
1250             A_paral  = (fNewMomentum.cross(A_ << 
1251             E2_abs   = std::sqrt(E2_total);   << 
1252             C_parl   = E2_parl / E2_abs;      << 
1253             C_perp   = E2_perp / E2_abs;      << 
1254                                                  724 
1255             fNewPolarization = C_parl * A_par << 725              }
1256           }                                   << 
1257           else                                << 
1258           {  // incident ray perpendicular    << 
1259             if(fRindex2 > fRindex1)           << 
1260             {                                 << 
1261               fNewPolarization = -fOldPolariz << 
1262             }                                 << 
1263             else                              << 
1264             {                                 << 
1265               fNewPolarization = fOldPolariza << 
1266             }                                 << 
1267           }                                   << 
1268         }                                     << 
1269       }                                       << 
1270       // NOT TIR: TRANSMISSION                << 
1271       else                                    << 
1272       {                                       << 
1273         inside  = !inside;                    << 
1274         through = true;                       << 
1275         fStatus = FresnelRefraction;          << 
1276                                               << 
1277         if(fSint1 > 0.0)                      << 
1278         {  // incident ray oblique            << 
1279           alpha        = cost1 - cost2 * (fRi << 
1280           fNewMomentum = (fOldMomentum + alph << 
1281           A_paral      = (fNewMomentum.cross( << 
1282           E2_abs       = std::sqrt(E2_total); << 
1283           C_parl       = E2_parl / E2_abs;    << 
1284           C_perp       = E2_perp / E2_abs;    << 
1285                                                  726 
1286           fNewPolarization = C_parl * A_paral << 727              OldMomentum = NewMomentum;
1287         }                                     << 728              OldPolarization = NewPolarization;
1288         else                                  << 
1289         {  // incident ray perpendicular      << 
1290           fNewMomentum     = fOldMomentum;    << 
1291           fNewPolarization = fOldPolarization << 
1292         }                                     << 
1293       }                                       << 
1294     }                                         << 
1295                                                  729 
1296     fOldMomentum     = fNewMomentum.unit();   << 730      }
1297     fOldPolarization = fNewPolarization.unit( << 
1298                                                  731 
1299     if(fStatus == FresnelRefraction)          << 732   } while (NewMomentum * theGlobalNormal < 0.0);
1300     {                                         << 733 }
1301       done = (fNewMomentum * fGlobalNormal <= << 
1302     }                                         << 
1303     else                                      << 
1304     {                                         << 
1305       done = (fNewMomentum * fGlobalNormal >= << 
1306     }                                         << 
1307     // Loop checking, 13-Aug-2015, Peter Gump << 
1308   } while(!done);                             << 
1309                                               << 
1310   if(inside && !swap)                         << 
1311   {                                           << 
1312     if(fFinish == polishedbackpainted || fFin << 
1313     {                                         << 
1314       G4double rand = G4UniformRand();        << 
1315       if(rand > fReflectivity + fTransmittanc << 
1316       {                                       << 
1317         DoAbsorption();                       << 
1318       }                                       << 
1319       else if(rand > fReflectivity)           << 
1320       {                                       << 
1321         fStatus          = Transmission;      << 
1322         fNewMomentum     = fOldMomentum;      << 
1323         fNewPolarization = fOldPolarization;  << 
1324       }                                       << 
1325       else                                    << 
1326       {                                       << 
1327         if(fStatus != FresnelRefraction)      << 
1328         {                                     << 
1329           fGlobalNormal = -fGlobalNormal;     << 
1330         }                                     << 
1331         else                                  << 
1332         {                                     << 
1333           swap = !swap;                       << 
1334           G4SwapPtr(fMaterial1, fMaterial2);  << 
1335           G4SwapObj(&fRindex1, &fRindex2);    << 
1336         }                                     << 
1337         if(fFinish == groundbackpainted)      << 
1338           fStatus = LambertianReflection;     << 
1339                                                  734 
1340         DoReflection();                       << 735 void G4OpBoundaryProcess::DielectricLUT()
                                                   >> 736 {
                                                   >> 737         G4int thetaIndex, phiIndex;
                                                   >> 738         G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
                                                   >> 739         G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
                                                   >> 740 
                                                   >> 741         theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) + 
                                                   >> 742                            (G4int(NoRINDEX)-G4int(groundbackpainted)));
                                                   >> 743 
                                                   >> 744         G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
                                                   >> 745         G4int phiIndexMax   = OpticalSurface->GetPhiIndexMax();
                                                   >> 746 
                                                   >> 747         do {
                                                   >> 748            if ( !G4BooleanRand(theReflectivity) ) // Not reflected, so Absorbed
                                                   >> 749               DoAbsorption();
                                                   >> 750            else {
                                                   >> 751               // Calculate Angle between Normal and Photon Momentum
                                                   >> 752               G4double anglePhotonToNormal = 
                                                   >> 753                                           OldMomentum.angle(-theGlobalNormal);
                                                   >> 754               // Round it to closest integer
                                                   >> 755               G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
                                                   >> 756 
                                                   >> 757               // Take random angles THETA and PHI, 
                                                   >> 758               // and see if below Probability - if not - Redo
                                                   >> 759               do {
                                                   >> 760                  thetaIndex = CLHEP::RandFlat::shootInt(thetaIndexMax-1);
                                                   >> 761                  phiIndex = CLHEP::RandFlat::shootInt(phiIndexMax-1);
                                                   >> 762                  // Find probability with the new indeces from LUT
                                                   >> 763                  AngularDistributionValue = OpticalSurface -> 
                                                   >> 764                    GetAngularDistributionValue(angleIncident,
                                                   >> 765                                                thetaIndex,
                                                   >> 766                                                phiIndex);
                                                   >> 767               } while ( !G4BooleanRand(AngularDistributionValue) );
                                                   >> 768 
                                                   >> 769               thetaRad = (-90 + 4*thetaIndex)*pi/180;
                                                   >> 770               phiRad = (-90 + 5*phiIndex)*pi/180;
                                                   >> 771               // Rotate Photon Momentum in Theta, then in Phi
                                                   >> 772               NewMomentum = -OldMomentum;
                                                   >> 773               PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
                                                   >> 774               if (PerpendicularVectorTheta.mag() > kCarTolerance ) {
                                                   >> 775                  PerpendicularVectorPhi = 
                                                   >> 776                                   PerpendicularVectorTheta.cross(NewMomentum);
                                                   >> 777               }
                                                   >> 778               else {
                                                   >> 779                  PerpendicularVectorTheta = NewMomentum.orthogonal();
                                                   >> 780                  PerpendicularVectorPhi =
                                                   >> 781                                   PerpendicularVectorTheta.cross(NewMomentum);
                                                   >> 782               }
                                                   >> 783               NewMomentum =
                                                   >> 784                  NewMomentum.rotate(anglePhotonToNormal-thetaRad,
                                                   >> 785                                     PerpendicularVectorTheta);
                                                   >> 786               NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
                                                   >> 787               // Rotate Polarization too:
                                                   >> 788               theFacetNormal = (NewMomentum - OldMomentum).unit();
                                                   >> 789               EdotN = OldPolarization * theFacetNormal;
                                                   >> 790               NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
                                                   >> 791            }
                                                   >> 792         } while (NewMomentum * theGlobalNormal <= 0.0);
                                                   >> 793 }
1341                                                  794 
1342         fGlobalNormal = -fGlobalNormal;       << 795 void G4OpBoundaryProcess::DielectricDielectric()
1343         fOldMomentum  = fNewMomentum;         << 796 {
                                                   >> 797   G4bool Inside = false;
                                                   >> 798   G4bool Swap = false;
1344                                                  799 
1345         goto leap;                            << 800   leap:
1346       }                                       << 801 
1347     }                                         << 802         G4bool Through = false;
1348   }                                           << 803   G4bool Done = false;
                                                   >> 804 
                                                   >> 805   do {
                                                   >> 806 
                                                   >> 807      if (Through) {
                                                   >> 808         Swap = !Swap;
                                                   >> 809         Through = false;
                                                   >> 810         theGlobalNormal = -theGlobalNormal;
                                                   >> 811         G4SwapPtr(Material1,Material2);
                                                   >> 812         G4SwapObj(&Rindex1,&Rindex2);
                                                   >> 813      }
                                                   >> 814   
                                                   >> 815      if ( theFinish == polished ) {
                                                   >> 816           theFacetNormal = theGlobalNormal;
                                                   >> 817      }
                                                   >> 818      else {
                                                   >> 819           theFacetNormal =
                                                   >> 820                        GetFacetNormal(OldMomentum,theGlobalNormal);
                                                   >> 821      }
                                                   >> 822 
                                                   >> 823      G4double PdotN = OldMomentum * theFacetNormal;
                                                   >> 824      G4double EdotN = OldPolarization * theFacetNormal;
                                                   >> 825 
                                                   >> 826      cost1 = - PdotN;
                                                   >> 827      if (std::abs(cost1) < 1.0-kCarTolerance){
                                                   >> 828         sint1 = std::sqrt(1.-cost1*cost1);
                                                   >> 829         sint2 = sint1*Rindex1/Rindex2;     // *** Snell's Law ***
                                                   >> 830      }
                                                   >> 831      else {
                                                   >> 832         sint1 = 0.0;
                                                   >> 833         sint2 = 0.0;
                                                   >> 834      }
                                                   >> 835 
                                                   >> 836      if (sint2 >= 1.0) {
                                                   >> 837 
                                                   >> 838         // Simulate total internal reflection
                                                   >> 839 
                                                   >> 840         if (Swap) Swap = !Swap;
                                                   >> 841 
                                                   >> 842               theStatus = TotalInternalReflection;
                                                   >> 843 
                                                   >> 844         if ( theModel == unified && theFinish != polished )
                                                   >> 845                  ChooseReflection();
                                                   >> 846 
                                                   >> 847         if ( theStatus == LambertianReflection ) {
                                                   >> 848      DoReflection();
                                                   >> 849         }
                                                   >> 850         else if ( theStatus == BackScattering ) {
                                                   >> 851      NewMomentum = -OldMomentum;
                                                   >> 852      NewPolarization = -OldPolarization;
                                                   >> 853         }
                                                   >> 854         else {
                                                   >> 855 
                                                   >> 856                  PdotN = OldMomentum * theFacetNormal;
                                                   >> 857      NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
                                                   >> 858      EdotN = OldPolarization * theFacetNormal;
                                                   >> 859      NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
                                                   >> 860 
                                                   >> 861         }
                                                   >> 862      }
                                                   >> 863      else if (sint2 < 1.0) {
                                                   >> 864 
                                                   >> 865         // Calculate amplitude for transmission (Q = P x N)
                                                   >> 866 
                                                   >> 867         if (cost1 > 0.0) {
                                                   >> 868            cost2 =  std::sqrt(1.-sint2*sint2);
                                                   >> 869         }
                                                   >> 870         else {
                                                   >> 871            cost2 = -std::sqrt(1.-sint2*sint2);
                                                   >> 872         }
                                                   >> 873 
                                                   >> 874         G4ThreeVector A_trans, A_paral, E1pp, E1pl;
                                                   >> 875         G4double E1_perp, E1_parl;
                                                   >> 876 
                                                   >> 877         if (sint1 > 0.0) {
                                                   >> 878            A_trans = OldMomentum.cross(theFacetNormal);
                                                   >> 879                  A_trans = A_trans.unit();
                                                   >> 880            E1_perp = OldPolarization * A_trans;
                                                   >> 881                  E1pp    = E1_perp * A_trans;
                                                   >> 882                  E1pl    = OldPolarization - E1pp;
                                                   >> 883                  E1_parl = E1pl.mag();
                                                   >> 884               }
                                                   >> 885         else {
                                                   >> 886            A_trans  = OldPolarization;
                                                   >> 887            // Here we Follow Jackson's conventions and we set the
                                                   >> 888            // parallel component = 1 in case of a ray perpendicular
                                                   >> 889            // to the surface
                                                   >> 890            E1_perp  = 0.0;
                                                   >> 891            E1_parl  = 1.0;
                                                   >> 892         }
                                                   >> 893 
                                                   >> 894               G4double s1 = Rindex1*cost1;
                                                   >> 895               G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
                                                   >> 896               G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
                                                   >> 897               G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
                                                   >> 898               G4double s2 = Rindex2*cost2*E2_total;
                                                   >> 899 
                                                   >> 900               G4double TransCoeff;
                                                   >> 901 
                                                   >> 902               if (theTransmittance > 0) TransCoeff = theTransmittance;
                                                   >> 903               else if (cost1 != 0.0) TransCoeff = s2/s1;
                                                   >> 904               else TransCoeff = 0.0;
                                                   >> 905 
                                                   >> 906         G4double E2_abs, C_parl, C_perp;
                                                   >> 907 
                                                   >> 908         if ( !G4BooleanRand(TransCoeff) ) {
                                                   >> 909 
                                                   >> 910            // Simulate reflection
                                                   >> 911 
                                                   >> 912                  if (Swap) Swap = !Swap;
                                                   >> 913 
                                                   >> 914      theStatus = FresnelReflection;
                                                   >> 915 
                                                   >> 916      if ( theModel == unified && theFinish != polished )
                                                   >> 917                  ChooseReflection();
                                                   >> 918 
                                                   >> 919      if ( theStatus == LambertianReflection ) {
                                                   >> 920         DoReflection();
                                                   >> 921      }
                                                   >> 922      else if ( theStatus == BackScattering ) {
                                                   >> 923         NewMomentum = -OldMomentum;
                                                   >> 924         NewPolarization = -OldPolarization;
                                                   >> 925      }
                                                   >> 926      else {
                                                   >> 927 
                                                   >> 928                     PdotN = OldMomentum * theFacetNormal;
                                                   >> 929               NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
                                                   >> 930 
                                                   >> 931               if (sint1 > 0.0) {   // incident ray oblique
                                                   >> 932 
                                                   >> 933            E2_parl   = Rindex2*E2_parl/Rindex1 - E1_parl;
                                                   >> 934            E2_perp   = E2_perp - E1_perp;
                                                   >> 935            E2_total  = E2_perp*E2_perp + E2_parl*E2_parl;
                                                   >> 936                        A_paral   = NewMomentum.cross(A_trans);
                                                   >> 937                        A_paral   = A_paral.unit();
                                                   >> 938            E2_abs    = std::sqrt(E2_total);
                                                   >> 939            C_parl    = E2_parl/E2_abs;
                                                   >> 940            C_perp    = E2_perp/E2_abs;
                                                   >> 941 
                                                   >> 942                        NewPolarization = C_parl*A_paral + C_perp*A_trans;
                                                   >> 943 
                                                   >> 944               }
                                                   >> 945 
                                                   >> 946               else {               // incident ray perpendicular
                                                   >> 947 
                                                   >> 948                  if (Rindex2 > Rindex1) {
                                                   >> 949               NewPolarization = - OldPolarization;
                                                   >> 950                  }
                                                   >> 951                  else {
                                                   >> 952                     NewPolarization =   OldPolarization;
                                                   >> 953                  }
                                                   >> 954 
                                                   >> 955               }
                                                   >> 956            }
                                                   >> 957         }
                                                   >> 958         else { // photon gets transmitted
                                                   >> 959 
                                                   >> 960            // Simulate transmission/refraction
                                                   >> 961 
                                                   >> 962      Inside = !Inside;
                                                   >> 963      Through = true;
                                                   >> 964      theStatus = FresnelRefraction;
                                                   >> 965 
                                                   >> 966            if (sint1 > 0.0) {      // incident ray oblique
                                                   >> 967 
                                                   >> 968         G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
                                                   >> 969         NewMomentum = OldMomentum + alpha*theFacetNormal;
                                                   >> 970         NewMomentum = NewMomentum.unit();
                                                   >> 971         PdotN = -cost2;
                                                   >> 972                     A_paral = NewMomentum.cross(A_trans);
                                                   >> 973                     A_paral = A_paral.unit();
                                                   >> 974         E2_abs     = std::sqrt(E2_total);
                                                   >> 975         C_parl     = E2_parl/E2_abs;
                                                   >> 976         C_perp     = E2_perp/E2_abs;
                                                   >> 977 
                                                   >> 978                     NewPolarization = C_parl*A_paral + C_perp*A_trans;
                                                   >> 979 
                                                   >> 980            }
                                                   >> 981            else {                  // incident ray perpendicular
                                                   >> 982 
                                                   >> 983         NewMomentum = OldMomentum;
                                                   >> 984         NewPolarization = OldPolarization;
                                                   >> 985 
                                                   >> 986            }
                                                   >> 987         }
                                                   >> 988      }
                                                   >> 989 
                                                   >> 990      OldMomentum = NewMomentum.unit();
                                                   >> 991      OldPolarization = NewPolarization.unit();
                                                   >> 992 
                                                   >> 993      if (theStatus == FresnelRefraction) {
                                                   >> 994         Done = (NewMomentum * theGlobalNormal <= 0.0);
                                                   >> 995      } 
                                                   >> 996      else {
                                                   >> 997         Done = (NewMomentum * theGlobalNormal >= 0.0);
                                                   >> 998      }
                                                   >> 999 
                                                   >> 1000   } while (!Done);
                                                   >> 1001 
                                                   >> 1002   if (Inside && !Swap) {
                                                   >> 1003           if( theFinish == polishedbackpainted ||
                                                   >> 1004               theFinish == groundbackpainted ) {
                                                   >> 1005 
                                                   >> 1006         if( !G4BooleanRand(theReflectivity) ) {
                                                   >> 1007     DoAbsorption();
                                                   >> 1008               }
                                                   >> 1009         else {
                                                   >> 1010     if (theStatus != FresnelRefraction ) {
                                                   >> 1011        theGlobalNormal = -theGlobalNormal;
                                                   >> 1012           }
                                                   >> 1013           else {
                                                   >> 1014        Swap = !Swap;
                                                   >> 1015        G4SwapPtr(Material1,Material2);
                                                   >> 1016        G4SwapObj(&Rindex1,&Rindex2);
                                                   >> 1017           }
                                                   >> 1018     if ( theFinish == groundbackpainted )
                                                   >> 1019           theStatus = LambertianReflection;
                                                   >> 1020 
                                                   >> 1021           DoReflection();
                                                   >> 1022 
                                                   >> 1023           theGlobalNormal = -theGlobalNormal;
                                                   >> 1024     OldMomentum = NewMomentum;
                                                   >> 1025 
                                                   >> 1026           goto leap;
                                                   >> 1027         }
                                                   >> 1028     }
                                                   >> 1029   }
1349 }                                                1030 }
1350                                                  1031 
1351 //....oooOO0OOooo........oooOO0OOooo........o << 1032 // GetMeanFreePath
1352 G4double G4OpBoundaryProcess::GetMeanFreePath << 1033 // ---------------
                                                   >> 1034 //
                                                   >> 1035 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
                                                   >> 1036                                               G4double ,
1353                                                  1037                                               G4ForceCondition* condition)
1354 {                                                1038 {
1355   *condition = Forced;                        << 1039   *condition = Forced;
1356   return DBL_MAX;                             << 1040 
                                                   >> 1041   return DBL_MAX;
1357 }                                                1042 }
1358                                                  1043 
1359 //....oooOO0OOooo........oooOO0OOooo........o << 1044 G4double G4OpBoundaryProcess::GetIncidentAngle() 
1360 G4double G4OpBoundaryProcess::GetIncidentAngl << 
1361 {                                                1045 {
1362   return pi - std::acos(fOldMomentum * fFacet << 1046     G4double PdotN = OldMomentum * theFacetNormal;
1363                         (fOldMomentum.mag() * << 1047     G4double magP= OldMomentum.mag();
                                                   >> 1048     G4double magN= theFacetNormal.mag();
                                                   >> 1049     G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
                                                   >> 1050 
                                                   >> 1051     return incidentangle;
1364 }                                                1052 }
1365                                                  1053 
1366 //....oooOO0OOooo........oooOO0OOooo........o << 
1367 G4double G4OpBoundaryProcess::GetReflectivity    1054 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1368                                                  1055                                               G4double E1_parl,
1369                                                  1056                                               G4double incidentangle,
1370                                               << 1057                                               G4double RealRindex,
1371                                               << 1058                                               G4double ImaginaryRindex)
1372 {                                                1059 {
1373   G4complex reflectivity, reflectivity_TE, re << 
1374   G4complex N1(fRindex1, 0.), N2(realRindex,  << 
1375   G4complex cosPhi;                           << 
1376                                                  1060 
1377   G4complex u(1., 0.);  // unit number 1      << 1061   G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
                                                   >> 1062   G4complex N(RealRindex, ImaginaryRindex);
                                                   >> 1063   G4complex CosPhi;
1378                                                  1064 
1379   G4complex numeratorTE;  // E1_perp=1 E1_par << 1065   G4complex u(1,0);           //unit number 1
1380   G4complex numeratorTM;  // E1_parl=1 E1_per << 1066 
                                                   >> 1067   G4complex numeratorTE;      // E1_perp=1 E1_parl=0 -> TE polarization
                                                   >> 1068   G4complex numeratorTM;      // E1_parl=1 E1_perp=0 -> TM polarization
1381   G4complex denominatorTE, denominatorTM;        1069   G4complex denominatorTE, denominatorTM;
1382   G4complex rTM, rTE;                            1070   G4complex rTM, rTE;
1383                                                  1071 
1384   G4MaterialPropertiesTable* MPT = fMaterial1 << 
1385   G4MaterialPropertyVector* ppR  = MPT->GetPr << 
1386   G4MaterialPropertyVector* ppI  = MPT->GetPr << 
1387   if(ppR && ppI)                              << 
1388   {                                           << 
1389     G4double rRindex = ppR->Value(fPhotonMome << 
1390     G4double iRindex = ppI->Value(fPhotonMome << 
1391     N1               = G4complex(rRindex, iRi << 
1392   }                                           << 
1393                                               << 
1394   // Following two equations, rTM and rTE, ar    1072   // Following two equations, rTM and rTE, are from: "Introduction To Modern
1395   // Optics" written by Fowles                   1073   // Optics" written by Fowles
1396   cosPhi = std::sqrt(u - ((std::sin(incidenta << 
1397                           (N1 * N1) / (N2 * N << 
1398                                                  1074 
1399   numeratorTE   = N1 * std::cos(incidentangle << 1075   CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));
1400   denominatorTE = N1 * std::cos(incidentangle << 1076 
1401   rTE           = numeratorTE / denominatorTE << 1077   numeratorTE   = std::cos(incidentangle) - N*CosPhi;
1402                                               << 1078   denominatorTE = std::cos(incidentangle) + N*CosPhi;
1403   numeratorTM   = N2 * std::cos(incidentangle << 1079   rTE = numeratorTE/denominatorTE;
1404   denominatorTM = N2 * std::cos(incidentangle << 1080 
1405   rTM           = numeratorTM / denominatorTM << 1081   numeratorTM   = N*std::cos(incidentangle) - CosPhi;
                                                   >> 1082   denominatorTM = N*std::cos(incidentangle) + CosPhi;
                                                   >> 1083   rTM = numeratorTM/denominatorTM;
1406                                                  1084 
1407   // This is my (PG) calculaton for reflectiv << 1085   // This is my calculaton for reflectivity on a metalic surface
1408   // depending on the fraction of TE and TM p    1086   // depending on the fraction of TE and TM polarization
1409   // when TE polarization, E1_parl=0 and E1_p    1087   // 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    1088   // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1411                                                  1089 
1412   reflectivity_TE = (rTE * conj(rTE)) * (E1_p << 1090   Reflectivity_TE =  (rTE*conj(rTE))*(E1_perp*E1_perp)
1413                     (E1_perp * E1_perp + E1_p << 1091                     / (E1_perp*E1_perp + E1_parl*E1_parl);
1414   reflectivity_TM = (rTM * conj(rTM)) * (E1_p << 1092   Reflectivity_TM =  (rTM*conj(rTM))*(E1_parl*E1_parl)
1415                     (E1_perp * E1_perp + E1_p << 1093                     / (E1_perp*E1_perp + E1_parl*E1_parl);
1416   reflectivity = reflectivity_TE + reflectivi << 1094   Reflectivity    = Reflectivity_TE + Reflectivity_TM;
1417                                               << 1095 
1418   do                                          << 1096   do {
1419   {                                           << 1097      if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1420     if(G4UniformRand() * real(reflectivity) > << 1098        {iTE = -1;}else{iTE = 1;}
1421     {                                         << 1099      if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1422       f_iTE = -1;                             << 1100        {iTM = -1;}else{iTM = 1;}
1423     }                                         << 1101   } while(iTE<0&&iTM<0);
1424     else                                      << 1102 
1425     {                                         << 1103   return real(Reflectivity);
1426       f_iTE = 1;                              << 
1427     }                                         << 
1428     if(G4UniformRand() * real(reflectivity) > << 
1429     {                                         << 
1430       f_iTM = -1;                             << 
1431     }                                         << 
1432     else                                      << 
1433     {                                         << 
1434       f_iTM = 1;                              << 
1435     }                                         << 
1436     // Loop checking, 13-Aug-2015, Peter Gump << 
1437   } while(f_iTE < 0 && f_iTM < 0);            << 
1438                                                  1104 
1439   return real(reflectivity);                  << 
1440 }                                                1105 }
1441                                                  1106 
1442 //....oooOO0OOooo........oooOO0OOooo........o << 
1443 void G4OpBoundaryProcess::CalculateReflectivi    1107 void G4OpBoundaryProcess::CalculateReflectivity()
1444 {                                                1108 {
1445   G4double realRindex = fRealRIndexMPV->Value << 1109   G4double RealRindex =
1446   G4double imaginaryRindex =                  << 1110            PropertyPointer1->Value(thePhotonMomentum);
1447     fImagRIndexMPV->Value(fPhotonMomentum, id << 1111   G4double ImaginaryRindex =
                                                   >> 1112            PropertyPointer2->Value(thePhotonMomentum);
1448                                                  1113 
1449   // calculate FacetNormal                       1114   // calculate FacetNormal
1450   if(fFinish == ground)                       << 1115   if ( theFinish == ground ) {
1451   {                                           << 1116      theFacetNormal =
1452     fFacetNormal = GetFacetNormal(fOldMomentu << 1117                GetFacetNormal(OldMomentum, theGlobalNormal);
1453   }                                           << 1118   } else {
1454   else                                        << 1119      theFacetNormal = theGlobalNormal;
1455   {                                           << 
1456     fFacetNormal = fGlobalNormal;             << 
1457   }                                              1120   }
1458                                                  1121 
1459   G4double cost1 = -fOldMomentum * fFacetNorm << 1122   G4double PdotN = OldMomentum * theFacetNormal;
1460   if(std::abs(cost1) < 1.0 - fCarTolerance)   << 1123   cost1 = -PdotN;
1461   {                                           << 1124 
1462     fSint1 = std::sqrt(1. - cost1 * cost1);   << 1125   if (std::abs(cost1) < 1.0 - kCarTolerance) {
1463   }                                           << 1126      sint1 = std::sqrt(1. - cost1*cost1);
1464   else                                        << 1127   } else {
1465   {                                           << 1128      sint1 = 0.0;
1466     fSint1 = 0.0;                             << 
1467   }                                              1129   }
1468                                                  1130 
1469   G4ThreeVector A_trans, A_paral, E1pp, E1pl;    1131   G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1470   G4double E1_perp, E1_parl;                     1132   G4double E1_perp, E1_parl;
1471                                                  1133 
1472   if(fSint1 > 0.0)                            << 1134   if (sint1 > 0.0 ) {
1473   {                                           << 1135      A_trans = OldMomentum.cross(theFacetNormal);
1474     A_trans = (fOldMomentum.cross(fFacetNorma << 1136      A_trans = A_trans.unit();
1475     E1_perp = fOldPolarization * A_trans;     << 1137      E1_perp = OldPolarization * A_trans;
1476     E1pp    = E1_perp * A_trans;              << 1138      E1pp    = E1_perp * A_trans;
1477     E1pl    = fOldPolarization - E1pp;        << 1139      E1pl    = OldPolarization - E1pp;
1478     E1_parl = E1pl.mag();                     << 1140      E1_parl = E1pl.mag();
1479   }                                           << 1141   }
1480   else                                        << 1142   else {
1481   {                                           << 1143      A_trans  = OldPolarization;
1482     A_trans = fOldPolarization;               << 1144      // Here we Follow Jackson's conventions and we set the
1483     // Here we Follow Jackson's conventions a << 1145      // parallel component = 1 in case of a ray perpendicular
1484     // component = 1 in case of a ray perpend << 1146      // to the surface
1485     E1_perp = 0.0;                            << 1147      E1_perp  = 0.0;
1486     E1_parl = 1.0;                            << 1148      E1_parl  = 1.0;
1487   }                                              1149   }
1488                                                  1150 
                                                   >> 1151   //calculate incident angle
1489   G4double incidentangle = GetIncidentAngle()    1152   G4double incidentangle = GetIncidentAngle();
1490                                                  1153 
1491   // calculate the reflectivity depending on  << 1154   //calculate the reflectivity depending on incident angle,
1492   // polarization and complex refractive      << 1155   //polarization and complex refractive
1493   fReflectivity = GetReflectivity(E1_perp, E1 << 
1494                                   imaginaryRi << 
1495 }                                             << 
1496                                                  1156 
1497 //....oooOO0OOooo........oooOO0OOooo........o << 1157   theReflectivity =
1498 G4bool G4OpBoundaryProcess::InvokeSD(const G4 << 1158              GetReflectivity(E1_perp, E1_parl, incidentangle,
1499 {                                             << 1159                                                  RealRindex, ImaginaryRindex);
1500   G4Step aStep = *pStep;                      << 
1501   aStep.AddTotalEnergyDeposit(fPhotonMomentum << 
1502                                               << 
1503   G4VSensitiveDetector* sd = aStep.GetPostSte << 
1504   if(sd != nullptr)                           << 
1505     return sd->Hit(&aStep);                   << 
1506   else                                        << 
1507     return false;                             << 
1508 }                                             << 
1509                                               << 
1510 //....oooOO0OOooo........oooOO0OOooo........o << 
1511 inline void G4OpBoundaryProcess::SetInvokeSD( << 
1512 {                                             << 
1513   fInvokeSD = flag;                           << 
1514   G4OpticalParameters::Instance()->SetBoundar << 
1515 }                                             << 
1516                                               << 
1517 //....oooOO0OOooo........oooOO0OOooo........o << 
1518 void G4OpBoundaryProcess::SetVerboseLevel(G4i << 
1519 {                                             << 
1520   verboseLevel = verbose;                     << 
1521   G4OpticalParameters::Instance()->SetBoundar << 
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 }                                                1160 }
1832                                                  1161