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