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