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