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