00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00027
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00072
00073 #include "G4ios.hh"
00074 #include "G4PhysicalConstants.hh"
00075 #include "G4OpProcessSubType.hh"
00076
00077 #include "G4OpBoundaryProcess.hh"
00078 #include "G4GeometryTolerance.hh"
00079
00081
00083
00085
00087
00088
00089
00090
00091
00093
00095
00096 G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
00097 G4ProcessType type)
00098 : G4VDiscreteProcess(processName, type)
00099 {
00100 if ( verboseLevel > 0) {
00101 G4cout << GetProcessName() << " is created " << G4endl;
00102 }
00103
00104 SetProcessSubType(fOpBoundary);
00105
00106 theStatus = Undefined;
00107 theModel = glisur;
00108 theFinish = polished;
00109 theReflectivity = 1.;
00110 theEfficiency = 0.;
00111 theTransmittance = 0.;
00112
00113 prob_sl = 0.;
00114 prob_ss = 0.;
00115 prob_bs = 0.;
00116
00117 PropertyPointer = NULL;
00118 PropertyPointer1 = NULL;
00119 PropertyPointer2 = NULL;
00120
00121 Material1 = NULL;
00122 Material2 = NULL;
00123
00124 OpticalSurface = NULL;
00125
00126 kCarTolerance = G4GeometryTolerance::GetInstance()
00127 ->GetSurfaceTolerance();
00128
00129 iTE = iTM = 0;
00130 thePhotonMomentum = 0.;
00131 Rindex1 = Rindex2 = cost1 = cost2 = sint1 = sint2 = 0.;
00132 }
00133
00134
00135
00136
00137
00139
00141
00142 G4OpBoundaryProcess::~G4OpBoundaryProcess(){}
00143
00145
00147
00148
00149
00150
00151 G4VParticleChange*
00152 G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00153 {
00154 theStatus = Undefined;
00155
00156 aParticleChange.Initialize(aTrack);
00157 aParticleChange.ProposeVelocity(aTrack.GetVelocity());
00158
00159 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
00160 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
00161
00162 if ( verboseLevel > 0 ) {
00163 G4cout << " Photon at Boundary! " << G4endl;
00164 G4VPhysicalVolume* thePrePV = pPreStepPoint->GetPhysicalVolume();
00165 G4VPhysicalVolume* thePostPV = pPostStepPoint->GetPhysicalVolume();
00166 if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
00167 if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
00168 }
00169
00170 if (pPostStepPoint->GetStepStatus() != fGeomBoundary){
00171 theStatus = NotAtBoundary;
00172 if ( verboseLevel > 0) BoundaryProcessVerbose();
00173 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00174 }
00175 if (aTrack.GetStepLength()<=kCarTolerance/2){
00176 theStatus = StepTooSmall;
00177 if ( verboseLevel > 0) BoundaryProcessVerbose();
00178 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00179 }
00180
00181 Material1 = pPreStepPoint -> GetMaterial();
00182 Material2 = pPostStepPoint -> GetMaterial();
00183
00184 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00185
00186 thePhotonMomentum = aParticle->GetTotalMomentum();
00187 OldMomentum = aParticle->GetMomentumDirection();
00188 OldPolarization = aParticle->GetPolarization();
00189
00190 if ( verboseLevel > 0 ) {
00191 G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
00192 G4cout << " Old Polarization: " << OldPolarization << G4endl;
00193 }
00194
00195 G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
00196
00197 G4Navigator* theNavigator =
00198 G4TransportationManager::GetTransportationManager()->
00199 GetNavigatorForTracking();
00200
00201 G4bool valid;
00202
00203
00204 theGlobalNormal =
00205 theNavigator->GetGlobalExitNormal(theGlobalPoint,&valid);
00206
00207 if (valid) {
00208 theGlobalNormal = -theGlobalNormal;
00209 }
00210 else
00211 {
00212 G4ExceptionDescription ed;
00213 ed << " G4OpBoundaryProcess/PostStepDoIt(): "
00214 << " The Navigator reports that it returned an invalid normal"
00215 << G4endl;
00216 G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
00217 EventMustBeAborted,ed,
00218 "Invalid Surface Normal - Geometry must return valid surface normal");
00219 }
00220
00221 if (OldMomentum * theGlobalNormal > 0.0) {
00222 #ifdef G4OPTICAL_DEBUG
00223 G4ExceptionDescription ed;
00224 ed << " G4OpBoundaryProcess/PostStepDoIt(): "
00225 << " theGlobalNormal points in a wrong direction. "
00226 << G4endl;
00227 ed << " The momentum of the photon arriving at interface (oldMomentum)"
00228 << " must exit the volume cross in the step. " << G4endl;
00229 ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
00230 ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
00231 ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
00232 ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
00233 ed << G4endl;
00234 G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
00235 EventMustBeAborted,
00236 ed,
00237 "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
00238 #else
00239 theGlobalNormal = -theGlobalNormal;
00240 #endif
00241 }
00242
00243 G4MaterialPropertiesTable* aMaterialPropertiesTable;
00244 G4MaterialPropertyVector* Rindex;
00245
00246 aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
00247 if (aMaterialPropertiesTable) {
00248 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
00249 }
00250 else {
00251 theStatus = NoRINDEX;
00252 if ( verboseLevel > 0) BoundaryProcessVerbose();
00253 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
00254 aParticleChange.ProposeTrackStatus(fStopAndKill);
00255 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00256 }
00257
00258 if (Rindex) {
00259 Rindex1 = Rindex->Value(thePhotonMomentum);
00260 }
00261 else {
00262 theStatus = NoRINDEX;
00263 if ( verboseLevel > 0) BoundaryProcessVerbose();
00264 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
00265 aParticleChange.ProposeTrackStatus(fStopAndKill);
00266 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00267 }
00268
00269 theReflectivity = 1.;
00270 theEfficiency = 0.;
00271 theTransmittance = 0.;
00272
00273 theModel = glisur;
00274 theFinish = polished;
00275
00276 G4SurfaceType type = dielectric_dielectric;
00277
00278 Rindex = NULL;
00279 OpticalSurface = NULL;
00280
00281 G4LogicalSurface* Surface = NULL;
00282
00283 Surface = G4LogicalBorderSurface::GetSurface
00284 (pPreStepPoint ->GetPhysicalVolume(),
00285 pPostStepPoint->GetPhysicalVolume());
00286
00287 if (Surface == NULL){
00288 G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
00289 ->GetMotherLogical() ==
00290 pPreStepPoint->GetPhysicalVolume()
00291 ->GetLogicalVolume());
00292 if(enteredDaughter){
00293 Surface = G4LogicalSkinSurface::GetSurface
00294 (pPostStepPoint->GetPhysicalVolume()->
00295 GetLogicalVolume());
00296 if(Surface == NULL)
00297 Surface = G4LogicalSkinSurface::GetSurface
00298 (pPreStepPoint->GetPhysicalVolume()->
00299 GetLogicalVolume());
00300 }
00301 else {
00302 Surface = G4LogicalSkinSurface::GetSurface
00303 (pPreStepPoint->GetPhysicalVolume()->
00304 GetLogicalVolume());
00305 if(Surface == NULL)
00306 Surface = G4LogicalSkinSurface::GetSurface
00307 (pPostStepPoint->GetPhysicalVolume()->
00308 GetLogicalVolume());
00309 }
00310 }
00311
00312 if (Surface) OpticalSurface =
00313 dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
00314
00315 if (OpticalSurface) {
00316
00317 type = OpticalSurface->GetType();
00318 theModel = OpticalSurface->GetModel();
00319 theFinish = OpticalSurface->GetFinish();
00320
00321 aMaterialPropertiesTable = OpticalSurface->
00322 GetMaterialPropertiesTable();
00323
00324 if (aMaterialPropertiesTable) {
00325
00326 if (theFinish == polishedbackpainted ||
00327 theFinish == groundbackpainted ) {
00328 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
00329 if (Rindex) {
00330 Rindex2 = Rindex->Value(thePhotonMomentum);
00331 }
00332 else {
00333 theStatus = NoRINDEX;
00334 if ( verboseLevel > 0) BoundaryProcessVerbose();
00335 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
00336 aParticleChange.ProposeTrackStatus(fStopAndKill);
00337 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00338 }
00339 }
00340
00341 PropertyPointer =
00342 aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
00343 PropertyPointer1 =
00344 aMaterialPropertiesTable->GetProperty("REALRINDEX");
00345 PropertyPointer2 =
00346 aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
00347
00348 iTE = 1;
00349 iTM = 1;
00350
00351 if (PropertyPointer) {
00352
00353 theReflectivity =
00354 PropertyPointer->Value(thePhotonMomentum);
00355
00356 } else if (PropertyPointer1 && PropertyPointer2) {
00357
00358 CalculateReflectivity();
00359
00360 }
00361
00362 PropertyPointer =
00363 aMaterialPropertiesTable->GetProperty("EFFICIENCY");
00364 if (PropertyPointer) {
00365 theEfficiency =
00366 PropertyPointer->Value(thePhotonMomentum);
00367 }
00368
00369 PropertyPointer =
00370 aMaterialPropertiesTable->GetProperty("TRANSMITTANCE");
00371 if (PropertyPointer) {
00372 theTransmittance =
00373 PropertyPointer->Value(thePhotonMomentum);
00374 }
00375
00376 if ( theModel == unified ) {
00377 PropertyPointer =
00378 aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
00379 if (PropertyPointer) {
00380 prob_sl =
00381 PropertyPointer->Value(thePhotonMomentum);
00382 } else {
00383 prob_sl = 0.0;
00384 }
00385
00386 PropertyPointer =
00387 aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
00388 if (PropertyPointer) {
00389 prob_ss =
00390 PropertyPointer->Value(thePhotonMomentum);
00391 } else {
00392 prob_ss = 0.0;
00393 }
00394
00395 PropertyPointer =
00396 aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
00397 if (PropertyPointer) {
00398 prob_bs =
00399 PropertyPointer->Value(thePhotonMomentum);
00400 } else {
00401 prob_bs = 0.0;
00402 }
00403 }
00404 }
00405 else if (theFinish == polishedbackpainted ||
00406 theFinish == groundbackpainted ) {
00407 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
00408 aParticleChange.ProposeTrackStatus(fStopAndKill);
00409 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00410 }
00411 }
00412
00413 if (type == dielectric_dielectric ) {
00414 if (theFinish == polished || theFinish == ground ) {
00415
00416 if (Material1 == Material2){
00417 theStatus = SameMaterial;
00418 if ( verboseLevel > 0) BoundaryProcessVerbose();
00419 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00420 }
00421 aMaterialPropertiesTable =
00422 Material2->GetMaterialPropertiesTable();
00423 if (aMaterialPropertiesTable)
00424 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
00425 if (Rindex) {
00426 Rindex2 = Rindex->Value(thePhotonMomentum);
00427 }
00428 else {
00429 theStatus = NoRINDEX;
00430 if ( verboseLevel > 0) BoundaryProcessVerbose();
00431 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
00432 aParticleChange.ProposeTrackStatus(fStopAndKill);
00433 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00434 }
00435 }
00436 }
00437
00438 if (type == dielectric_metal) {
00439
00440 DielectricMetal();
00441
00442
00443
00444
00445
00446
00447
00448 }
00449 else if (type == dielectric_LUT) {
00450
00451 DielectricLUT();
00452
00453 }
00454 else if (type == dielectric_dielectric) {
00455
00456 if ( theFinish == polishedbackpainted ||
00457 theFinish == groundbackpainted ) {
00458 DielectricDielectric();
00459 }
00460 else {
00461 if ( !G4BooleanRand(theReflectivity) ) {
00462 DoAbsorption();
00463 }
00464 else {
00465 if ( theFinish == polishedfrontpainted ) {
00466 DoReflection();
00467 }
00468 else if ( theFinish == groundfrontpainted ) {
00469 theStatus = LambertianReflection;
00470 DoReflection();
00471 }
00472 else {
00473 DielectricDielectric();
00474 }
00475 }
00476 }
00477 }
00478 else {
00479
00480 G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
00481 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00482
00483 }
00484
00485 NewMomentum = NewMomentum.unit();
00486 NewPolarization = NewPolarization.unit();
00487
00488 if ( verboseLevel > 0) {
00489 G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
00490 G4cout << " New Polarization: " << NewPolarization << G4endl;
00491 BoundaryProcessVerbose();
00492 }
00493
00494 aParticleChange.ProposeMomentumDirection(NewMomentum);
00495 aParticleChange.ProposePolarization(NewPolarization);
00496
00497 if ( theStatus == FresnelRefraction ) {
00498 G4MaterialPropertyVector* groupvel =
00499 Material2->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
00500 G4double finalVelocity = groupvel->Value(thePhotonMomentum);
00501 aParticleChange.ProposeVelocity(finalVelocity);
00502 }
00503
00504 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00505 }
00506
00507 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
00508 {
00509 if ( theStatus == Undefined )
00510 G4cout << " *** Undefined *** " << G4endl;
00511 if ( theStatus == FresnelRefraction )
00512 G4cout << " *** FresnelRefraction *** " << G4endl;
00513 if ( theStatus == FresnelReflection )
00514 G4cout << " *** FresnelReflection *** " << G4endl;
00515 if ( theStatus == TotalInternalReflection )
00516 G4cout << " *** TotalInternalReflection *** " << G4endl;
00517 if ( theStatus == LambertianReflection )
00518 G4cout << " *** LambertianReflection *** " << G4endl;
00519 if ( theStatus == LobeReflection )
00520 G4cout << " *** LobeReflection *** " << G4endl;
00521 if ( theStatus == SpikeReflection )
00522 G4cout << " *** SpikeReflection *** " << G4endl;
00523 if ( theStatus == BackScattering )
00524 G4cout << " *** BackScattering *** " << G4endl;
00525 if ( theStatus == PolishedLumirrorAirReflection )
00526 G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
00527 if ( theStatus == PolishedLumirrorGlueReflection )
00528 G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
00529 if ( theStatus == PolishedAirReflection )
00530 G4cout << " *** PolishedAirReflection *** " << G4endl;
00531 if ( theStatus == PolishedTeflonAirReflection )
00532 G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
00533 if ( theStatus == PolishedTiOAirReflection )
00534 G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
00535 if ( theStatus == PolishedTyvekAirReflection )
00536 G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
00537 if ( theStatus == PolishedVM2000AirReflection )
00538 G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
00539 if ( theStatus == PolishedVM2000GlueReflection )
00540 G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
00541 if ( theStatus == EtchedLumirrorAirReflection )
00542 G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
00543 if ( theStatus == EtchedLumirrorGlueReflection )
00544 G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
00545 if ( theStatus == EtchedAirReflection )
00546 G4cout << " *** EtchedAirReflection *** " << G4endl;
00547 if ( theStatus == EtchedTeflonAirReflection )
00548 G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
00549 if ( theStatus == EtchedTiOAirReflection )
00550 G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
00551 if ( theStatus == EtchedTyvekAirReflection )
00552 G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
00553 if ( theStatus == EtchedVM2000AirReflection )
00554 G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
00555 if ( theStatus == EtchedVM2000GlueReflection )
00556 G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
00557 if ( theStatus == GroundLumirrorAirReflection )
00558 G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
00559 if ( theStatus == GroundLumirrorGlueReflection )
00560 G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
00561 if ( theStatus == GroundAirReflection )
00562 G4cout << " *** GroundAirReflection *** " << G4endl;
00563 if ( theStatus == GroundTeflonAirReflection )
00564 G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
00565 if ( theStatus == GroundTiOAirReflection )
00566 G4cout << " *** GroundTiOAirReflection *** " << G4endl;
00567 if ( theStatus == GroundTyvekAirReflection )
00568 G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
00569 if ( theStatus == GroundVM2000AirReflection )
00570 G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
00571 if ( theStatus == GroundVM2000GlueReflection )
00572 G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
00573 if ( theStatus == Absorption )
00574 G4cout << " *** Absorption *** " << G4endl;
00575 if ( theStatus == Detection )
00576 G4cout << " *** Detection *** " << G4endl;
00577 if ( theStatus == NotAtBoundary )
00578 G4cout << " *** NotAtBoundary *** " << G4endl;
00579 if ( theStatus == SameMaterial )
00580 G4cout << " *** SameMaterial *** " << G4endl;
00581 if ( theStatus == StepTooSmall )
00582 G4cout << " *** StepTooSmall *** " << G4endl;
00583 if ( theStatus == NoRINDEX )
00584 G4cout << " *** NoRINDEX *** " << G4endl;
00585 }
00586
00587 G4ThreeVector
00588 G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
00589 const G4ThreeVector& Normal ) const
00590 {
00591 G4ThreeVector FacetNormal;
00592
00593 if (theModel == unified || theModel == LUT) {
00594
00595
00596
00597
00598
00599
00600
00601 G4double alpha;
00602
00603 G4double sigma_alpha = 0.0;
00604 if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
00605
00606 if (sigma_alpha == 0.0) return FacetNormal = Normal;
00607
00608 G4double f_max = std::min(1.0,4.*sigma_alpha);
00609
00610 do {
00611 do {
00612 alpha = G4RandGauss::shoot(0.0,sigma_alpha);
00613 } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
00614
00615 G4double phi = G4UniformRand()*twopi;
00616
00617 G4double SinAlpha = std::sin(alpha);
00618 G4double CosAlpha = std::cos(alpha);
00619 G4double SinPhi = std::sin(phi);
00620 G4double CosPhi = std::cos(phi);
00621
00622 G4double unit_x = SinAlpha * CosPhi;
00623 G4double unit_y = SinAlpha * SinPhi;
00624 G4double unit_z = CosAlpha;
00625
00626 FacetNormal.setX(unit_x);
00627 FacetNormal.setY(unit_y);
00628 FacetNormal.setZ(unit_z);
00629
00630 G4ThreeVector tmpNormal = Normal;
00631
00632 FacetNormal.rotateUz(tmpNormal);
00633 } while (Momentum * FacetNormal >= 0.0);
00634 }
00635 else {
00636
00637 G4double polish = 1.0;
00638 if (OpticalSurface) polish = OpticalSurface->GetPolish();
00639
00640 if (polish < 1.0) {
00641 do {
00642 G4ThreeVector smear;
00643 do {
00644 smear.setX(2.*G4UniformRand()-1.0);
00645 smear.setY(2.*G4UniformRand()-1.0);
00646 smear.setZ(2.*G4UniformRand()-1.0);
00647 } while (smear.mag()>1.0);
00648 smear = (1.-polish) * smear;
00649 FacetNormal = Normal + smear;
00650 } while (Momentum * FacetNormal >= 0.0);
00651 FacetNormal = FacetNormal.unit();
00652 }
00653 else {
00654 FacetNormal = Normal;
00655 }
00656 }
00657 return FacetNormal;
00658 }
00659
00660 void G4OpBoundaryProcess::DielectricMetal()
00661 {
00662 G4int n = 0;
00663
00664 do {
00665
00666 n++;
00667
00668 if( !G4BooleanRand(theReflectivity) && n == 1 ) {
00669
00670
00671
00672
00673 DoAbsorption();
00674
00675 break;
00676
00677 }
00678 else {
00679
00680 if (PropertyPointer1 && PropertyPointer2) {
00681 if ( n > 1 ) {
00682 CalculateReflectivity();
00683 if ( !G4BooleanRand(theReflectivity) ) {
00684 DoAbsorption();
00685 break;
00686 }
00687 }
00688 }
00689
00690 if ( theModel == glisur || theFinish == polished ) {
00691
00692 DoReflection();
00693
00694 } else {
00695
00696 if ( n == 1 ) ChooseReflection();
00697
00698 if ( theStatus == LambertianReflection ) {
00699 DoReflection();
00700 }
00701 else if ( theStatus == BackScattering ) {
00702 NewMomentum = -OldMomentum;
00703 NewPolarization = -OldPolarization;
00704 }
00705 else {
00706
00707 if(theStatus==LobeReflection){
00708 if ( PropertyPointer1 && PropertyPointer2 ){
00709 } else {
00710 theFacetNormal =
00711 GetFacetNormal(OldMomentum,theGlobalNormal);
00712 }
00713 }
00714
00715 G4double PdotN = OldMomentum * theFacetNormal;
00716 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
00717 G4double EdotN = OldPolarization * theFacetNormal;
00718
00719 G4ThreeVector A_trans, A_paral;
00720
00721 if (sint1 > 0.0 ) {
00722 A_trans = OldMomentum.cross(theFacetNormal);
00723 A_trans = A_trans.unit();
00724 } else {
00725 A_trans = OldPolarization;
00726 }
00727 A_paral = NewMomentum.cross(A_trans);
00728 A_paral = A_paral.unit();
00729
00730 if(iTE>0&&iTM>0) {
00731 NewPolarization =
00732 -OldPolarization + (2.*EdotN)*theFacetNormal;
00733 } else if (iTE>0) {
00734 NewPolarization = -A_trans;
00735 } else if (iTM>0) {
00736 NewPolarization = -A_paral;
00737 }
00738
00739 }
00740
00741 }
00742
00743 OldMomentum = NewMomentum;
00744 OldPolarization = NewPolarization;
00745
00746 }
00747
00748 } while (NewMomentum * theGlobalNormal < 0.0);
00749 }
00750
00751 void G4OpBoundaryProcess::DielectricLUT()
00752 {
00753 G4int thetaIndex, phiIndex;
00754 G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
00755 G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
00756
00757 theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) +
00758 (G4int(NoRINDEX)-G4int(groundbackpainted)));
00759
00760 G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
00761 G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
00762
00763 do {
00764 if ( !G4BooleanRand(theReflectivity) )
00765 DoAbsorption();
00766 else {
00767
00768 G4double anglePhotonToNormal =
00769 OldMomentum.angle(-theGlobalNormal);
00770
00771 G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
00772
00773
00774
00775 do {
00776 thetaIndex = CLHEP::RandFlat::shootInt(thetaIndexMax-1);
00777 phiIndex = CLHEP::RandFlat::shootInt(phiIndexMax-1);
00778
00779 AngularDistributionValue = OpticalSurface ->
00780 GetAngularDistributionValue(angleIncident,
00781 thetaIndex,
00782 phiIndex);
00783 } while ( !G4BooleanRand(AngularDistributionValue) );
00784
00785 thetaRad = (-90 + 4*thetaIndex)*pi/180;
00786 phiRad = (-90 + 5*phiIndex)*pi/180;
00787
00788 NewMomentum = -OldMomentum;
00789 PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
00790 if (PerpendicularVectorTheta.mag() > kCarTolerance ) {
00791 PerpendicularVectorPhi =
00792 PerpendicularVectorTheta.cross(NewMomentum);
00793 }
00794 else {
00795 PerpendicularVectorTheta = NewMomentum.orthogonal();
00796 PerpendicularVectorPhi =
00797 PerpendicularVectorTheta.cross(NewMomentum);
00798 }
00799 NewMomentum =
00800 NewMomentum.rotate(anglePhotonToNormal-thetaRad,
00801 PerpendicularVectorTheta);
00802 NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
00803
00804 theFacetNormal = (NewMomentum - OldMomentum).unit();
00805 EdotN = OldPolarization * theFacetNormal;
00806 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
00807 }
00808 } while (NewMomentum * theGlobalNormal <= 0.0);
00809 }
00810
00811 void G4OpBoundaryProcess::DielectricDielectric()
00812 {
00813 G4bool Inside = false;
00814 G4bool Swap = false;
00815
00816 leap:
00817
00818 G4bool Through = false;
00819 G4bool Done = false;
00820
00821 do {
00822
00823 if (Through) {
00824 Swap = !Swap;
00825 Through = false;
00826 theGlobalNormal = -theGlobalNormal;
00827 G4SwapPtr(Material1,Material2);
00828 G4SwapObj(&Rindex1,&Rindex2);
00829 }
00830
00831 if ( theFinish == polished ) {
00832 theFacetNormal = theGlobalNormal;
00833 }
00834 else {
00835 theFacetNormal =
00836 GetFacetNormal(OldMomentum,theGlobalNormal);
00837 }
00838
00839 G4double PdotN = OldMomentum * theFacetNormal;
00840 G4double EdotN = OldPolarization * theFacetNormal;
00841
00842 cost1 = - PdotN;
00843 if (std::abs(cost1) < 1.0-kCarTolerance){
00844 sint1 = std::sqrt(1.-cost1*cost1);
00845 sint2 = sint1*Rindex1/Rindex2;
00846 }
00847 else {
00848 sint1 = 0.0;
00849 sint2 = 0.0;
00850 }
00851
00852 if (sint2 >= 1.0) {
00853
00854
00855
00856 if (Swap) Swap = !Swap;
00857
00858 theStatus = TotalInternalReflection;
00859
00860 if ( theModel == unified && theFinish != polished )
00861 ChooseReflection();
00862
00863 if ( theStatus == LambertianReflection ) {
00864 DoReflection();
00865 }
00866 else if ( theStatus == BackScattering ) {
00867 NewMomentum = -OldMomentum;
00868 NewPolarization = -OldPolarization;
00869 }
00870 else {
00871
00872 PdotN = OldMomentum * theFacetNormal;
00873 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
00874 EdotN = OldPolarization * theFacetNormal;
00875 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
00876
00877 }
00878 }
00879 else if (sint2 < 1.0) {
00880
00881
00882
00883 if (cost1 > 0.0) {
00884 cost2 = std::sqrt(1.-sint2*sint2);
00885 }
00886 else {
00887 cost2 = -std::sqrt(1.-sint2*sint2);
00888 }
00889
00890 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
00891 G4double E1_perp, E1_parl;
00892
00893 if (sint1 > 0.0) {
00894 A_trans = OldMomentum.cross(theFacetNormal);
00895 A_trans = A_trans.unit();
00896 E1_perp = OldPolarization * A_trans;
00897 E1pp = E1_perp * A_trans;
00898 E1pl = OldPolarization - E1pp;
00899 E1_parl = E1pl.mag();
00900 }
00901 else {
00902 A_trans = OldPolarization;
00903
00904
00905
00906 E1_perp = 0.0;
00907 E1_parl = 1.0;
00908 }
00909
00910 G4double s1 = Rindex1*cost1;
00911 G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
00912 G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
00913 G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
00914 G4double s2 = Rindex2*cost2*E2_total;
00915
00916 G4double TransCoeff;
00917
00918 if (theTransmittance > 0) TransCoeff = theTransmittance;
00919 else if (cost1 != 0.0) TransCoeff = s2/s1;
00920 else TransCoeff = 0.0;
00921
00922 G4double E2_abs, C_parl, C_perp;
00923
00924 if ( !G4BooleanRand(TransCoeff) ) {
00925
00926
00927
00928 if (Swap) Swap = !Swap;
00929
00930 theStatus = FresnelReflection;
00931
00932 if ( theModel == unified && theFinish != polished )
00933 ChooseReflection();
00934
00935 if ( theStatus == LambertianReflection ) {
00936 DoReflection();
00937 }
00938 else if ( theStatus == BackScattering ) {
00939 NewMomentum = -OldMomentum;
00940 NewPolarization = -OldPolarization;
00941 }
00942 else {
00943
00944 PdotN = OldMomentum * theFacetNormal;
00945 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
00946
00947 if (sint1 > 0.0) {
00948
00949 E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
00950 E2_perp = E2_perp - E1_perp;
00951 E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
00952 A_paral = NewMomentum.cross(A_trans);
00953 A_paral = A_paral.unit();
00954 E2_abs = std::sqrt(E2_total);
00955 C_parl = E2_parl/E2_abs;
00956 C_perp = E2_perp/E2_abs;
00957
00958 NewPolarization = C_parl*A_paral + C_perp*A_trans;
00959
00960 }
00961
00962 else {
00963
00964 if (Rindex2 > Rindex1) {
00965 NewPolarization = - OldPolarization;
00966 }
00967 else {
00968 NewPolarization = OldPolarization;
00969 }
00970
00971 }
00972 }
00973 }
00974 else {
00975
00976
00977
00978 Inside = !Inside;
00979 Through = true;
00980 theStatus = FresnelRefraction;
00981
00982 if (sint1 > 0.0) {
00983
00984 G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
00985 NewMomentum = OldMomentum + alpha*theFacetNormal;
00986 NewMomentum = NewMomentum.unit();
00987 PdotN = -cost2;
00988 A_paral = NewMomentum.cross(A_trans);
00989 A_paral = A_paral.unit();
00990 E2_abs = std::sqrt(E2_total);
00991 C_parl = E2_parl/E2_abs;
00992 C_perp = E2_perp/E2_abs;
00993
00994 NewPolarization = C_parl*A_paral + C_perp*A_trans;
00995
00996 }
00997 else {
00998
00999 NewMomentum = OldMomentum;
01000 NewPolarization = OldPolarization;
01001
01002 }
01003 }
01004 }
01005
01006 OldMomentum = NewMomentum.unit();
01007 OldPolarization = NewPolarization.unit();
01008
01009 if (theStatus == FresnelRefraction) {
01010 Done = (NewMomentum * theGlobalNormal <= 0.0);
01011 }
01012 else {
01013 Done = (NewMomentum * theGlobalNormal >= 0.0);
01014 }
01015
01016 } while (!Done);
01017
01018 if (Inside && !Swap) {
01019 if( theFinish == polishedbackpainted ||
01020 theFinish == groundbackpainted ) {
01021
01022 if( !G4BooleanRand(theReflectivity) ) {
01023 DoAbsorption();
01024 }
01025 else {
01026 if (theStatus != FresnelRefraction ) {
01027 theGlobalNormal = -theGlobalNormal;
01028 }
01029 else {
01030 Swap = !Swap;
01031 G4SwapPtr(Material1,Material2);
01032 G4SwapObj(&Rindex1,&Rindex2);
01033 }
01034 if ( theFinish == groundbackpainted )
01035 theStatus = LambertianReflection;
01036
01037 DoReflection();
01038
01039 theGlobalNormal = -theGlobalNormal;
01040 OldMomentum = NewMomentum;
01041
01042 goto leap;
01043 }
01044 }
01045 }
01046 }
01047
01048
01049
01050
01051 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
01052 G4double ,
01053 G4ForceCondition* condition)
01054 {
01055 *condition = Forced;
01056
01057 return DBL_MAX;
01058 }
01059
01060 G4double G4OpBoundaryProcess::GetIncidentAngle()
01061 {
01062 G4double PdotN = OldMomentum * theFacetNormal;
01063 G4double magP= OldMomentum.mag();
01064 G4double magN= theFacetNormal.mag();
01065 G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
01066
01067 return incidentangle;
01068 }
01069
01070 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
01071 G4double E1_parl,
01072 G4double incidentangle,
01073 G4double RealRindex,
01074 G4double ImaginaryRindex)
01075 {
01076
01077 G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
01078 G4complex N(RealRindex, ImaginaryRindex);
01079 G4complex CosPhi;
01080
01081 G4complex u(1,0);
01082
01083 G4complex numeratorTE;
01084 G4complex numeratorTM;
01085 G4complex denominatorTE, denominatorTM;
01086 G4complex rTM, rTE;
01087
01088
01089
01090
01091 CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));
01092
01093 numeratorTE = std::cos(incidentangle) - N*CosPhi;
01094 denominatorTE = std::cos(incidentangle) + N*CosPhi;
01095 rTE = numeratorTE/denominatorTE;
01096
01097 numeratorTM = N*std::cos(incidentangle) - CosPhi;
01098 denominatorTM = N*std::cos(incidentangle) + CosPhi;
01099 rTM = numeratorTM/denominatorTM;
01100
01101
01102
01103
01104
01105
01106 Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
01107 / (E1_perp*E1_perp + E1_parl*E1_parl);
01108 Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
01109 / (E1_perp*E1_perp + E1_parl*E1_parl);
01110 Reflectivity = Reflectivity_TE + Reflectivity_TM;
01111
01112 do {
01113 if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
01114 {iTE = -1;}else{iTE = 1;}
01115 if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
01116 {iTM = -1;}else{iTM = 1;}
01117 } while(iTE<0&&iTM<0);
01118
01119 return real(Reflectivity);
01120
01121 }
01122
01123 void G4OpBoundaryProcess::CalculateReflectivity()
01124 {
01125 G4double RealRindex =
01126 PropertyPointer1->Value(thePhotonMomentum);
01127 G4double ImaginaryRindex =
01128 PropertyPointer2->Value(thePhotonMomentum);
01129
01130
01131 if ( theFinish == ground ) {
01132 theFacetNormal =
01133 GetFacetNormal(OldMomentum, theGlobalNormal);
01134 } else {
01135 theFacetNormal = theGlobalNormal;
01136 }
01137
01138 G4double PdotN = OldMomentum * theFacetNormal;
01139 cost1 = -PdotN;
01140
01141 if (std::abs(cost1) < 1.0 - kCarTolerance) {
01142 sint1 = std::sqrt(1. - cost1*cost1);
01143 } else {
01144 sint1 = 0.0;
01145 }
01146
01147 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
01148 G4double E1_perp, E1_parl;
01149
01150 if (sint1 > 0.0 ) {
01151 A_trans = OldMomentum.cross(theFacetNormal);
01152 A_trans = A_trans.unit();
01153 E1_perp = OldPolarization * A_trans;
01154 E1pp = E1_perp * A_trans;
01155 E1pl = OldPolarization - E1pp;
01156 E1_parl = E1pl.mag();
01157 }
01158 else {
01159 A_trans = OldPolarization;
01160
01161
01162
01163 E1_perp = 0.0;
01164 E1_parl = 1.0;
01165 }
01166
01167
01168 G4double incidentangle = GetIncidentAngle();
01169
01170
01171
01172
01173 theReflectivity =
01174 GetReflectivity(E1_perp, E1_parl, incidentangle,
01175 RealRindex, ImaginaryRindex);
01176 }