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