Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-10 07:50:29 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 ////////////////////////////////////////////////////////////////////////
0027 // Optical Photon Boundary Process Class Implementation
0028 ////////////////////////////////////////////////////////////////////////
0029 //
0030 // File:        G4OpBoundaryProcess.cc
0031 // Description: Discrete Process -- reflection/refraction at
0032 //                                  optical interfaces
0033 // Version:     1.1
0034 // Created:     1997-06-18
0035 // Modified:    1998-05-25 - Correct parallel component of polarization
0036 //                           (thanks to: Stefano Magni + Giovanni Pieri)
0037 //              1998-05-28 - NULL Rindex pointer before reuse
0038 //                           (thanks to: Stefano Magni)
0039 //              1998-06-11 - delete *sint1 in oblique reflection
0040 //                           (thanks to: Giovanni Pieri)
0041 //              1998-06-19 - move from GetLocalExitNormal() to the new 
0042 //                           method: GetLocalExitNormal(&valid) to get
0043 //                           the surface normal in all cases
0044 //              1998-11-07 - NULL OpticalSurface pointer before use
0045 //                           comparison not sharp for: std::abs(cost1) < 1.0
0046 //                           remove sin1, sin2 in lines 556,567
0047 //                           (thanks to Stefano Magni)
0048 //              1999-10-10 - Accommodate changes done in DoAbsorption by
0049 //                           changing logic in DielectricMetal
0050 //              2001-10-18 - avoid Linux (gcc-2.95.2) warning about variables
0051 //                           might be used uninitialized in this function
0052 //                           moved E2_perp, E2_parl and E2_total out of 'if'
0053 //              2003-11-27 - Modified line 168-9 to reflect changes made to
0054 //                           G4OpticalSurface class ( by Fan Lei)
0055 //              2004-02-02 - Set theStatus = Undefined at start of DoIt
0056 //              2005-07-28 - add G4ProcessType to constructor
0057 //              2006-11-04 - add capability of calculating the reflectivity
0058 //                           off a metal surface by way of a complex index 
0059 //                           of refraction - Thanks to Sehwook Lee and John 
0060 //                           Hauptman (Dept. of Physics - Iowa State Univ.)
0061 //              2009-11-10 - add capability of simulating surface reflections
0062 //                           with Look-Up-Tables (LUT) containing measured
0063 //                           optical reflectance for a variety of surface
0064 //                           treatments - Thanks to Martin Janecek and
0065 //                           William Moses (Lawrence Berkeley National Lab.)
0066 //              2013-06-01 - add the capability of simulating the transmission
0067 //                           of a dichronic filter
0068 //              2017-02-24 - add capability of simulating surface reflections
0069 //                           with Look-Up-Tables (LUT) developed in DAVIS
0070 //
0071 // Author:      Peter Gumplinger
0072 //      adopted from work by Werner Keil - April 2/96
0073 // mail:        gum@triumf.ca
0074 //
0075 ////////////////////////////////////////////////////////////////////////
0076 
0077 #include "G4ios.hh"
0078 #include "G4PhysicalConstants.hh"
0079 #include "G4OpProcessSubType.hh"
0080 #include "InstrumentedG4OpBoundaryProcess.hh"
0081 #include "G4GeometryTolerance.hh"
0082 #include "G4VSensitiveDetector.hh"
0083 #include "G4ParallelWorldProcess.hh"
0084 #include "G4SystemOfUnits.hh"
0085 #include <csignal>
0086 
0087 
0088 #ifdef WITH_PMTFASTSIM
0089 #include "CustomART.h"
0090 #endif
0091 
0092 #include "SLOG.hh"
0093 #include "U4OpticalSurface.h"
0094 #include "U4OpBoundaryProcessStatus.h"
0095 #include "U4MaterialPropertiesTable.h"
0096 
0097 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0098 
0099 #include "scuda.h"
0100 #include "squad.h"
0101 #include "SEvt.hh"
0102 #include "SSys.hh"
0103 //#include "U4PhotonInfo.h"  // TODO: should be STrackInfo<spho> now
0104 const int InstrumentedG4OpBoundaryProcess::PIDX = SSys::getenvint("PIDX", -1) ; 
0105 
0106 #ifdef WITH_CUSTOM4
0107 #include "C4TrackInfo.h"
0108 #include "C4Pho.h"
0109 #else
0110 #include "STrackInfo.h"
0111 #include "spho.h"
0112 #endif
0113 #endif
0114 
0115 
0116 
0117 #include "U4UniformRand.h"
0118 NP* U4UniformRand::UU = nullptr ;  
0119 // UU gets set by U4Recorder::saveOrLoadStates when doing single photon reruns
0120 
0121 
0122 
0123 
0124 #ifdef DEBUG_TAG
0125 #include "U4Stack.h"
0126 #include "SEvt.hh"
0127 
0128 const plog::Severity InstrumentedG4OpBoundaryProcess::LEVEL = SLOG::EnvLevel("InstrumentedG4OpBoundaryProcess", "DEBUG" ); 
0129 
0130 const bool InstrumentedG4OpBoundaryProcess::FLOAT  = getenv("InstrumentedG4OpBoundaryProcess_FLOAT") != nullptr ;
0131 
0132 //void InstrumentedG4OpBoundaryProcess::ResetNumberOfInteractionLengthLeft(){ G4VProcess::ResetNumberOfInteractionLengthLeft(); }
0133 void InstrumentedG4OpBoundaryProcess::ResetNumberOfInteractionLengthLeft()
0134 {
0135     G4double u0 = G4UniformRand() ; 
0136 
0137     bool burn_enabled = true ; 
0138     if(burn_enabled)
0139     {
0140         int u0_idx = U4UniformRand::Find(u0, SEvt::UU ) ;  
0141         int u0_idx_delta = u0_idx - m_u0_idx ; 
0142  
0143         LOG(LEVEL) 
0144             << U4UniformRand::Desc(u0, SEvt::UU ) 
0145             << " u0_idx " << u0_idx 
0146             << " u0_idx_delta " << u0_idx_delta 
0147             << " BOP.RESET " 
0148             ; 
0149 
0150         int uu_burn = SEvt::UU_BURN ? SEvt::UU_BURN->ifind2D<int>(u0_idx, 0, 1 ) : -1  ; 
0151 
0152         if( uu_burn > 0 )
0153         {
0154             u0 = U4UniformRand::Burn(uu_burn);
0155             u0_idx = U4UniformRand::Find(u0, SEvt::UU ) ;  
0156             u0_idx_delta = u0_idx - m_u0_idx ;   
0157 
0158             LOG(LEVEL) 
0159                 << U4UniformRand::Desc(u0, SEvt::UU ) 
0160                 << " u0_idx " << u0_idx 
0161                 << " u0_idx_delta " << u0_idx_delta 
0162                 << " after uu_burn " << uu_burn
0163                 << " BOP.RESET " 
0164                 ; 
0165         } 
0166 
0167         m_u0 = u0 ; 
0168         m_u0_idx = u0_idx ; 
0169         m_u0_idx_delta = u0_idx_delta ; 
0170     }
0171 
0172 #ifndef PRODUCTION
0173     SEvt::AddTag( 1, U4Stack_BoundaryDiscreteReset, u0 );  
0174 #endif
0175 
0176     if(FLOAT)
0177     {   
0178         float f = -1.f*std::log( float(u0) ) ;   
0179         theNumberOfInteractionLengthLeft = f ; 
0180     }   
0181     else
0182     {   
0183         theNumberOfInteractionLengthLeft = -1.*G4Log(u0) ;   
0184     }   
0185     theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft; 
0186 
0187 }
0188 #endif
0189 
0190 #ifdef WITH_PMTFASTSIM
0191 double InstrumentedG4OpBoundaryProcess::getU0() const
0192 {
0193     return m_u0 ; 
0194 }
0195 int InstrumentedG4OpBoundaryProcess::getU0_idx() const
0196 {
0197     return m_u0_idx ; 
0198 }
0199 const double* InstrumentedG4OpBoundaryProcess::getRecoveredNormal() const 
0200 {
0201     return (const double*)&theRecoveredNormal ;
0202 }
0203 
0204 /**
0205 InstrumentedG4OpBoundaryProcess::getCustomStatus
0206 --------------------------------------------------
0207 
0208 HMM: actually makes more sense to hold a more 
0209 general status within InstrumentedG4OpBoundaryProcess 
0210 which describes whether the custom_boundary was used or not 
0211 and which is informative on why custom boundary used/not-used.
0212 That is more relevant than the R/T/A which is already described
0213 by the standard theStatus. 
0214 
0215 **/
0216 
0217 char InstrumentedG4OpBoundaryProcess::getCustomStatus() const 
0218 {
0219     return theCustomStatus ; 
0220     //return m_custom_boundary ? m_custom_boundary->customStatus : '-' ; 
0221 }
0222 void InstrumentedG4OpBoundaryProcess::Save(const char* fold) // static
0223 {
0224     //CustomBoundary<JPMT>::Save(fold); 
0225 }
0226 
0227 #endif
0228 
0229 
0230 InstrumentedG4OpBoundaryProcess::InstrumentedG4OpBoundaryProcess(
0231 #ifdef WITH_PMTFASTSIM
0232     const IPMTAccessor* accessor,
0233 #endif
0234     const G4String& processName, 
0235     G4ProcessType type
0236     ) 
0237     : 
0238     G4VDiscreteProcess(processName, type)
0239 #ifdef WITH_PMTFASTSIM
0240     ,SOpBoundaryProcess(processName.c_str())
0241 #endif
0242     ,theCustomStatus('U')
0243 #ifdef WITH_PMTFASTSIM
0244     ,m_custom_art(new CustomART(
0245                   accessor,
0246                   theTransmittance,
0247                   theReflectivity,
0248                   theEfficiency,
0249                   theGlobalPoint,
0250                   OldMomentum,
0251                   OldPolarization,
0252                   theRecoveredNormal,
0253                   thePhotonMomentum))
0254     ,m_u0(-1.)
0255     ,m_u0_idx(-1)
0256 #endif
0257     ,PostStepDoIt_count(-1)
0258 {
0259     LOG(LEVEL) << " processName " << GetProcessName()  ; 
0260 
0261     SetProcessSubType(fOpBoundary);
0262 
0263     theStatus = Undefined;
0264     theModel = glisur;
0265     theFinish = polished;
0266     theReflectivity =  1.;
0267     theEfficiency   =  0.;
0268     theTransmittance = 0.;
0269 
0270     theSurfaceRoughness = 0.;
0271 
0272     prob_sl = 0.;
0273     prob_ss = 0.;
0274     prob_bs = 0.;
0275 
0276     PropertyPointer  = NULL;
0277     PropertyPointer1 = NULL;
0278     PropertyPointer2 = NULL;
0279 
0280     Material1 = NULL;
0281     Material2 = NULL;
0282 
0283     theRecoveredNormal.set(0.,0.,0.); 
0284 
0285     OpticalSurface = NULL;
0286 
0287     kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 
0288 
0289     iTE = iTM = 0;
0290     thePhotonMomentum = 0.;
0291     Rindex1 = Rindex2 = 1.;
0292     cost1 = cost2 = sint1 = sint2 = 0.;
0293 
0294     idx = idy = 0;
0295     DichroicVector = NULL;
0296 
0297     fInvokeSD = true;
0298 }
0299 
0300 InstrumentedG4OpBoundaryProcess::~InstrumentedG4OpBoundaryProcess(){}
0301 
0302 
0303 /**
0304 InstrumentedG4OpBoundaryProcess::PostStepDoIt
0305 -----------------------------------------------
0306 
0307 Wrapper to help with the spagetti mush  
0308 
0309 **/
0310 
0311 G4VParticleChange* InstrumentedG4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
0312 {
0313 #ifdef WITH_PMTFASTSIM
0314     PostStepDoIt_count += 1 ;  
0315     spho* label = STrackInfo<spho>::GetRef(&aTrack) ; 
0316 
0317     LOG(LEVEL) 
0318         << "[ "
0319         << " PostStepDoIt_count " << PostStepDoIt_count
0320         << " label.desc " << label->desc()
0321         ; 
0322 #endif
0323     G4VParticleChange* change = PostStepDoIt_(aTrack, aStep) ; 
0324 
0325 #ifdef WITH_PMTFASTSIM
0326     LOG(LEVEL) << "] " << PostStepDoIt_count  ; 
0327 #endif
0328     return change ; 
0329 }
0330 
0331 /**
0332 InstrumentedG4OpBoundaryProcess::PostStepDoIt_
0333 -------------------------------------------------
0334 
0335 UNBELIEVABLY ABYSMAL CODING EVEN FOR GEANT4 
0336 
0337 
0338 head
0339     Material1, Material2, thePhotonMomentum, OldMomentum, OldPolarization
0340 theGlobalNormal
0341     note the flip 
0342 Rindex1  
0343     from Material1
0344 Defaults
0345     theReflectivity=1, theEfficiency=0, theTransmittance=0
0346 Surface
0347     border or skin  
0348 opsu
0349     detect if CustomART enabled based on OpticalSurface name
0350 OpticalSurface
0351     type, theModel, theFinish
0352 
0353     OpticalSurface.mpt
0354 
0355         OpticalSurface.mpt.backpainted
0356             Rindex2 from OpticalSurface.mpt   
0357             OpticalSurface with finish:polishedbackpainted/groundbackpainted
0358             require mpt with RINDEX to avoid NoRINDEX fStopAndKill
0359 
0360         OpticalSurface.mpt.theReflectivity,theTransmittance,theEfficiency
0361             sets theReflectivity
0362             sets theTransmittance
0363             sets theEfficiency
0364 
0365         OpticalSurface.mpt.CustomART
0366             custom:true call CustomART  
0367          
0368         OpticalSurface.mpt.unified
0369             theModel:unified sets prob_sl/prob_ss/prob_bs 
0370 
0371         OpticalSurface.no-mpt.backpainted
0372             fStopAndKill with ProposeLocalEnergyDeposit(thePhotonMomentum)  
0373 
0374 didi.polished/ground
0375     get Material2.mpt (not OpticalSurface.mpt)
0376     get Material2.mpt.Rindex2 otherwise NoRINDEX fStopAndKill
0377 
0378 type_switch
0379     type_switch.dime
0380         DielectricMetal
0381     type_switch.didi
0382         type_switch.didi.backpainted
0383              DielectricDielectric
0384         type_switch.didi.not_backpainted 
0385 
0386              rand-3-way depending on 
0387 
0388              theReflectivity 
0389              theReflectivity+theTransmittance         
0390 
0391                     R            R+T
0392              +------|-------------|----------------+
0393               R:refl    T:trans         A:abs               
0394 
0395              type_switch.didi.not_backpainted.A
0396                  DoAbsorption 
0397 
0398              type_switch.didi.not_backpainted.T
0399                  set theStatus=Transmission 
0400                  [UNNATURAL UNCHANGED DIR, POL]
0401 
0402                  NB: theTransmittance default is zero, so this needs a 
0403                  TRANSMITTANCE property to get it to happen
0404 
0405              type_switch.didi.not_backpainted.R
0406                  theFinish:(polished/ground)frontpainted
0407                      ground: set theStatus=LambertianReflection 
0408                      DoReflection
0409                  not-frontpainted
0410                      DielectricDielectric  
0411 
0412 tail
0413     aParticleChange NewMomentum, NewPolarization after .unit()
0414 
0415 
0416 **/
0417 
0418 G4VParticleChange* InstrumentedG4OpBoundaryProcess::PostStepDoIt_(const G4Track& aTrack, const G4Step& aStep)
0419 {
0420     //[head
0421 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0422     // formerly used U4PhotonInfo::GetIndex to pick up the label set 
0423     // in U4Recorder::PreUserTrackingAction_Optical for initially unlabelled input photons
0424     // pidx = U4PhotonInfo::GetIndex(&aTrack);    // TODO: now needs to be STrackInfo<spho>
0425 #ifdef WITH_CUSTOM4
0426     C4Pho* label = C4TrackInfo<C4Pho>::GetRef(&aTrack); 
0427     pidx = label->id ; 
0428 #else
0429     spho* label = STrackInfo<spho>::GetRef(&aTrack); 
0430     pidx = label->id ; 
0431 #endif
0432     pidx_dump = pidx == PIDX ; 
0433 #endif
0434     theStatus = Undefined;
0435     theCustomStatus = 'B' ; 
0436 
0437     aParticleChange.Initialize(aTrack);
0438     aParticleChange.ProposeVelocity(aTrack.GetVelocity());
0439 
0440     // Get hyperStep from  G4ParallelWorldProcess
0441     //  NOTE: PostSetpDoIt of this process should be
0442     //        invoked after G4ParallelWorldProcess!
0443 
0444     const G4Step* pStep = &aStep;
0445     const G4Step* hStep = G4ParallelWorldProcess::GetHyperStep();
0446         
0447     if (hStep) pStep = hStep;
0448     G4bool isOnBoundary = (pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
0449 
0450     if (isOnBoundary) 
0451     {
0452         Material1 = pStep->GetPreStepPoint()->GetMaterial();
0453         Material2 = pStep->GetPostStepPoint()->GetMaterial();
0454     } 
0455     else 
0456     {
0457         theStatus = NotAtBoundary;
0458         if ( verboseLevel > 0) BoundaryProcessVerbose();
0459         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0460     }
0461 
0462     G4VPhysicalVolume* thePrePV  = pStep->GetPreStepPoint() ->GetPhysicalVolume(); 
0463     G4VPhysicalVolume* thePostPV = pStep->GetPostStepPoint()->GetPhysicalVolume(); 
0464     G4bool haveEnteredDaughter= (thePostPV->GetMotherLogical() == thePrePV ->GetLogicalVolume()); // SCB
0465 
0466     LOG(LEVEL)
0467         << " PostStepDoIt_count " << PostStepDoIt_count
0468         << " thePrePV " << ( thePrePV ? thePrePV->GetName() : "-" )
0469         << " thePostPV " << ( thePostPV ? thePostPV->GetName() : "-" )
0470         << " haveEnteredDaughter " << haveEnteredDaughter
0471         << " kCarTolerance/2 " << kCarTolerance/2
0472         ;
0473 
0474     if (aTrack.GetStepLength()<=kCarTolerance/2)
0475     {
0476         theStatus = StepTooSmall;
0477         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0478     }
0479 
0480     const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0481     thePhotonMomentum = aParticle->GetTotalMomentum();
0482     OldMomentum       = aParticle->GetMomentumDirection();
0483     OldPolarization   = aParticle->GetPolarization();
0484 
0485     //]head
0486 
0487     //[theGlobalNorml
0488     theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
0489 
0490     G4bool valid;
0491     // Use the new method for Exit Normal in global coordinates,
0492     // which provides the normal more reliably.
0493     // ID of Navigator which limits step
0494 
0495     G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
0496     std::vector<G4Navigator*>::iterator iNav =
0497                 G4TransportationManager::GetTransportationManager()->
0498                                          GetActiveNavigatorsIterator();
0499     theGlobalExitNormal =
0500                    (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
0501 
0502     // theGlobalExitNormal is already oriented by G4Navigator to point from vol1 -> vol2 
0503     // so try to undo that flip by G4Navigator in order to recover the original geometry 
0504     // normal that is independent of G4Track direction
0505 
0506     theRecoveredNormal = ( haveEnteredDaughter ? -1. : 1. )* theGlobalExitNormal  ; 
0507 
0508 
0509 
0510     theGlobalNormal = theGlobalExitNormal ;  
0511 
0512 
0513 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0514     {
0515         quad2& prd = SEvt::Get_ECPU()->current_prd ; 
0516         prd.q0.f.x = theGlobalNormal.x() ; 
0517         prd.q0.f.y = theGlobalNormal.y() ; 
0518         prd.q0.f.z = theGlobalNormal.z() ; 
0519 
0520         // TRY USING PRE->POST POSITION CHANGE TO GET THE PRD DISTANCE ? 
0521         G4ThreeVector theGlobalPoint_pre = pStep->GetPreStepPoint()->GetPosition();
0522         G4ThreeVector theGlobalPoint_delta = theGlobalPoint - theGlobalPoint_pre  ;  
0523         prd.q0.f.w = theGlobalPoint_delta.mag() ; 
0524 
0525         // HMM: PRD intersect identity ? how to mimic what Opticks does ? 
0526     }
0527 #endif
0528 
0529     if (valid) 
0530     {
0531         theGlobalNormal = -theGlobalNormal;
0532     }
0533     else 
0534     {
0535         G4ExceptionDescription ed;
0536         ed << " InstrumentedG4OpBoundaryProcess/PostStepDoIt(): "
0537            << " The Navigator reports that it returned an invalid normal"
0538            << G4endl;
0539         G4Exception("InstrumentedG4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
0540                       EventMustBeAborted,ed,
0541                       "Invalid Surface Normal - Geometry must return valid surface normal");
0542     }
0543 
0544 
0545 
0546     if (OldMomentum * theGlobalNormal > 0.0) 
0547     {
0548 #ifdef G4OPTICAL_DEBUG
0549         G4ExceptionDescription ed;
0550         ed << " InstrumentedG4OpBoundaryProcess/PostStepDoIt(): "
0551            << " theGlobalNormal points in a wrong direction. "
0552            << G4endl;
0553         ed << "    The momentum of the photon arriving at interface (oldMomentum)"
0554            << " must exit the volume cross in the step. " << G4endl;
0555         ed << "  So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
0556         ed << "  >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
0557         ed << "     Old Momentum  (during step)     = " << OldMomentum << G4endl;
0558         ed << "     Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
0559         ed << G4endl;
0560         G4Exception("InstrumentedG4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
0561                        EventMustBeAborted,  // Or JustWarning to see if it happens repeatedbly on one ray
0562                        ed,
0563                       "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
0564 #else
0565         theGlobalNormal = -theGlobalNormal;
0566 #endif
0567     }
0568     //]theGlobalNorml
0569 
0570 
0571     //[Rindex1
0572     G4MaterialPropertiesTable* aMaterialPropertiesTable;
0573     G4MaterialPropertyVector* Rindex;
0574 
0575     aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
0576     if (aMaterialPropertiesTable) 
0577     {
0578         Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX);
0579     }
0580     else 
0581     {
0582         theStatus = NoRINDEX;
0583         if ( verboseLevel > 0) BoundaryProcessVerbose();
0584         aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
0585         aParticleChange.ProposeTrackStatus(fStopAndKill);
0586         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0587     }
0588 
0589     if (Rindex) 
0590     {
0591         Rindex1 = Rindex->Value(thePhotonMomentum);
0592     }
0593     else 
0594     {
0595         theStatus = NoRINDEX;
0596         if ( verboseLevel > 0) BoundaryProcessVerbose();
0597         aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
0598         aParticleChange.ProposeTrackStatus(fStopAndKill);
0599         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0600     }
0601     //]Rindex1
0602     // SCB: UNBELIEVABLY POOR CODING STYLE : POINTLESS REPETITION OF BLOCKS OF CODE
0603 
0604     //[Defaults
0605     theReflectivity =  1.;
0606     theEfficiency   =  0.;
0607     theTransmittance = 0.;
0608 
0609     theSurfaceRoughness = 0.;
0610 
0611     theModel = glisur;
0612     theFinish = polished;
0613 
0614     G4SurfaceType type = dielectric_dielectric;
0615 
0616     //]Defaults
0617 
0618     //[Surface
0619     Rindex = NULL;
0620     OpticalSurface = NULL;
0621 
0622     G4LogicalSurface* Surface = NULL;
0623     Surface = G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
0624     if (Surface == NULL)
0625     {
0626         G4bool enteredDaughter= (thePostPV->GetMotherLogical() == thePrePV ->GetLogicalVolume());
0627         if(enteredDaughter)
0628         {
0629             Surface = G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
0630             if(Surface == NULL) Surface = G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
0631         }
0632         else 
0633         {
0634             Surface = G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
0635             if(Surface == NULL) Surface = G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
0636         }
0637     }
0638     //]Surface
0639 
0640 
0641     if (Surface) OpticalSurface = dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty()); 
0642     LOG(LEVEL) 
0643         << " PostStepDoIt_count " << PostStepDoIt_count
0644         << " Surface " << ( Surface ? Surface->GetName() : "-" ) 
0645         << " OpticalSurface " << ( OpticalSurface ? OpticalSurface->GetName() : "-" ) 
0646         ; 
0647 
0648     //[OpticalSurface
0649     if (OpticalSurface) 
0650     {
0651         const char* OpticalSurfaceName = OpticalSurface->GetName().c_str() ;  // GetName by ref, so not transient 
0652         LOG(LEVEL) 
0653             << " PostStepDoIt_count " << PostStepDoIt_count << " "
0654             << " OpticalSurfaceName " << OpticalSurfaceName
0655             << U4OpticalSurface::Desc(OpticalSurface) 
0656             ; 
0657 
0658         type      = OpticalSurface->GetType();
0659         theModel  = OpticalSurface->GetModel();
0660         theFinish = OpticalSurface->GetFinish();
0661 
0662         aMaterialPropertiesTable = OpticalSurface->GetMaterialPropertiesTable(); 
0663 
0664         //[OpticalSurface.mpt
0665         if (aMaterialPropertiesTable) 
0666         {
0667             /*
0668             LOG(LEVEL) 
0669                 << " PostStepDoIt_count " << PostStepDoIt_count << " " 
0670                 << U4MaterialPropertiesTable::Desc(aMaterialPropertiesTable) 
0671                 ; 
0672             */
0673 
0674             //[OpticalSurface.mpt.backpainted
0675             if (theFinish == polishedbackpainted || theFinish == groundbackpainted ) 
0676             {
0677                 Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX);
0678                 if (Rindex) 
0679                 {
0680                     Rindex2 = Rindex->Value(thePhotonMomentum);
0681                 }
0682                 else 
0683                 {
0684                     theStatus = NoRINDEX;
0685                     if ( verboseLevel > 0) BoundaryProcessVerbose();
0686                     aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
0687                     aParticleChange.ProposeTrackStatus(fStopAndKill);
0688                     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0689                 }
0690             }
0691             //]OpticalSurface.mpt.backpainted
0692 
0693             //[OpticalSurface.mpt.theReflectivity,theTransmittance,theEfficiency
0694             PropertyPointer  = aMaterialPropertiesTable->GetProperty(kREFLECTIVITY); 
0695             PropertyPointer1 = aMaterialPropertiesTable->GetProperty(kREALRINDEX); 
0696             PropertyPointer2 = aMaterialPropertiesTable->GetProperty(kIMAGINARYRINDEX); 
0697 
0698             /*
0699             LOG(LEVEL) 
0700                 << " PostStepDoIt_count " << PostStepDoIt_count
0701                 << " PropertyPointer.kREFLECTIVITY " << PropertyPointer
0702                 << " PropertyPointer1.kREALRINDEX " << PropertyPointer1 
0703                 << " PropertyPointer2.kIMAGINARYRINDEX " << PropertyPointer2
0704                 ;
0705             */
0706 
0707             iTE = 1;
0708             iTM = 1;
0709 
0710             if (PropertyPointer) 
0711             {
0712                 theReflectivity = PropertyPointer->Value(thePhotonMomentum); 
0713             }
0714             else if (PropertyPointer1 && PropertyPointer2) 
0715             {
0716                 CalculateReflectivity();  // sets theReflectivity 
0717             }
0718 
0719 
0720             PropertyPointer = aMaterialPropertiesTable->GetProperty(kEFFICIENCY); 
0721             if (PropertyPointer) 
0722             {
0723                 theEfficiency = PropertyPointer->Value(thePhotonMomentum); 
0724             }
0725 
0726             PropertyPointer = aMaterialPropertiesTable->GetProperty(kTRANSMITTANCE); 
0727             if (PropertyPointer) 
0728             {
0729                 theTransmittance = PropertyPointer->Value(thePhotonMomentum); 
0730             }
0731             //]OpticalSurface.mpt.theReflectivity,theTransmittance,theEfficiency
0732 
0733             //[OpticalSurface.mpt.CustomBoundary
0734 #ifdef WITH_PMTFASTSIM
0735             char osn = OpticalSurfaceName[0] ; 
0736             if(  osn == '@' || osn == '#' )  // only customize specially named OpticalSurfaces 
0737             {
0738                 if( m_custom_art->local_z(aTrack) < 0. ) // lower hemi : No customization, standard boundary  
0739                 {
0740                     theCustomStatus = 'Z' ; 
0741                 }
0742                 else if( osn == '@') //  upper hemi with name starting @ : MultiFilm ART transmit thru into PMT
0743                 {
0744                     theCustomStatus = 'Y' ;
0745      
0746                     m_custom_art->doIt(aTrack, aStep) ;  // calculate theReflectivity theTransmittance theEfficiency 
0747 
0748                     type = dielectric_dielectric ;  
0749                     theModel = glisur ; 
0750                     theFinish = polished ;  
0751                     // guide thru the below jungle : only when custom handling is triggered 
0752                 }
0753                 else if( osn == '#' ) // upper hemi with name starting # : Traditional Detection at photocathode
0754                 {
0755                     theCustomStatus = 'Y' ;
0756 
0757                     type = dielectric_metal ; 
0758                     theModel = glisur ; 
0759                     theReflectivity = 0. ; 
0760                     theTransmittance = 0. ; 
0761                     theEfficiency = 1. ; 
0762                 }
0763             }
0764             else
0765             {
0766                 theCustomStatus = 'X' ; 
0767             }
0768 #else
0769             theCustomStatus = 'X' ; 
0770 #endif
0771             //]OpticalSurface.mpt.CustomBoundary
0772 
0773             LOG(LEVEL)
0774                 << " PostStepDoIt_count " << PostStepDoIt_count
0775                 << " theReflectivity " << theReflectivity
0776                 << " theEfficiency " << theEfficiency
0777                 << " theTransmittance " << theTransmittance
0778                 << " theCustomStatus " << theCustomStatus
0779                 ; 
0780 
0781             if (aMaterialPropertiesTable->ConstPropertyExists("SURFACEROUGHNESS"))
0782                 theSurfaceRoughness = aMaterialPropertiesTable-> GetConstProperty(kSURFACEROUGHNESS); 
0783 
0784             //[OpticalSurface.mpt.unified
0785             if ( theModel == unified ) 
0786             {
0787                 PropertyPointer = aMaterialPropertiesTable->GetProperty(kSPECULARLOBECONSTANT); 
0788                 prob_sl = PropertyPointer ? PropertyPointer->Value(thePhotonMomentum) : 0.0 ; 
0789 
0790                 PropertyPointer = aMaterialPropertiesTable->GetProperty(kSPECULARSPIKECONSTANT); 
0791                 prob_ss = PropertyPointer ? PropertyPointer->Value(thePhotonMomentum) : 0.0 ; 
0792 
0793                 PropertyPointer = aMaterialPropertiesTable->GetProperty(kBACKSCATTERCONSTANT); 
0794                 prob_bs = PropertyPointer ? PropertyPointer->Value(thePhotonMomentum) : 0.0 ; 
0795             }
0796             //]OpticalSurface.mpt.unified
0797 
0798 
0799 
0800 
0801 
0802 
0803 
0804         }
0805         //]OpticalSurface.mpt
0806         //[OpticalSurface.no-mpt.backpainted
0807         else if (theFinish == polishedbackpainted || theFinish == groundbackpainted ) 
0808         {
0809             aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
0810             aParticleChange.ProposeTrackStatus(fStopAndKill);
0811             return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0812         }
0813         //]OpticalSurface.no-mpt.backpainted
0814     }
0815     //]OpticalSurface
0816 
0817     LOG(LEVEL) 
0818         << " PostStepDoIt_count " << PostStepDoIt_count
0819         << " after OpticalSurface if "  
0820         ;
0821 
0822 
0823     //[didi.polished/ground
0824     if (type == dielectric_dielectric ) 
0825     {
0826         if (theFinish == polished || theFinish == ground ) 
0827         {
0828             if (Material1 == Material2)
0829             {
0830                 theStatus = SameMaterial;
0831                 if ( verboseLevel > 0) BoundaryProcessVerbose();
0832                 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0833             }
0834             aMaterialPropertiesTable = Material2->GetMaterialPropertiesTable(); 
0835 
0836             if (aMaterialPropertiesTable) Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX); 
0837             if (Rindex)
0838             {
0839                 Rindex2 = Rindex->Value(thePhotonMomentum);
0840             }
0841             else 
0842             {
0843                 theStatus = NoRINDEX;
0844                 if ( verboseLevel > 0) BoundaryProcessVerbose();
0845                 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
0846                 aParticleChange.ProposeTrackStatus(fStopAndKill);
0847                 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0848             }
0849         }
0850         LOG(LEVEL) 
0851             << " PostStepDoIt_count " << PostStepDoIt_count
0852             << " didi.polished/ground " 
0853             << " Rindex2 " << Rindex2 
0854             ;  
0855     }
0856     //]didi.polished/ground
0857 
0858 
0859     //[type_switch 
0860 #ifdef WITH_PMTFASTSIM
0861     if( theCustomStatus == 'Y' )
0862     {
0863         G4double rand = G4UniformRand();
0864 
0865         G4double A = 1. - (theReflectivity + theTransmittance) ;  
0866 
0867         if ( rand < A )  // HMM: more normally rand > theReflectivity + theTransmittance 
0868         {
0869             DoAbsorption();   // theStatus is set to Detection/Absorption depending on a random and theEfficiency  
0870         }
0871         else
0872         {
0873             DielectricDielectric();
0874         }
0875     }
0876     else 
0877 #endif
0878     if (type == dielectric_metal) 
0879     {
0880         //[type_switch.dime
0881         DielectricMetal();
0882         //]type_switch.dime
0883     }
0884     else if (type == dielectric_LUT) 
0885     {
0886         DielectricLUT();
0887     }
0888     else if (type == dielectric_LUTDAVIS) 
0889     {
0890         DielectricLUTDAVIS();
0891     }
0892     else if (type == dielectric_dichroic) 
0893     {
0894         DielectricDichroic();
0895     }
0896     else if (type == dielectric_dielectric) 
0897     {
0898         //[type_switch.didi
0899         if ( theFinish == polishedbackpainted || theFinish == groundbackpainted ) 
0900         {
0901             //[type_switch.didi.backpainted
0902             DielectricDielectric();
0903             //]type_switch.didi.backpainted
0904         }
0905         else
0906         {   
0907             //[type_switch.didi.not_backpainted
0908             G4double rand = G4UniformRand();
0909             LOG(LEVEL) 
0910                 << " PostStepDoIt_count " << PostStepDoIt_count
0911                 << " didi.rand " << U4UniformRand::Desc(rand) 
0912                 << " theReflectivity " << std::setw(10) << std::fixed << std::setprecision(4) << theReflectivity
0913                 << " theTransmittance " << std::setw(10) << std::fixed << std::setprecision(4) << theTransmittance
0914                 ; 
0915 #ifndef PRODUCTION
0916 #ifdef DEBUG_TAG
0917             SEvt::AddTag(1, U4Stack_BoundaryBurn_SurfaceReflectTransmitAbsorb, rand );  
0918 #endif
0919 #endif
0920 
0921 #ifndef PRODUCTION
0922 #ifdef DEBUG_PIDX
0923             // SCB theReflectivity default is 1. so  "rand > theReflectivity" always false
0924             //     meaning that "rand" is always a burn doing nothing.  
0925             //     Unless a surface is associated which changes theReflectivity to be less than 1. 
0926             //     which can make the "rand" actually control the Reflect-OR-Absorb-OR-Transmit decision.  
0927             //
0928             LOG_IF(LEVEL, pidx_dump) 
0929                 << " DiDi0 " 
0930                 << " pidx " << std::setw(6) << pidx 
0931                 << " rand " << std::setw(10) << std::fixed << std::setprecision(5) << rand  
0932                 << " theReflectivity " << std::setw(10) << std::fixed << std::setprecision(4) << theReflectivity
0933                 << " rand > theReflectivity  " << (rand > theReflectivity )
0934                 ;
0935 #endif
0936 #endif
0937             if ( rand > theReflectivity )
0938             {
0939                 if (rand > theReflectivity + theTransmittance) 
0940                 {
0941                     //[type_switch.didi.not_backpainted.A
0942                     DoAbsorption();
0943                     //]type_switch.didi.not_backpainted.A
0944                 }
0945                 else 
0946                 {
0947                     //[type_switch.didi.not_backpainted.T
0948                     theStatus = Transmission;
0949                     NewMomentum = OldMomentum;
0950                     NewPolarization = OldPolarization;
0951                     LOG(LEVEL) << " mystifying Transmission with Mom and Pol unchanged ? " ; 
0952                     //]type_switch.didi.not_backpainted.T
0953                 }
0954             }
0955             else 
0956             {
0957                 //[type_switch.didi.not_backpainted.R
0958                 if ( theFinish == polishedfrontpainted ) 
0959                 {
0960                     DoReflection();
0961                 }
0962                 else if ( theFinish == groundfrontpainted ) 
0963                 {
0964                     theStatus = LambertianReflection;
0965                     DoReflection();
0966                 }
0967                 else 
0968                 {
0969                     DielectricDielectric();
0970                 }
0971                 //]type_switch.didi.not_backpainted.R
0972             }
0973             //]type_switch.didi.not_backpainted
0974         }   
0975         //[type_switch.didi
0976     }
0977     else 
0978     {
0979         //[type_switch.illegal
0980         G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
0981         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0982         //]type_switch.illegal
0983     }
0984     //]type switch 
0985 
0986     //[tail 
0987     NewMomentum = NewMomentum.unit();
0988     NewPolarization = NewPolarization.unit();
0989 
0990     aParticleChange.ProposeMomentumDirection(NewMomentum);
0991     aParticleChange.ProposePolarization(NewPolarization);
0992 
0993     if ( theStatus == FresnelRefraction || theStatus == Transmission ) 
0994     {
0995         G4MaterialPropertyVector* groupvel = Material2->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
0996         G4double finalVelocity = groupvel->Value(thePhotonMomentum);
0997         aParticleChange.ProposeVelocity(finalVelocity);
0998     }
0999     if ( theStatus == Detection && fInvokeSD ) InvokeSD(pStep);
1000 
1001     // ]tail 
1002     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
1003 }
1004 
1005 
1006 /**
1007 InstrumentedG4OpBoundaryProcess::GetFacetNormal
1008 ------------------------------------------------
1009 
1010 
1011 This function code alpha to a random value taken from the
1012 distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
1013 for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
1014 is a gaussian distribution with mean 0 and standard deviation
1015 sigma_alpha.  
1016 
1017 
1018 **/
1019 
1020 G4ThreeVector InstrumentedG4OpBoundaryProcess::GetFacetNormal(
1021      const G4ThreeVector& Momentum, const G4ThreeVector&  Normal ) const 
1022 {
1023     G4ThreeVector FacetNormal;
1024 
1025     if (theModel == unified || theModel == LUT || theModel== DAVIS) 
1026     {
1027         G4double alpha;
1028         G4double sigma_alpha = 0.0;
1029         if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
1030 
1031         if (sigma_alpha == 0.0) return FacetNormal = Normal;
1032 
1033         G4double f_max = std::min(1.0,4.*sigma_alpha);
1034 
1035         G4double phi, SinAlpha, CosAlpha, SinPhi, CosPhi, unit_x, unit_y, unit_z;
1036         G4ThreeVector tmpNormal;
1037 
1038         do 
1039         {
1040             do 
1041             {
1042                  alpha = G4RandGauss::shoot(0.0,sigma_alpha);
1043                  // Loop checking, 13-Aug-2015, Peter Gumplinger
1044             } 
1045             while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
1046 
1047             phi = G4UniformRand()*twopi;
1048 
1049             SinAlpha = std::sin(alpha);
1050             CosAlpha = std::cos(alpha);
1051             SinPhi = std::sin(phi);
1052             CosPhi = std::cos(phi);
1053 
1054             unit_x = SinAlpha * CosPhi;
1055             unit_y = SinAlpha * SinPhi;
1056             unit_z = CosAlpha;
1057 
1058             FacetNormal.setX(unit_x);
1059             FacetNormal.setY(unit_y);
1060             FacetNormal.setZ(unit_z);
1061 
1062             tmpNormal = Normal;
1063 
1064             FacetNormal.rotateUz(tmpNormal);
1065               // Loop checking, 13-Aug-2015, Peter Gumplinger
1066         } 
1067         while (Momentum * FacetNormal >= 0.0);
1068     }
1069     else 
1070     {
1071         G4double  polish = OpticalSurface ? OpticalSurface->GetPolish() : 1.0 ;
1072         if (polish < 1.0) 
1073         {
1074             do 
1075             {
1076                 G4ThreeVector smear;
1077                 do 
1078                 {
1079                     smear.setX(2.*G4UniformRand()-1.0);
1080                     smear.setY(2.*G4UniformRand()-1.0);
1081                     smear.setZ(2.*G4UniformRand()-1.0);
1082                     // Loop checking, 13-Aug-2015, Peter Gumplinger
1083                 } 
1084                 while (smear.mag()>1.0);
1085                 smear = (1.-polish) * smear;
1086                 FacetNormal = Normal + smear;
1087                 // Loop checking, 13-Aug-2015, Peter Gumplinger
1088             } 
1089             while (Momentum * FacetNormal >= 0.0);
1090             FacetNormal = FacetNormal.unit();
1091         }
1092         else 
1093         {
1094             FacetNormal = Normal;
1095         }
1096     }
1097     return FacetNormal;
1098 }
1099 
1100 
1101 /**
1102 InstrumentedG4OpBoundaryProcess::DielectricMetal
1103 --------------------------------------------------
1104 
1105 Changes::
1106 
1107    theStatus
1108    NewMomentum
1109    NewPolarization
1110 
1111 **/
1112 
1113 void InstrumentedG4OpBoundaryProcess::DielectricMetal()
1114 {
1115     LOG(LEVEL) 
1116        << " PostStepDoIt_count " << PostStepDoIt_count
1117        ;
1118 
1119     G4int n = 0;
1120     G4double rand, PdotN, EdotN;
1121     G4ThreeVector A_trans, A_paral;
1122 
1123     do 
1124     {
1125         n++;
1126 
1127         rand = G4UniformRand();
1128 
1129         LOG(LEVEL) 
1130             << " PostStepDoIt_count " << PostStepDoIt_count
1131             << " do-while n " << n 
1132             << " rand " << U4UniformRand::Desc(rand) 
1133             << " theReflectivity " << theReflectivity
1134             << " theTransmittance " << theTransmittance
1135             ;
1136 
1137 #ifndef PRODUCTION
1138 #ifdef DEBUG_TAG
1139         SEvt::AddTag(1, U4Stack_BoundaryDiMeReflectivity, rand); 
1140 #endif
1141 #endif
1142 
1143 
1144         if ( rand > theReflectivity && n == 1 ) 
1145         {
1146             if (rand > theReflectivity + theTransmittance) 
1147             {
1148                 DoAbsorption();
1149             } 
1150             else 
1151             {
1152                 theStatus = Transmission;
1153                 NewMomentum = OldMomentum;
1154                 NewPolarization = OldPolarization;
1155             }
1156             LOG(LEVEL) << " rand > theReflectivity && n == 1  break " ; 
1157             break;
1158         }
1159         else 
1160         {
1161             if (PropertyPointer1 && PropertyPointer2) 
1162             {
1163                 if ( n > 1 ) 
1164                 {
1165                     CalculateReflectivity();
1166                     if ( !G4BooleanRand_theReflectivity(theReflectivity) ) 
1167                     {
1168                         DoAbsorption();
1169                         break;
1170                     }
1171                 }
1172             }
1173 
1174             if ( theModel == glisur || theFinish == polished ) 
1175             {
1176                 DoReflection();
1177             } 
1178             else 
1179             {
1180                 if ( n == 1 ) ChooseReflection();
1181                                                                                 
1182                 if ( theStatus == LambertianReflection ) 
1183                 {
1184                     DoReflection();
1185                 }
1186                 else if ( theStatus == BackScattering ) 
1187                 {
1188                     NewMomentum = -OldMomentum;
1189                     NewPolarization = -OldPolarization;
1190                 }
1191                 else 
1192                 {
1193                     if(theStatus==LobeReflection)
1194                     {
1195                         if ( PropertyPointer1 && PropertyPointer2 )
1196                         {
1197                         } 
1198                         else 
1199                         {
1200                             theFacetNormal = GetFacetNormal(OldMomentum,theGlobalNormal);
1201                         }
1202                     }
1203 
1204                     PdotN = OldMomentum * theFacetNormal;
1205                     NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1206                     EdotN = OldPolarization * theFacetNormal;
1207 
1208                     if (sint1 > 0.0 ) 
1209                     {
1210                         A_trans = OldMomentum.cross(theFacetNormal);
1211                         A_trans = A_trans.unit();
1212                     } 
1213                     else 
1214                     {
1215                         A_trans  = OldPolarization;
1216                     }
1217                     A_paral   = NewMomentum.cross(A_trans);
1218                     A_paral   = A_paral.unit();
1219 
1220                     if(iTE>0&&iTM>0) 
1221                     {
1222                         NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1223                     } 
1224                     else if (iTE>0) 
1225                     {
1226                         NewPolarization = -A_trans;
1227                     } 
1228                     else if (iTM>0) 
1229                     {
1230                         NewPolarization = -A_paral;
1231                     }
1232                 }
1233             }
1234 
1235             //LOG(LEVEL) << " dime: Old<-New ? UGLY " ; 
1236             OldMomentum = NewMomentum;
1237             OldPolarization = NewPolarization;
1238         }
1239         //LOG(LEVEL) << " while check " ; 
1240           // Loop checking, 13-Aug-2015, Peter Gumplinger
1241     } while (NewMomentum * theGlobalNormal < 0.0);
1242     //LOG(LEVEL) << " dime: after while  " ; 
1243 }
1244 
1245 void InstrumentedG4OpBoundaryProcess::DielectricLUT()
1246 {
1247     G4int thetaIndex, phiIndex;
1248     G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
1249     G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
1250 
1251     theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) + (G4int(NoRINDEX)-G4int(groundbackpainted))); 
1252 
1253     G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
1254     G4int phiIndexMax   = OpticalSurface->GetPhiIndexMax();
1255 
1256     G4double rand;
1257 
1258     do 
1259     {
1260         rand = G4UniformRand();
1261         if ( rand > theReflectivity ) 
1262         {
1263             if (rand > theReflectivity + theTransmittance) 
1264             {
1265                 DoAbsorption();
1266             } 
1267             else 
1268             {
1269                 theStatus = Transmission;
1270                 NewMomentum = OldMomentum;
1271                 NewPolarization = OldPolarization;
1272             }
1273             break;
1274         }
1275         else 
1276         {
1277             // Calculate Angle between Normal and Photon Momentum
1278             G4double anglePhotonToNormal = 
1279                                           OldMomentum.angle(-theGlobalNormal);
1280             // Round it to closest integer
1281             G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
1282 
1283             // Take random angles THETA and PHI, 
1284             // and see if below Probability - if not - Redo
1285             do 
1286             {
1287                 thetaIndex = G4RandFlat::shootInt(thetaIndexMax-1);
1288                 phiIndex = G4RandFlat::shootInt(phiIndexMax-1);
1289                 // Find probability with the new indeces from LUT
1290                 AngularDistributionValue = OpticalSurface -> 
1291                    GetAngularDistributionValue(angleIncident,
1292                                                thetaIndex,
1293                                                phiIndex);
1294                 // Loop checking, 13-Aug-2015, Peter Gumplinger
1295             } 
1296             while ( !G4BooleanRand(AngularDistributionValue) );
1297 
1298             thetaRad = (-90 + 4*thetaIndex)*pi/180;
1299             phiRad = (-90 + 5*phiIndex)*pi/180;
1300             // Rotate Photon Momentum in Theta, then in Phi
1301             NewMomentum = -OldMomentum;
1302 
1303             PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
1304             if (PerpendicularVectorTheta.mag() < kCarTolerance )
1305                           PerpendicularVectorTheta = NewMomentum.orthogonal();
1306             NewMomentum =
1307                  NewMomentum.rotate(anglePhotonToNormal-thetaRad,
1308                                     PerpendicularVectorTheta);
1309             PerpendicularVectorPhi = 
1310                                   PerpendicularVectorTheta.cross(NewMomentum);
1311             NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
1312 
1313             // Rotate Polarization too:
1314             theFacetNormal = (NewMomentum - OldMomentum).unit();
1315             EdotN = OldPolarization * theFacetNormal;
1316             NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1317         }
1318         // Loop checking, 13-Aug-2015, Peter Gumplinger
1319     } 
1320     while (NewMomentum * theGlobalNormal <= 0.0);
1321 }
1322 
1323 void InstrumentedG4OpBoundaryProcess::DielectricLUTDAVIS()
1324 {
1325     G4int angindex, random, angleIncident;
1326     G4double ReflectivityValue, elevation, azimuth, EdotN;
1327     G4double anglePhotonToNormal;
1328 
1329     G4int LUTbin = OpticalSurface->GetLUTbins();
1330     G4double rand = G4UniformRand();
1331     do 
1332     {
1333         anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
1334         angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
1335 
1336         ReflectivityValue = OpticalSurface -> GetReflectivityLUTValue(angleIncident);
1337 
1338         if ( rand > ReflectivityValue ) 
1339         {
1340             if ( theEfficiency > 0 ) 
1341             {
1342                 DoAbsorption();
1343                 break;
1344             }
1345             else 
1346             {
1347                 theStatus = Transmission;
1348 
1349                 if (angleIncident <= 0.01) 
1350                 {
1351                     NewMomentum = OldMomentum;
1352                     break;
1353                 }
1354 
1355                 do 
1356                 {
1357                     random   = G4RandFlat::shootInt(1,LUTbin+1);
1358                     angindex = (((random*2)-1))+angleIncident*LUTbin*2 + 3640000;
1359 
1360                     azimuth  = OpticalSurface -> GetAngularDistributionValueLUT(angindex-1);
1361                     elevation= OpticalSurface -> GetAngularDistributionValueLUT(angindex);
1362 
1363                 } while ( elevation == 0 && azimuth == 0);
1364 
1365                 NewMomentum = -OldMomentum;
1366 
1367                 G4ThreeVector v = theGlobalNormal.cross(-NewMomentum);
1368                 G4ThreeVector vNorm = v/v.mag();
1369                 G4ThreeVector u = vNorm.cross(theGlobalNormal);
1370 
1371                 u = u *= (std::sin(elevation) * std::cos(azimuth));
1372                 v = vNorm *= (std::sin(elevation) * std::sin(azimuth));
1373                 G4ThreeVector w = theGlobalNormal *= (std::cos(elevation));
1374                 NewMomentum = G4ThreeVector(u+v+w);
1375 
1376                 // Rotate Polarization too:
1377                 theFacetNormal = (NewMomentum - OldMomentum).unit();
1378                 EdotN = OldPolarization * theFacetNormal;
1379                 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1380             }
1381         }
1382         else 
1383         {
1384             theStatus = LobeReflection;
1385             if (angleIncident == 0) 
1386             {
1387                 NewMomentum = -OldMomentum;
1388                 break;
1389             }
1390 
1391             do 
1392             {
1393                 random   = G4RandFlat::shootInt(1,LUTbin+1);
1394                 angindex = (((random*2)-1))+(angleIncident-1)*LUTbin*2;
1395 
1396                 azimuth   = OpticalSurface -> GetAngularDistributionValueLUT(angindex-1);
1397                 elevation = OpticalSurface -> GetAngularDistributionValueLUT(angindex);
1398             } 
1399             while (elevation == 0 && azimuth == 0);
1400 
1401             NewMomentum = -OldMomentum;
1402 
1403             G4ThreeVector v     = theGlobalNormal.cross(-NewMomentum);
1404             G4ThreeVector vNorm = v/v.mag();
1405             G4ThreeVector u     = vNorm.cross(theGlobalNormal);
1406 
1407             u = u *= (std::sin(elevation) * std::cos(azimuth));
1408             v = vNorm *= (std::sin(elevation) * std::sin(azimuth));
1409             G4ThreeVector w = theGlobalNormal*=(std::cos(elevation));
1410 
1411             NewMomentum = G4ThreeVector(u+v+w);
1412 
1413             // Rotate Polarization too: (needs revision)
1414             NewPolarization = OldPolarization;
1415         }
1416     } 
1417     while (NewMomentum * theGlobalNormal <= 0.0);
1418 }
1419 
1420 void InstrumentedG4OpBoundaryProcess::DielectricDichroic()
1421 {
1422     // Calculate Angle between Normal and Photon Momentum
1423     G4double anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
1424 
1425     // Round it to closest integer
1426     G4double angleIncident = std::floor(180/pi*anglePhotonToNormal+0.5);
1427 
1428     if (!DichroicVector) 
1429     {
1430         if (OpticalSurface) DichroicVector = OpticalSurface->GetDichroicVector();
1431     }
1432 
1433     if (DichroicVector) 
1434     {
1435         G4double wavelength = h_Planck*c_light/thePhotonMomentum;
1436         theTransmittance = DichroicVector->Value(wavelength/nm,angleIncident,idx,idy)*perCent; 
1437 
1438 //      G4cout << "wavelength: " << std::floor(wavelength/nm) 
1439 //                               << "nm" << G4endl;
1440 //      G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
1441 //      G4cout << "Transmittance: " 
1442 //             << std::floor(theTransmittance/perCent) << "%" << G4endl;
1443     } 
1444     else 
1445     {
1446         G4ExceptionDescription ed;
1447         ed << " InstrumentedG4OpBoundaryProcess/DielectricDichroic(): "
1448            << " The dichroic surface has no G4Physics2DVector"
1449            << G4endl;
1450         G4Exception("InstrumentedG4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
1451                     FatalException,ed,
1452                     "A dichroic surface must have an associated G4Physics2DVector");
1453     }
1454 
1455     if ( !G4BooleanRand(theTransmittance) ) 
1456     {   // Not transmitted, so reflect
1457 
1458         if ( theModel == glisur || theFinish == polished ) 
1459         {
1460             DoReflection();
1461         } 
1462         else 
1463         {
1464             ChooseReflection();
1465             if ( theStatus == LambertianReflection ) 
1466             {
1467                 DoReflection();
1468             } 
1469             else if ( theStatus == BackScattering ) 
1470             {
1471                 NewMomentum = -OldMomentum;
1472                 NewPolarization = -OldPolarization;
1473             } 
1474             else 
1475             {
1476                 G4double PdotN, EdotN;
1477                 do 
1478                 {
1479                     if (theStatus==LobeReflection)
1480                         theFacetNormal = GetFacetNormal(OldMomentum,theGlobalNormal);
1481                     PdotN = OldMomentum * theFacetNormal;
1482                     NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1483                     // Loop checking, 13-Aug-2015, Peter Gumplinger
1484                 } while (NewMomentum * theGlobalNormal <= 0.0);
1485                 EdotN = OldPolarization * theFacetNormal;
1486                 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1487             }
1488         }
1489     } 
1490     else 
1491     {
1492         theStatus = Dichroic;
1493         NewMomentum = OldMomentum;
1494         NewPolarization = OldPolarization;
1495     }
1496 }
1497 
1498 
1499 /**
1500 InstrumentedG4OpBoundaryProcess::DielectricDielectric
1501 ------------------------------------------------------
1502 
1503 Sets::
1504 
1505    NewMomentum
1506    NewPolarization
1507    theStatus : TotalInternalReflection/FresnelReflection/FresnelRefraction
1508 
1509 Needs:
1510 
1511 1. theFinish = polished, to use theGlobalNormal
1512 2. theTransmittance > 0 to avoid using the default TransCoeff (so had to add theCustomStatus)
1513 
1514 **/
1515 
1516 
1517 void InstrumentedG4OpBoundaryProcess::DielectricDielectric()
1518 {
1519     G4bool Inside = false;
1520     G4bool Swap = false;
1521 
1522     G4bool SurfaceRoughnessCriterionPass = 1;
1523     if (theSurfaceRoughness != 0. && Rindex1 > Rindex2) 
1524     {
1525         G4double wavelength = h_Planck*c_light/thePhotonMomentum;
1526         G4double SurfaceRoughnessCriterion = std::exp(-std::pow((4*pi*theSurfaceRoughness*Rindex1*cost1/wavelength),2)); 
1527         SurfaceRoughnessCriterionPass = G4BooleanRand(SurfaceRoughnessCriterion); 
1528     }
1529 
1530     leap:
1531 
1532     G4bool Through = false;
1533     G4bool Done = false;
1534 
1535     G4double PdotN, EdotN;
1536 
1537     G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1538     G4double E1_perp, E1_parl;
1539     G4double s1, s2, E2_perp, E2_parl, E2_total, TransCoeff;
1540     G4double E2_abs, C_parl, C_perp;
1541     G4double alpha;
1542 
1543     do 
1544     {
1545         if (Through) 
1546         {
1547             Swap = !Swap;
1548             Through = false;
1549             theGlobalNormal = -theGlobalNormal;
1550             G4SwapPtr(Material1,Material2);
1551             G4SwapObj(&Rindex1,&Rindex2);
1552         }
1553 
1554         theFacetNormal = theFinish == polished ? theGlobalNormal : GetFacetNormal(OldMomentum,theGlobalNormal) ; 
1555 
1556         PdotN = OldMomentum * theFacetNormal;
1557         EdotN = OldPolarization * theFacetNormal;
1558 
1559         cost1 = - PdotN;
1560         if (std::abs(cost1) < 1.0-kCarTolerance)
1561         {
1562             sint1 = std::sqrt(1.-cost1*cost1);
1563             sint2 = sint1*Rindex1/Rindex2;     // *** Snell's Law ***
1564         }
1565         else 
1566         {
1567             sint1 = 0.0;
1568             sint2 = 0.0;
1569         }
1570 
1571 #ifdef DEBUG_PIDX
1572         LOG_IF(LEVEL, pidx_dump)
1573             << "DiDi.pidx " << std::setw(4) << pidx  
1574             << " PIDX " << std::setw(4) << PIDX
1575             << " OldMomentum (" 
1576             << " " << std::setw(10) << std::fixed << std::setprecision(5) << OldMomentum.x()
1577             << " " << std::setw(10) << std::fixed << std::setprecision(5) << OldMomentum.y()
1578             << " " << std::setw(10) << std::fixed << std::setprecision(5) << OldMomentum.z()
1579             << ")" 
1580             << " OldPolarization (" 
1581             << " " << std::setw(10) << std::fixed << std::setprecision(5) << OldPolarization.x()
1582             << " " << std::setw(10) << std::fixed << std::setprecision(5) << OldPolarization.y()
1583             << " " << std::setw(10) << std::fixed << std::setprecision(5) << OldPolarization.z()
1584             << ")" 
1585             << " cost1 " << std::setw(10) << std::fixed << std::setprecision(5) << cost1
1586             << " Rindex1 " << std::setw(10) << std::fixed << std::setprecision(5) << Rindex1
1587             << " Rindex2 " << std::setw(10) << std::fixed << std::setprecision(5) << Rindex2
1588             << " sint1 " << std::setw(10) << std::fixed << std::setprecision(5) << sint1 
1589             << " sint2 " << std::setw(10) << std::fixed << std::setprecision(5) << sint2 
1590             ;    
1591 #endif
1592 
1593         if (sint2 >= 1.0) 
1594         {
1595             //[didi.tir
1596             // Simulate total internal reflection
1597             if (Swap) Swap = !Swap;
1598 
1599             theStatus = TotalInternalReflection;
1600             if ( !SurfaceRoughnessCriterionPass ) theStatus = LambertianReflection; 
1601 
1602             if ( theModel == unified && theFinish != polished ) ChooseReflection(); 
1603 
1604             if ( theStatus == LambertianReflection ) 
1605             {
1606                 DoReflection();
1607             }
1608             else if ( theStatus == BackScattering ) 
1609             {
1610                 NewMomentum = -OldMomentum;
1611                 NewPolarization = -OldPolarization;
1612             }
1613             else 
1614             {
1615                 PdotN = OldMomentum * theFacetNormal;
1616                 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1617                 EdotN = OldPolarization * theFacetNormal;
1618                 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1619             }
1620             //]didi.tir
1621         }
1622         else if (sint2 < 1.0) 
1623         {
1624 
1625             // Calculate amplitude for transmission (Q = P x N)
1626 
1627             if (cost1 > 0.0) 
1628             {
1629                 cost2 =  std::sqrt(1.-sint2*sint2);
1630             }
1631             else 
1632             {
1633                 cost2 = -std::sqrt(1.-sint2*sint2);
1634             }
1635 
1636             if (sint1 > 0.0)
1637             {
1638 #ifdef DEBUG_PIDX
1639                 if(pidx_dump) printf("//DiDi pidx %6d : sint1 > 0 \n", pidx );  
1640 #endif
1641 
1642                 A_trans = OldMomentum.cross(theFacetNormal);
1643                 A_trans = A_trans.unit();
1644                 E1_perp = OldPolarization * A_trans;
1645                 E1pp    = E1_perp * A_trans;
1646                 E1pl    = OldPolarization - E1pp;
1647                 E1_parl = E1pl.mag();
1648             }
1649             else 
1650             {
1651 #ifdef DEBUG_PIDX
1652                 LOG_IF(LEVEL, pidx_dump)
1653                    << " pidx " << pidx 
1654                    << " NOT:sint1 > 0 : JACKSON NORMAL INCIDENCE "
1655                    ;
1656 #endif
1657 
1658                 A_trans  = OldPolarization;
1659                 // Here we Follow Jackson's conventions and we set the
1660                 // parallel component = 1 in case of a ray perpendicular
1661                 // to the surface
1662                 E1_perp  = 0.0;
1663                 E1_parl  = 1.0;
1664             }
1665 
1666             s1 = Rindex1*cost1;
1667             E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
1668             E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
1669             E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1670             s2 = Rindex2*cost2*E2_total;
1671 
1672             if (theTransmittance > 0 || theCustomStatus == 'Y') TransCoeff = theTransmittance;
1673             else if (cost1 != 0.0) TransCoeff = s2/s1;
1674             else TransCoeff = 0.0;
1675 
1676 #ifdef DEBUG_PIDX
1677             LOG_IF(LEVEL, pidx_dump)
1678                  << " pidx " << pidx 
1679                  << " TransCoeff " << TransCoeff
1680                  ;
1681 #endif
1682             if ( !G4BooleanRand_TransCoeff(TransCoeff) ) 
1683             {
1684                 // Simulate reflection
1685                 if (Swap) Swap = !Swap;
1686 
1687                 theStatus = FresnelReflection;
1688 
1689                 if ( !SurfaceRoughnessCriterionPass ) theStatus = LambertianReflection; 
1690 
1691                 if ( theModel == unified && theFinish != polished ) ChooseReflection(); 
1692 
1693                 if ( theStatus == LambertianReflection ) 
1694                 {
1695                     DoReflection();
1696                 }
1697                 else if ( theStatus == BackScattering ) 
1698                 {
1699                     NewMomentum = -OldMomentum;
1700                     NewPolarization = -OldPolarization;
1701                 }
1702                 else 
1703                 {
1704                     PdotN = OldMomentum * theFacetNormal;
1705                     NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1706 
1707                     if (sint1 > 0.0) 
1708                     {   // incident ray oblique
1709 
1710                         E2_parl   = Rindex2*E2_parl/Rindex1 - E1_parl;
1711                         E2_perp   = E2_perp - E1_perp;
1712                         E2_total  = E2_perp*E2_perp + E2_parl*E2_parl;
1713                         A_paral   = NewMomentum.cross(A_trans);
1714                         A_paral   = A_paral.unit();
1715                         E2_abs    = std::sqrt(E2_total);
1716                         C_parl    = E2_parl/E2_abs;
1717                         C_perp    = E2_perp/E2_abs;
1718 
1719                         NewPolarization = C_parl*A_paral + C_perp*A_trans;
1720                     }
1721                     else 
1722                     {     // incident ray perpendicular
1723 
1724                         if (Rindex2 > Rindex1) 
1725                         {
1726                             NewPolarization = - OldPolarization;
1727                         }
1728                         else 
1729                         {
1730                             NewPolarization =   OldPolarization;
1731                         }
1732                     }
1733                 }
1734             }
1735             else 
1736             {   // photon gets transmitted
1737                 // Simulate transmission/refraction
1738 #ifdef DEBUG_PIDX
1739                 LOG_IF(LEVEL, pidx_dump) 
1740                     << " pidx " << pidx 
1741                     << " TRANSMIT "
1742                     ; 
1743 #endif
1744                 Inside = !Inside;
1745                 Through = true;
1746                 theStatus = FresnelRefraction;
1747 
1748                 if (sint1 > 0.0) 
1749                 {      // incident ray oblique
1750 
1751                     alpha = cost1 - cost2*(Rindex2/Rindex1);
1752                     NewMomentum = OldMomentum + alpha*theFacetNormal;
1753                     NewMomentum = NewMomentum.unit();
1754 //                      PdotN = -cost2;
1755                     A_paral = NewMomentum.cross(A_trans);
1756                     A_paral = A_paral.unit();
1757                     E2_abs     = std::sqrt(E2_total);
1758                     C_parl     = E2_parl/E2_abs;
1759                     C_perp     = E2_perp/E2_abs;
1760 
1761                     NewPolarization = C_parl*A_paral + C_perp*A_trans;
1762                 }
1763                 else  
1764                 {    // incident ray perpendicular
1765                     NewMomentum = OldMomentum;
1766                     NewPolarization = OldPolarization;
1767                 }
1768 #ifdef DEBUG_PIDX
1769                 LOG_IF(LEVEL, pidx_dump)
1770                      << " pidx " << pidx 
1771                      << " TRANSMIT "
1772                      << " NewMom "
1773                      <<  "(" 
1774                      << " " << NewMomentum.x()
1775                      << " " << NewMomentum.y()
1776                      << " " << NewMomentum.z()
1777                      << ")" 
1778                      << " NewPol "
1779                      <<  "(" 
1780                      << " " << NewPolarization.x()
1781                      << " " << NewPolarization.y()
1782                      << " " << NewPolarization.z()
1783                      << ")" 
1784                      ;
1785                 //
1786 #endif
1787             }
1788         }
1789 
1790         OldMomentum = NewMomentum.unit();
1791         OldPolarization = NewPolarization.unit();
1792 
1793         if (theStatus == FresnelRefraction) 
1794         {
1795             Done = (NewMomentum * theGlobalNormal <= 0.0);
1796         } 
1797         else 
1798         {
1799             Done = (NewMomentum * theGlobalNormal >= -kCarTolerance);
1800         }
1801 
1802         // Loop checking, 13-Aug-2015, Peter Gumplinger
1803     } while (!Done);
1804 
1805     if (Inside && !Swap) 
1806     {
1807         if( theFinish == polishedbackpainted ||  theFinish == groundbackpainted ) 
1808         {
1809 
1810             G4double rand = G4UniformRand();
1811             if ( rand > theReflectivity ) 
1812             {
1813                 if (rand > theReflectivity + theTransmittance) 
1814                 {
1815                     DoAbsorption();
1816                 } 
1817                 else 
1818                 {
1819                     theStatus = Transmission;
1820                     NewMomentum = OldMomentum;
1821                     NewPolarization = OldPolarization;
1822                 }
1823             }
1824             else 
1825             {
1826                 if (theStatus != FresnelRefraction )
1827                 {
1828                     theGlobalNormal = -theGlobalNormal;
1829                 }
1830                 else 
1831                 {
1832                     Swap = !Swap;
1833                     G4SwapPtr(Material1,Material2);
1834                     G4SwapObj(&Rindex1,&Rindex2);
1835                 }
1836                 if ( theFinish == groundbackpainted ) theStatus = LambertianReflection; 
1837 
1838                 DoReflection();
1839 
1840                 theGlobalNormal = -theGlobalNormal;
1841                 OldMomentum = NewMomentum;
1842 
1843                 goto leap;
1844             }
1845         }
1846     }
1847 }
1848 
1849 
1850 /**
1851 InstrumentedG4OpBoundaryProcess::GetIncidentAngle
1852 ---------------------------------------------------
1853 
1854 SCB: HUH, both magP and magN should be 1 ? 
1855 SCB; also why bother to acos, what you really need is cos_angle and sin_angle
1856 
1857 **/
1858 
1859 G4double InstrumentedG4OpBoundaryProcess::GetIncidentAngle() 
1860 {
1861     G4double PdotN = OldMomentum * theFacetNormal;
1862     G4double magP= OldMomentum.mag();
1863     G4double magN= theFacetNormal.mag();
1864 
1865     G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
1866 
1867     return incidentangle;
1868 }
1869 
1870 /**
1871 InstrumentedG4OpBoundaryProcess::GetReflectivity
1872 --------------------------------------------------
1873 
1874 Consumes randoms in do-while to set iTE, iTM
1875 
1876 **/
1877 
1878 G4double InstrumentedG4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1879                                               G4double E1_parl,
1880                                               G4double incidentangle,
1881                                               G4double RealRindex,
1882                                               G4double ImaginaryRindex)
1883 {
1884     G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1885     G4complex N1(Rindex1, 0), N2(RealRindex, ImaginaryRindex);
1886     G4complex CosPhi;
1887 
1888     G4complex u(1,0);           //unit number 1
1889 
1890     G4complex numeratorTE;      // E1_perp=1 E1_parl=0 -> TE polarization
1891     G4complex numeratorTM;      // E1_parl=1 E1_perp=0 -> TM polarization
1892     G4complex denominatorTE, denominatorTM;
1893     G4complex rTM, rTE;
1894 
1895     G4MaterialPropertiesTable* aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable(); 
1896     G4MaterialPropertyVector* aPropertyPointerR = aMaterialPropertiesTable->GetProperty(kREALRINDEX); 
1897     G4MaterialPropertyVector* aPropertyPointerI = aMaterialPropertiesTable->GetProperty(kIMAGINARYRINDEX); 
1898 
1899     if (aPropertyPointerR && aPropertyPointerI) 
1900     {
1901         G4double RRindex = aPropertyPointerR->Value(thePhotonMomentum);
1902         G4double IRindex = aPropertyPointerI->Value(thePhotonMomentum);
1903         N1 = G4complex(RRindex,IRindex);
1904     }
1905 
1906     // Following two equations, rTM and rTE, are from: "Introduction To Modern
1907     // Optics" written by Fowles
1908 
1909     CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))*(N1*N1)/(N2*N2)));
1910 
1911     numeratorTE   = N1*std::cos(incidentangle) - N2*CosPhi;
1912     denominatorTE = N1*std::cos(incidentangle) + N2*CosPhi;
1913     rTE = numeratorTE/denominatorTE;
1914 
1915     numeratorTM   = N2*std::cos(incidentangle) - N1*CosPhi;
1916     denominatorTM = N2*std::cos(incidentangle) + N1*CosPhi;
1917     rTM = numeratorTM/denominatorTM;
1918 
1919     // This is my calculaton for reflectivity on a metalic surface
1920     // depending on the fraction of TE and TM polarization
1921     // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1922     // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1923 
1924     Reflectivity_TE =  (rTE*conj(rTE))*(E1_perp*E1_perp)
1925                      / (E1_perp*E1_perp + E1_parl*E1_parl);
1926     Reflectivity_TM =  (rTM*conj(rTM))*(E1_parl*E1_parl)
1927                     / (E1_perp*E1_perp + E1_parl*E1_parl);
1928     Reflectivity    = Reflectivity_TE + Reflectivity_TM;
1929 
1930     do {
1931         if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1932         {iTE = -1;}else{iTE = 1;}
1933         if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1934         {iTM = -1;}else{iTM = 1;}
1935         // Loop checking, 13-Aug-2015, Peter Gumplinger
1936     } while(iTE<0&&iTM<0);
1937 
1938     return real(Reflectivity);
1939 }
1940 
1941 /**
1942 InstrumentedG4OpBoundaryProcess::CalculateReflectivity
1943 ---------------------------------------------------------
1944 
1945 UGLY: 
1946 
1947 * depends on PropertyPointer1, PropertyPointer2 for Real and Imag Rindex
1948 * sets theReflectivity 
1949 
1950 **/
1951 
1952 void InstrumentedG4OpBoundaryProcess::CalculateReflectivity()
1953 {
1954     G4double RealRindex = PropertyPointer1->Value(thePhotonMomentum); 
1955     G4double ImaginaryRindex = PropertyPointer2->Value(thePhotonMomentum); 
1956 
1957     theFacetNormal = theFinish == ground ? GetFacetNormal(OldMomentum, theGlobalNormal) : theGlobalNormal ;
1958 
1959     G4double PdotN = OldMomentum * theFacetNormal;
1960     cost1 = -PdotN;
1961     sint1 = (std::abs(cost1) < 1.0 - kCarTolerance) ? std::sqrt(1. - cost1*cost1) : 0.0 ;
1962 
1963     G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1964     G4double E1_perp, E1_parl;
1965 
1966     if (sint1 > 0.0 ) 
1967     {
1968         A_trans = OldMomentum.cross(theFacetNormal);
1969         A_trans = A_trans.unit();
1970         E1_perp = OldPolarization * A_trans;
1971         E1pp    = E1_perp * A_trans;
1972         E1pl    = OldPolarization - E1pp;
1973         E1_parl = E1pl.mag();
1974     }
1975     else 
1976     {
1977         A_trans  = OldPolarization;
1978         // Here we Follow Jackson's conventions and we set the
1979         // parallel component = 1 in case of a ray perpendicular
1980         // to the surface
1981         E1_perp  = 0.0;
1982         E1_parl  = 1.0;
1983     }
1984 
1985     //calculate incident angle
1986     G4double incidentangle = GetIncidentAngle();
1987 
1988     //calculate the reflectivity depending on incident angle,
1989     //polarization and complex refractive
1990 
1991     theReflectivity = GetReflectivity(E1_perp, E1_parl, incidentangle, RealRindex, ImaginaryRindex); 
1992 }
1993 
1994 
1995 /**
1996 InstrumentedG4OpBoundaryProcess::DoAbsorption
1997 ----------------------------------------------
1998 
1999 Sets: 
2000 
2001 * theStatus to Detection/Absorption depending on one random and theEfficiency 
2002 * sets NewMomentum, NewPolarization to the Old 
2003 * aParticleChange 
2004 
2005 **/
2006 
2007 void InstrumentedG4OpBoundaryProcess::DoAbsorption()
2008 {
2009     LOG(LEVEL) 
2010         << " PostStepDoIt_count " << PostStepDoIt_count
2011         << " theEfficiency " << theEfficiency
2012         ;
2013 
2014     bool detect = G4BooleanRand_theEfficiency(theEfficiency) ; 
2015     theStatus = detect ? Detection : Absorption ; 
2016 
2017     NewMomentum = OldMomentum;
2018     NewPolarization = OldPolarization;
2019 
2020     aParticleChange.ProposeLocalEnergyDeposit(detect ? thePhotonMomentum : 0.0);
2021     aParticleChange.ProposeTrackStatus(fStopAndKill);
2022 }
2023 
2024 
2025 /**
2026 InstrumentedG4OpBoundaryProcess::DoReflection
2027 -----------------------------------------------
2028 
2029 Sets:: 
2030 
2031    NewMomentum
2032    NewPolarization
2033    theFacetNormal
2034    theStatus (unless already LambertianReflection)
2035 
2036 
2037 theStatus:LambertianReflection
2038     NewMomentum Lambertian sampled around theGlobalNormal, changes theFacetNormal  
2039 
2040 theFinish:ground
2041     sets theStatus:LobeReflection
2042 
2043 
2044 
2045            B
2046           /|\
2047          / | \
2048         /  |  \
2049        /   |   \
2050       + - -+- - +      
2051        \   |   /
2052        OLD |  NEW
2053          \ | / 
2054           \|/
2055   ---------A------------
2056 
2057     A->B = NEW - OLD
2058 
2059 **/
2060 
2061 void InstrumentedG4OpBoundaryProcess::DoReflection()
2062 {
2063     if( theStatus == LambertianReflection ) 
2064     {
2065         NewMomentum = U4LambertianRand(theGlobalNormal);
2066         theFacetNormal = (NewMomentum - OldMomentum).unit();
2067     }
2068     else if( theFinish == ground ) 
2069     {
2070         theStatus = LobeReflection;
2071         if (PropertyPointer1 && PropertyPointer2)
2072         {
2073         } 
2074         else 
2075         {
2076             theFacetNormal = GetFacetNormal(OldMomentum,theGlobalNormal); 
2077         }
2078         G4double PdotN = OldMomentum * theFacetNormal;
2079         NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
2080    }
2081    else 
2082    {
2083        theStatus = SpikeReflection;
2084        theFacetNormal = theGlobalNormal;
2085        G4double PdotN = OldMomentum * theFacetNormal;
2086        NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
2087    }
2088    G4double EdotN = OldPolarization * theFacetNormal;
2089    NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
2090 }
2091 
2092 
2093 /**
2094 InstrumentedG4OpBoundaryProcess::ChooseReflection
2095 --------------------------------------------------
2096 
2097 Sets theStatus to SpikeReflection/LobeReflection/BackScattering/LambertianReflection
2098 depending on one rand and prob_ss/prob_sl/prob_bs. 
2099 
2100 But prob_ss/prob_sl/prob_bs all default to zero and require unified model and appropriate 
2101 properties in OpticalSurface.mpt to change from zero. So this will typically return LambertianReflection
2102 
2103 For SpikeReflection theFacetNormal is set to theGlobalNormal
2104 **/
2105 
2106 void InstrumentedG4OpBoundaryProcess::ChooseReflection()
2107 {
2108     G4double rand = G4UniformRand();
2109 #ifndef PRODUCTION
2110 #ifdef DEBUG_TAG
2111     SEvt::AddTag( 1, U4Stack_ChooseReflection, rand ); 
2112 #endif
2113 #endif
2114 
2115     if ( rand >= 0.0 && rand < prob_ss ) 
2116     {
2117         theStatus = SpikeReflection;
2118         theFacetNormal = theGlobalNormal;
2119     }
2120     else if ( rand >= prob_ss && rand <= prob_ss+prob_sl) 
2121     {
2122         theStatus = LobeReflection;
2123     }
2124     else if ( rand > prob_ss+prob_sl && rand < prob_ss+prob_sl+prob_bs ) 
2125     {
2126         theStatus = BackScattering;
2127     }
2128     else 
2129     {
2130         theStatus = LambertianReflection;
2131     }
2132 
2133     LOG(LEVEL) 
2134         << " theStatus " << U4OpBoundaryProcessStatus::Name(theStatus) 
2135         << " prob_ss " << prob_ss 
2136         << " prob_sl " << prob_sl 
2137         << " prob_bs " << prob_bs 
2138         ; 
2139 }
2140 
2141 
2142 
2143 
2144 G4bool InstrumentedG4OpBoundaryProcess::IsApplicable(const G4ParticleDefinition& aParticleType)
2145 {
2146     return ( &aParticleType == G4OpticalPhoton::OpticalPhoton() );
2147 }
2148 
2149 G4OpBoundaryProcessStatus InstrumentedG4OpBoundaryProcess::GetStatus() const
2150 {
2151     return theStatus;
2152 }
2153 
2154 void InstrumentedG4OpBoundaryProcess::SetInvokeSD(G4bool flag)
2155 {
2156     fInvokeSD = flag;
2157 }
2158 
2159 
2160 G4bool InstrumentedG4OpBoundaryProcess::G4BooleanRand(const G4double prob) const
2161 {
2162     G4double u = G4UniformRand() ; 
2163     G4bool ret = u < prob  ; 
2164 #ifdef DEBUG_PIDX
2165     LOG_IF(LEVEL, pidx_dump)
2166          << " pidx " << pidx
2167          << " prob " << prob 
2168          << " ret " << ret 
2169          ;   
2170 #endif
2171     return ret ; 
2172 }
2173 
2174 G4bool InstrumentedG4OpBoundaryProcess::G4BooleanRand_TransCoeff(const G4double prob) const
2175 {
2176     G4double u = G4UniformRand() ; 
2177     G4bool ret = u < prob  ; 
2178 
2179     LOG(LEVEL) 
2180         << U4UniformRand::Desc(u, SEvt::UU)
2181         << " TransCoeff " << prob 
2182         << " DECISION " << ( ret ? "T" : "R" ) 
2183         ; 
2184 #ifndef PRODUCTION
2185 #ifdef DEBUG_TAG
2186     SEvt::AddTag(1, U4Stack_BoundaryDiDiTransCoeff, u ); 
2187 #endif   
2188 #endif
2189 
2190 #ifndef PRODUCTION
2191 #ifdef DEBUG_PIDX
2192     LOG_IF(LEVEL, pidx_dump)
2193          << " pidx " << pidx
2194          << " prob " << prob 
2195          << " ret " << ret 
2196          ;   
2197 #endif
2198 #endif
2199     return ret ; 
2200 }
2201 
2202 G4bool InstrumentedG4OpBoundaryProcess::G4BooleanRand_theEfficiency(const G4double prob) const
2203 {
2204     G4double u = G4UniformRand() ; 
2205     G4bool ret = u < prob  ; 
2206 #ifndef PRODUCTION
2207 #ifdef DEBUG_TAG
2208     SEvt::AddTag(1, U4Stack_AbsorptionEffDetect, u ); 
2209 #endif   
2210 #endif   
2211 #ifndef PRODUCTION
2212 #ifdef DEBUG_PIDX
2213     LOG_IF(LEVEL, pidx_dump)
2214          << " pidx " << pidx
2215          << " prob " << prob 
2216          << " ret " << ret 
2217          ;   
2218 #endif
2219 #endif
2220     return ret ; 
2221 }
2222 
2223 G4bool InstrumentedG4OpBoundaryProcess::G4BooleanRand_theReflectivity(const G4double prob) const
2224 {
2225     G4double u = G4UniformRand() ; 
2226     G4bool ret = u < prob ; 
2227 #ifdef DEBUG_PIDX
2228     LOG_IF(LEVEL, pidx_dump)
2229         << " pidx " << pidx 
2230         << " prob " << prob
2231         << " u " << u 
2232         << " ret " << ret
2233         ;
2234 
2235 #endif
2236     return ret  ; 
2237 }
2238 
2239 
2240 G4bool InstrumentedG4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
2241 {
2242     G4Step aStep = *pStep;
2243     aStep.AddTotalEnergyDeposit(thePhotonMomentum);
2244     G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
2245     return sd ? sd->Hit(&aStep) : false ;
2246 }
2247 
2248 G4double InstrumentedG4OpBoundaryProcess::GetMeanFreePath(const G4Track& , G4double , G4ForceCondition* condition) 
2249 {
2250     *condition = Forced;
2251     return DBL_MAX;
2252 }
2253 
2254 
2255 void InstrumentedG4OpBoundaryProcess::BoundaryProcessVerbose() const
2256 {
2257         if ( theStatus == Undefined )
2258                 G4cout << " *** Undefined *** " << G4endl;
2259         if ( theStatus == Transmission )
2260                 G4cout << " *** Transmission *** " << G4endl;
2261         if ( theStatus == FresnelRefraction )
2262                 G4cout << " *** FresnelRefraction *** " << G4endl;
2263         if ( theStatus == FresnelReflection )
2264                 G4cout << " *** FresnelReflection *** " << G4endl;
2265         if ( theStatus == TotalInternalReflection )
2266                 G4cout << " *** TotalInternalReflection *** " << G4endl;
2267         if ( theStatus == LambertianReflection )
2268                 G4cout << " *** LambertianReflection *** " << G4endl;
2269         if ( theStatus == LobeReflection )
2270                 G4cout << " *** LobeReflection *** " << G4endl;
2271         if ( theStatus == SpikeReflection )
2272                 G4cout << " *** SpikeReflection *** " << G4endl;
2273         if ( theStatus == BackScattering )
2274                 G4cout << " *** BackScattering *** " << G4endl;
2275         if ( theStatus == PolishedLumirrorAirReflection )
2276                 G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
2277         if ( theStatus == PolishedLumirrorGlueReflection )
2278                 G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
2279         if ( theStatus == PolishedAirReflection )
2280                 G4cout << " *** PolishedAirReflection *** " << G4endl;
2281         if ( theStatus == PolishedTeflonAirReflection )
2282                 G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
2283         if ( theStatus == PolishedTiOAirReflection )
2284                 G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
2285         if ( theStatus == PolishedTyvekAirReflection )
2286                 G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
2287         if ( theStatus == PolishedVM2000AirReflection )
2288                 G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
2289         if ( theStatus == PolishedVM2000GlueReflection )
2290                 G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
2291         if ( theStatus == EtchedLumirrorAirReflection )
2292                 G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
2293         if ( theStatus == EtchedLumirrorGlueReflection )
2294                 G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
2295         if ( theStatus == EtchedAirReflection )
2296                 G4cout << " *** EtchedAirReflection *** " << G4endl;
2297         if ( theStatus == EtchedTeflonAirReflection )
2298                 G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
2299         if ( theStatus == EtchedTiOAirReflection )
2300                 G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
2301         if ( theStatus == EtchedTyvekAirReflection )
2302                 G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
2303         if ( theStatus == EtchedVM2000AirReflection )
2304                 G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
2305         if ( theStatus == EtchedVM2000GlueReflection )
2306                 G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
2307         if ( theStatus == GroundLumirrorAirReflection )
2308                 G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
2309         if ( theStatus == GroundLumirrorGlueReflection )
2310                 G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
2311         if ( theStatus == GroundAirReflection )
2312                 G4cout << " *** GroundAirReflection *** " << G4endl;
2313         if ( theStatus == GroundTeflonAirReflection )
2314                 G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
2315         if ( theStatus == GroundTiOAirReflection )
2316                 G4cout << " *** GroundTiOAirReflection *** " << G4endl;
2317         if ( theStatus == GroundTyvekAirReflection )
2318                 G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
2319         if ( theStatus == GroundVM2000AirReflection )
2320                 G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
2321         if ( theStatus == GroundVM2000GlueReflection )
2322                 G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
2323         if ( theStatus == Absorption )
2324                 G4cout << " *** Absorption *** " << G4endl;
2325         if ( theStatus == Detection )
2326                 G4cout << " *** Detection *** " << G4endl;
2327         if ( theStatus == NotAtBoundary )
2328                 G4cout << " *** NotAtBoundary *** " << G4endl;
2329         if ( theStatus == SameMaterial )
2330                 G4cout << " *** SameMaterial *** " << G4endl;
2331         if ( theStatus == StepTooSmall )
2332                 G4cout << " *** StepTooSmall *** " << G4endl;
2333         if ( theStatus == NoRINDEX )
2334                 G4cout << " *** NoRINDEX *** " << G4endl;
2335         if ( theStatus == Dichroic )
2336                 G4cout << " *** Dichroic Transmission *** " << G4endl;
2337 }
2338 
2339
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!