G4OpBoundaryProcess.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00027 // Optical Photon Boundary Process Class Implementation
00029 //
00030 // File:        G4OpBoundaryProcess.cc
00031 // Description: Discrete Process -- reflection/refraction at
00032 //                                  optical interfaces
00033 // Version:     1.1
00034 // Created:     1997-06-18
00035 // Modified:    1998-05-25 - Correct parallel component of polarization
00036 //                           (thanks to: Stefano Magni + Giovanni Pieri)
00037 //              1998-05-28 - NULL Rindex pointer before reuse
00038 //                           (thanks to: Stefano Magni)
00039 //              1998-06-11 - delete *sint1 in oblique reflection
00040 //                           (thanks to: Giovanni Pieri)
00041 //              1998-06-19 - move from GetLocalExitNormal() to the new 
00042 //                           method: GetLocalExitNormal(&valid) to get
00043 //                           the surface normal in all cases
00044 //              1998-11-07 - NULL OpticalSurface pointer before use
00045 //                           comparison not sharp for: std::abs(cost1) < 1.0
00046 //                           remove sin1, sin2 in lines 556,567
00047 //                           (thanks to Stefano Magni)
00048 //              1999-10-10 - Accommodate changes done in DoAbsorption by
00049 //                           changing logic in DielectricMetal
00050 //              2001-10-18 - avoid Linux (gcc-2.95.2) warning about variables
00051 //                           might be used uninitialized in this function
00052 //                           moved E2_perp, E2_parl and E2_total out of 'if'
00053 //              2003-11-27 - Modified line 168-9 to reflect changes made to
00054 //                           G4OpticalSurface class ( by Fan Lei)
00055 //              2004-02-02 - Set theStatus = Undefined at start of DoIt
00056 //              2005-07-28 - add G4ProcessType to constructor
00057 //              2006-11-04 - add capability of calculating the reflectivity
00058 //                           off a metal surface by way of a complex index 
00059 //                           of refraction - Thanks to Sehwook Lee and John 
00060 //                           Hauptman (Dept. of Physics - Iowa State Univ.)
00061 //              2009-11-10 - add capability of simulating surface reflections
00062 //                           with Look-Up-Tables (LUT) containing measured
00063 //                           optical reflectance for a variety of surface
00064 //                           treatments - Thanks to Martin Janecek and
00065 //                           William Moses (Lawrence Berkeley National Lab.)
00066 //
00067 // Author:      Peter Gumplinger
00068 //              adopted from work by Werner Keil - April 2/96
00069 // mail:        gum@triumf.ca
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 // Class Implementation
00083 
00085         // Operators
00087 
00088 // G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
00089 // {
00090 // }
00091 
00093         // Constructors
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 // G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
00135 // {
00136 // }
00137 
00139         // Destructors
00141 
00142 G4OpBoundaryProcess::~G4OpBoundaryProcess(){}
00143 
00145         // Methods
00147 
00148 // PostStepDoIt
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         //  Use the new method for Exit Normal in global coordinates,
00203         //    which provides the normal more reliably. 
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,  // Or JustWarning to see if it happens repeatedbly on one ray
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           // Uncomment the following lines if you wish to have 
00443           //         Transmission instead of Absorption
00444           // if (theStatus == Absorption) {
00445           //    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
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         /* This function code alpha to a random value taken from the
00596            distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
00597            for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
00598            is a gaussian distribution with mean 0 and standard deviation
00599            sigma_alpha.  */
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              // Comment out DoAbsorption and uncomment theStatus = Absorption;
00671              // if you wish to have Transmission instead of Absorption
00672 
00673              DoAbsorption();
00674              // theStatus = Absorption;
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) ) // Not reflected, so Absorbed
00765               DoAbsorption();
00766            else {
00767               // Calculate Angle between Normal and Photon Momentum
00768               G4double anglePhotonToNormal = 
00769                                           OldMomentum.angle(-theGlobalNormal);
00770               // Round it to closest integer
00771               G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
00772 
00773               // Take random angles THETA and PHI, 
00774               // and see if below Probability - if not - Redo
00775               do {
00776                  thetaIndex = CLHEP::RandFlat::shootInt(thetaIndexMax-1);
00777                  phiIndex = CLHEP::RandFlat::shootInt(phiIndexMax-1);
00778                  // Find probability with the new indeces from LUT
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               // Rotate Photon Momentum in Theta, then in Phi
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               // Rotate Polarization too:
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;     // *** Snell's Law ***
00846            }
00847            else {
00848               sint1 = 0.0;
00849               sint2 = 0.0;
00850            }
00851 
00852            if (sint2 >= 1.0) {
00853 
00854               // Simulate total internal reflection
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               // Calculate amplitude for transmission (Q = P x N)
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                  // Here we Follow Jackson's conventions and we set the
00904                  // parallel component = 1 in case of a ray perpendicular
00905                  // to the surface
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                  // Simulate reflection
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) {   // incident ray oblique
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 {               // incident ray perpendicular
00963 
00964                        if (Rindex2 > Rindex1) {
00965                           NewPolarization = - OldPolarization;
00966                        }
00967                        else {
00968                           NewPolarization =   OldPolarization;
00969                        }
00970 
00971                     }
00972                  }
00973               }
00974               else { // photon gets transmitted
00975 
00976                  // Simulate transmission/refraction
00977 
00978                  Inside = !Inside;
00979                  Through = true;
00980                  theStatus = FresnelRefraction;
00981 
00982                  if (sint1 > 0.0) {      // incident ray oblique
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 {                  // incident ray perpendicular
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 // GetMeanFreePath
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);           //unit number 1
01082 
01083   G4complex numeratorTE;      // E1_perp=1 E1_parl=0 -> TE polarization
01084   G4complex numeratorTM;      // E1_parl=1 E1_perp=0 -> TM polarization
01085   G4complex denominatorTE, denominatorTM;
01086   G4complex rTM, rTE;
01087 
01088   // Following two equations, rTM and rTE, are from: "Introduction To Modern
01089   // Optics" written by Fowles
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   // This is my calculaton for reflectivity on a metalic surface
01102   // depending on the fraction of TE and TM polarization
01103   // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
01104   // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
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   // calculate FacetNormal
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      // Here we Follow Jackson's conventions and we set the
01161      // parallel component = 1 in case of a ray perpendicular
01162      // to the surface
01163      E1_perp  = 0.0;
01164      E1_parl  = 1.0;
01165   }
01166 
01167   //calculate incident angle
01168   G4double incidentangle = GetIncidentAngle();
01169 
01170   //calculate the reflectivity depending on incident angle,
01171   //polarization and complex refractive
01172 
01173   theReflectivity =
01174              GetReflectivity(E1_perp, E1_parl, incidentangle,
01175                                                  RealRindex, ImaginaryRindex);
01176 }

Generated on Mon May 27 17:49:09 2013 for Geant4 by  doxygen 1.4.7