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