eic-opticks/u4/CustomBoundary.h

0001 #pragma once
0002 /**
0003 CustomBoundary.h : Adding Multi-Layer TMM ARTD to G4OpBoundaryProcess
0004 =======================================================================
0005
0006 Alternative names::
0007
0008     CustomSurface.h
0009     TMMSurface.h
0010     MultiLayerSurface.h
0011     MultiFilmSurface.h
0012
0013 Incorporating "Surface" in the name seems more appropriate than Boundary
0014 as its acting as a "smart" surface from Geant4 point of view.
0015
0016
0017 This extends G4OpBoundaryProcess providing custom calculation of ARTD (Absorb,
0018 Reflect, Transmit, Detect) probabilities.  This can be used for example to add
0019 multi-layer TMM calculation of ARTD based on layer refractive indices and
0020 thicknesses.
0021
0022 For example this could implement coherent interference effects for a PMT
0023 on the boundary between Pyrex and Vacuum where thin layers of ARC
0024 (anti-reflection coating) and PHC (photocathode) are located::
0025
0026          +---------------pmt-Pyrex----------------+
0027          |                                        |
0028          |                                        |
0029          |                                        |
0030          |       +~inner~Vacuum~~~~~~~~~~~+       |
0031          |       !                        !       |
0032          |       !                        !       |
0033          |       !                        !       |
0034          |       !                        !       |
0035          |       !                        !       |
0036          |       +                        +       |
0037          |       |                        |       |
0038          |       |                        |       |
0039          |       |                        |       |
0040          |       |                        |       |
0041          |       |                        |       |
0042          |       +------------------------+       |
0043          |                                        |
0044          |                                        |
0045          |                                        |
0046          +----------------------------------------+
0047
0048 The parameter names follow those of G4OpBoundaryProcess apart from
0049 *theRecoveredNormal*.  *theRecoveredNormal* is required to be the outwards
0050 normal for the inner volume in global frame. Absolutely no G4Track dependent
0051 flips must be applied to this normal, it must be purely an outwards geometry
0052 normal. NB this means that theGlobalExitNormal must be corrected, removing a
0053 flip done by G4Navigator based on "enteredDaughter".
0054
0055 For example, from a modified G4OpBoundaryProcess::
0056
0057      418     G4bool haveEnteredDaughter= (thePostPV->GetMotherLogical() == thePrePV ->GetLogicalVolume()); // SCB
0058      419
0059      ...
0060      449     G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
0061      450     std::vector<G4Navigator*>::iterator iNav =
0062      451                 G4TransportationManager::GetTransportationManager()->
0063      452                                          GetActiveNavigatorsIterator();
0064      453     theGlobalExitNormal =
0065      454                    (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
0066      455
0067      456     // theGlobalExitNormal is already oriented by G4Navigator to point from vol1 -> vol2
0068      457     // so undo that flip by G4Navigator in order to recover the original geometry
0069      458     // normal that is independent of G4Track direction
0070      459
0071      460     theRecoveredNormal = ( haveEnteredDaughter ? -1. : 1. )* theGlobalExitNormal  ;
0072      464     theGlobalNormal = theGlobalExitNormal ;
0073
0074 CAUTION: theGlobalNormal is compromised (by being flipped) multiple times
0075 within G4OpBoundaryProcess so it is vital that a separate *theGlobalExitNormal*
0076 is used which is never compromised.
0077
0078 Requiring *theRecoveredNormal* allows the *doIt* to be implemented without
0079 having to transform into local frame, because it provides an absolute
0080 orientation. If this approach turns out to be inconvenient another way of
0081 establishing an absolute orientation is to take a compromised normal (that has
0082 been flipped depending on G4Track an unknown number of times) and transform
0083 that into the local frame and then fix the orientation based on the sign of the
0084 z-component of the local normal::
0085
0086     G4ThreeVector localNormal = transform.TransformAxis(theGlobalNormal);
0087     bool outwards = localNormal.z() > 0. ; // as always upper hemi of PMT in local frame
0088     G4ThreeVector surface_normal = (outwards ? 1. : -1.)*localNormal ;
0089
0090 The disadvantage of using a local normal is that it then makes it necessary to
0091 transform momentum and polarization into local frame and transform back to
0092 global frame.  Conversely using *theRecoveredNormal* which is an absolute
0093 global normal allows the whole calculation to be done in global frame, avoiding
0094 all that transforming.
0095
0096
0097 CustomBoundary::maybe_doIt
0098     for OpticalSurfaceName starting with '@' local_z is used to
0099     check for upper half (eg upper hemisphere of ellipsoidal PMT) in which
0100     case the doIt is run. The local_z is obtained from theGlobalPoint transformed
0101     into the frame of the G4Track.
0102
0103 CustomBoundary::doIt
0104     performs its calculations purely in global frame,
0105     using the dot product between *theRecoveredNormal* and *OldMomentum* to provide the
0106     orientation of the photon with respect to the boundary and the *minus_cos_theta*
0107     where theta is the angle of incidence.
0108
0109 template type J
0110     required to have enumerations : RINDEX, KINDEX, L0, L1, L2, L3, DEFAULT_CAT
0111     and methods to access layer refractive indices and thicknesses for various
0112     categories of sensors (eg PMTs) and layers with arguments including energy_eV
0113
0114     * currently required to have argumentless ctor
0115
0116
0117 WIP : polarization comparison sysrap/tests/stmm_vs_sboundary_test.sh
0118 ----------------------------------------------------------------------
0119
0120 The polarization on reflection copied from junoPMTOpticalModel is for TIR
0121 (not ordinary Fresnel Reflection)
0122
0123 WIP : Can standard DielectricDielectric be used with TMM calc ?
0124 --------------------------------------------------------------------------------
0125
0126 Trying to reuse the standard mom/pol boundary impl together with
0127 the overidden calculation of the below, done in CustomART::
0128
0129    theTransmittance
0130    theReflectivity
0131    theEfficiency
0132
0133 * will require restriction to specific theModel and theFinish
0134   values to avoid the code wandering off the needed path
0135
0136 WIP : Detector Specific property rindex + thickness + qe access
0137 ------------------------------------------------------------------
0138
0139 * maybe make template type J accept an MPT MaterialPropertiesTable in its ctor, enabling
0140   it to be implemented purely from a bunch of custom properties loading into the
0141   OpticalSurface : then it would be more convenient for the instance of J to
0142   be passed in as a parameter to CustomBoundary.
0143
0144 **/
0145
0146
0147 #include "G4ThreeVector.hh"
0148 #include "Randomize.hh"
0149
0150 #include "SLOG.hh"
0151 #include "JPMT.h"
0152 #include "Layr.h"
0153 #include "U4Touchable.h"
0154
0155 #include "CustomStatus.h"
0156
0157
0158 #ifdef DEBUG_TAG
0159 #include "U4UniformRand.h"
0160 #include "SPhoton_Debug.h"
0161 template<> std::vector<SPhoton_Debug<'B'>> SPhoton_Debug<'B'>::record = {} ;
0162 #endif
0163
0164
0165
0166
0167
0168 template<typename J>
0169 struct CustomBoundary
0170 {
0171     static const constexpr plog::Severity LEVEL = debug ;
0172     static void Save(const char* fold);
0173
0174     J*     parameter_accessor ;
0175     int    count ;
0176     double zlocal ;
0177     double lposcost ;
0178
0179     G4ThreeVector& NewMomentum ;
0180     G4ThreeVector& NewPolarization ;
0181     G4ParticleChange& aParticleChange ;
0182     G4OpBoundaryProcessStatus& theStatus ;
0183
0184     const G4ThreeVector& theGlobalPoint ;
0185     const G4ThreeVector& OldMomentum ;
0186     const G4ThreeVector& OldPolarization ;
0187     const G4ThreeVector& theRecoveredNormal ;
0188     const G4double& thePhotonMomentum ;
0189
0190     CustomBoundary(
0191         G4ThreeVector& NewMomentum,
0192         G4ThreeVector& NewPolarization,
0193         G4ParticleChange& aParticleChange,
0194         G4OpBoundaryProcessStatus& theStatus,
0195         const G4ThreeVector& theGlobalPoint,
0196         const G4ThreeVector& OldMomentum,
0197         const G4ThreeVector& OldPolarization,
0198         const G4ThreeVector& theRecoveredNormal,
0199         const G4double& thePhotonMomentum
0200       );
0201
0202     char maybe_doIt(const char* OpticalSurfaceName, const G4Track& aTrack, const G4Step& aStep) ;
0203     char doIt(const G4Track& aTrack, const G4Step& aStep );
0204 };
0205
0206
0207 template<typename J>
0208 inline void CustomBoundary<J>::Save(const char* fold) // static
0209 {
0210     SPhoton_Debug<'B'>::Save(fold);
0211 }
0212
0213 template<typename J>
0214 inline CustomBoundary<J>::CustomBoundary(
0215           G4ThreeVector& NewMomentum_,
0216           G4ThreeVector& NewPolarization_,
0217           G4ParticleChange& aParticleChange_,
0218           G4OpBoundaryProcessStatus& theStatus_,
0219     const G4ThreeVector& theGlobalPoint_,
0220     const G4ThreeVector& OldMomentum_,
0221     const G4ThreeVector& OldPolarization_,
0222     const G4ThreeVector& theRecoveredNormal_,
0223     const G4double&      thePhotonMomentum_
0224     )
0225     :
0226     parameter_accessor(new J),
0227     count(0),
0228     zlocal(-1.),
0229     lposcost(-2.),
0230     NewMomentum(NewMomentum_),
0231     NewPolarization(NewPolarization_),
0232     aParticleChange(aParticleChange_),
0233     theStatus(theStatus_),
0234     theGlobalPoint(theGlobalPoint_),
0235     OldMomentum(OldMomentum_),
0236     OldPolarization(OldPolarization_),
0237     theRecoveredNormal(theRecoveredNormal_),
0238     thePhotonMomentum(thePhotonMomentum_)
0239 {
0240 }
0241
0242 template<typename J>
0243 inline char CustomBoundary<J>::maybe_doIt(const char* OpticalSurfaceName, const G4Track& aTrack, const G4Step& aStep )
0244 {
0245     if( OpticalSurfaceName == nullptr || OpticalSurfaceName[0] != '@') return 'N' ;  // N:OpticalSurfaceName does not start with '@'
0246
0247     const G4AffineTransform& transform = aTrack.GetTouchable()->GetHistory()->GetTopTransform();
0248     G4ThreeVector localPoint = transform.TransformPoint(theGlobalPoint);
0249     zlocal = localPoint.z() ;
0250     lposcost = localPoint.cosTheta() ;
0251
0252     if( zlocal <= 0. ) return 'Z' ;    // Z:-ve Z not triggered
0253
0254     char status = doIt(aTrack, aStep) ;
0255     assert( strchr("ARTD", status) != nullptr ) ;
0256
0257     return status ;
0258 }
0259
0260
0261 template<typename J>
0262 inline char CustomBoundary<J>::doIt(const G4Track& aTrack, const G4Step& aStep )
0263 {
0264     const G4ThreeVector& surface_normal = theRecoveredNormal ;
0265     const G4ThreeVector& direction      = OldMomentum ;
0266     const G4ThreeVector& polarization   = OldPolarization ;
0267
0268     G4double minus_cos_theta = direction*surface_normal ;
0269     G4double orientation = minus_cos_theta < 0. ? 1. : -1.  ;
0270     G4ThreeVector oriented_normal = orientation*surface_normal ;
0271
0272     G4double energy = thePhotonMomentum ;
0273     G4double wavelength = twopi*hbarc/energy ;
0274
0275     G4double energy_eV = energy/eV ;
0276     G4double wavelength_nm = wavelength/nm ;
0277
0278     const G4VTouchable* touch = aTrack.GetTouchable();
0279     int replicaNumber = U4Touchable::ReplicaNumber(touch);  // aka pmtid
0280     // TODO: add J API to access pmtcat+qe from replicaNumber
0281     int pmtcat = J::DEFAULT_CAT ;
0282     double _qe = 0.0 ;
0283
0284
0285     StackSpec<double,4> spec ;
0286     spec.ls[0].d = 0. ;
0287     spec.ls[1].d = parameter_accessor->get_thickness_nm( pmtcat, J::L1 );
0288     spec.ls[2].d = parameter_accessor->get_thickness_nm( pmtcat, J::L2 );
0289     spec.ls[3].d = 0. ;
0290
0291     spec.ls[0].nr = parameter_accessor->get_rindex( pmtcat, J::L0, J::RINDEX, energy_eV );
0292     spec.ls[0].ni = parameter_accessor->get_rindex( pmtcat, J::L0, J::KINDEX, energy_eV );
0293
0294     spec.ls[1].nr = parameter_accessor->get_rindex( pmtcat, J::L1, J::RINDEX, energy_eV );
0295     spec.ls[1].ni = parameter_accessor->get_rindex( pmtcat, J::L1, J::KINDEX, energy_eV );
0296
0297     spec.ls[2].nr = parameter_accessor->get_rindex( pmtcat, J::L2, J::RINDEX, energy_eV );
0298     spec.ls[2].ni = parameter_accessor->get_rindex( pmtcat, J::L2, J::KINDEX, energy_eV );
0299
0300     spec.ls[3].nr = parameter_accessor->get_rindex( pmtcat, J::L3, J::RINDEX, energy_eV );
0301     spec.ls[3].ni = parameter_accessor->get_rindex( pmtcat, J::L3, J::KINDEX, energy_eV );
0302
0303
0304     Stack<double,4> stack(      wavelength_nm, minus_cos_theta, spec );
0305
0306     // NB stack is flipped for minus_cos_theta > 0. so:
0307     //
0308     //    stack.ll[0] always incident side
0309     //    stack.ll[3] always transmission side
0310
0311     const double _ni = stack.ll[0].n.real() ;
0312     const double _si = stack.ll[0].st.real() ;
0313     const double _ci = stack.ll[0].ct.real() ;
0314
0315     const double _nt = stack.ll[3].n.real() ;
0316     const double _ct = stack.ll[3].ct.real() ;
0317     const double eta = _ni/_nt ;
0318
0319     Stack<double,4> stackNormal(wavelength_nm, -1.            , spec );
0320     // stackNormal is not flipped (as minus_cos_theta is fixed at -1.)
0321     //      presumably this is due to _qe definition
0322
0323
0324     // E_s2 : S-vs-P power fraction : signs make no difference as squared
0325     double E_s2 = _si > 0. ? (polarization*direction.cross(oriented_normal))/_si : 0. ;
0326     E_s2 *= E_s2;      // E_s2 matches E1_perp*E1_perp see sysrap/tests/stmm_vs_sboundary_test.cc
0327
0328     LOG(LEVEL)
0329         << " count " << count
0330         << " replicaNumber " << replicaNumber
0331         << " _si " << std::fixed << std::setw(10) << std::setprecision(5) << _si
0332         << " oriented_normal " << oriented_normal
0333         << " polarization*direction.cross(oriented_normal) " << polarization*direction.cross(oriented_normal)
0334         << " E_s2 " << std::fixed << std::setw(10) << std::setprecision(5) << E_s2
0335         ;
0336
0337     double fT_s = stack.art.T_s ;
0338     double fT_p = stack.art.T_p ;
0339     double fR_s = stack.art.R_s ;
0340     double fR_p = stack.art.R_p ;
0341     double one = 1.0 ;
0342     double T = fT_s*E_s2 + fT_p*(one-E_s2);  // matched with TransCoeff see sysrap/tests/stmm_vs_sboundary_test.cc
0343     double R = fR_s*E_s2 + fR_p*(one-E_s2);
0344     double A = one - (T+R);
0345
0346
0347    LOG(LEVEL)
0348         << " count " << count
0349         << " E_s2 " << std::fixed << std::setw(10) << std::setprecision(5) << E_s2
0350         << " fT_s " << std::fixed << std::setw(10) << std::setprecision(5) << fT_s
0351         << " 1-E_s2 " << std::fixed << std::setw(10) << std::setprecision(5) << (1.-E_s2)
0352         << " fT_p " << std::fixed << std::setw(10) << std::setprecision(5) << fT_p
0353         << " T " << std::fixed << std::setw(10) << std::setprecision(5) << T
0354         ;
0355
0356     LOG(LEVEL)
0357         << " count " << count
0358         << " E_s2 " << std::fixed << std::setw(10) << std::setprecision(5) << E_s2
0359         << " fR_s " << std::fixed << std::setw(10) << std::setprecision(5) << fR_s
0360         << " 1-E_s2 " << std::fixed << std::setw(10) << std::setprecision(5) << (1.-E_s2)
0361         << " fR_p " << std::fixed << std::setw(10) << std::setprecision(5) << fR_p
0362         << " R " << std::fixed << std::setw(10) << std::setprecision(5) << R
0363         << " A " << std::fixed << std::setw(10) << std::setprecision(5) << A
0364         ;
0365
0366
0367     double fT_n = stackNormal.art.T ;
0368     double fR_n = stackNormal.art.R ;
0369     double An = one - (fT_n+fR_n);
0370     double D = _qe/An;
0371
0372     LOG_IF(error, D > 1.)
0373          << " ERR: D > 1. : " << D
0374          << " _qe " << _qe
0375          << " An " << An
0376          ;
0377
0378     double u0 = G4UniformRand();
0379     double u1 = G4UniformRand();
0380
0381     /**
0382        0         A         A+R         1
0383        |---------+----------+----------|  u0
0384           D/A         R          T
0385           u1
0386     **/
0387
0388     char status =  u0 < A ? ( u1 < D ? 'D' : 'A' ) : ( u0 < A + R ? 'R' : 'T' ) ;
0389
0390
0391
0392 #ifdef DEBUG_TAG
0393     int u0_idx = U4UniformRand::Find(u0, SEvt::UU);
0394     int u1_idx = U4UniformRand::Find(u1, SEvt::UU);
0395
0396     LOG(LEVEL)
0397          << " u0 " << U4UniformRand::Desc(u0)
0398          << " u0_idx " << u0_idx
0399          << " A "   << std::setw(10) << std::fixed << std::setprecision(4) << A
0400          << " A+R " << std::setw(10) << std::fixed << std::setprecision(4) << (A+R)
0401          << " T "   << std::setw(10) << std::fixed << std::setprecision(4) << T
0402          << " status "
0403          << status
0404          << " DECISION "
0405          ;
0406     LOG(LEVEL)
0407          << " u1 " << U4UniformRand::Desc(u1)
0408          << " u1_idx " << u1_idx
0409          << " D " << std::setw(10) << std::fixed << std::setprecision(4) << D
0410          ;
0411 #endif
0412
0413
0414     NewMomentum = OldMomentum ;
0415     NewPolarization = OldPolarization ;
0416
0417     // The below reflect/refract is a copy of junoPMTOpticalModel
0418     //
0419     //    NewMomentum
0420     //        looks OK
0421     //
0422     //    NewPolarization
0423     //        looks to be guess(aka wrong)
0424     //        which does not follow along from the above S vs P fractions
0425     //
0426     //        polarization should depend on S and P fractions (think Brewsters angle)
0427     //        polarization vector does not reflect like the momentum vector does
0428     //
0429     // comparing with G4OpBoundaryProcess in sysrap/tests/stmm_vs_sboundary_test.cc
0430     //
0431
0432     if( status == 'R' )
0433     {
0434         theStatus = FresnelReflection ;
0435
0436         NewMomentum   -= 2.*(OldMomentum*oriented_normal)*oriented_normal ;
0437         NewPolarization -= 2.*(OldPolarization*oriented_normal)*oriented_normal ;
0438
0439        // THIS IS NOW G4OpBoundaryProcess CHANGES POLARIZATION FOR TOTAL INTERNAL REFLECTION,
0440        // NOT FRESNEL REFLECTION
0441
0442     }
0443     else if( status == 'T' )
0444     {
0445         theStatus = FresnelRefraction ;
0446         NewMomentum = eta*OldMomentum + (eta*_ci - _ct)*oriented_normal ;
0447
0448         // NewMomentum = NewMomentum.unit() ;
0449         // No need to normalize that is inherent. Derived in qsim.h:propagate_at_boundary.
0450         // This formula matches that used in qsim.h:propagate_at_boundary
0451         // The bracket is flipped compared with junoPMTOpticalModel.
0452         // TODO: check the oriented_normal sign convention, is it really flipped.
0453
0454         NewPolarization = (OldPolarization-(OldPolarization*NewMomentum)*NewMomentum).unit();
0455
0456         // JPOM expression:: pol = (pol-(pol*dir)*dir).unit();  // using old pol and new mom
0457         // The above just restates the JPOM expression, it does not match G4OpBoundaryProcess/qsim.h/sboundary.h
0458
0459
0460     }
0461     else if(status == 'A' || status == 'D')
0462     {
0463         theStatus = status == 'D' ? Detection : Absorption ;
0464
0465         aParticleChange.ProposeLocalEnergyDeposit(status == 'D' ? thePhotonMomentum : 0.0) ;
0466         aParticleChange.ProposeTrackStatus(fStopAndKill) ;
0467     }
0468
0469
0470
0471 #ifdef DEBUG_TAG
0472     SPhoton_Debug<'B'> dbg ;
0473
0474     dbg.pos = theGlobalPoint ;
0475     dbg.time = aTrack.GetLocalTime();  // just for debug, the only use of aTrack
0476
0477     dbg.mom = NewMomentum ;
0478     dbg.iindex = 0 ;
0479
0480     dbg.pol = NewPolarization ;
0481     dbg.wavelength = wavelength_nm ;
0482
0483     dbg.nrm = theRecoveredNormal ;
0484     dbg.spare = 0. ;
0485
0486     dbg.u0 = u0 ;
0487     dbg.x1 = 0. ;
0488     dbg.x2 = 0. ;
0489     dbg.u0_idx = 0 ;
0490
0491     LOG(LEVEL)
0492        << " time " << dbg.time
0493        << " dbg.Count " << dbg.Count()
0494        << " dbg.Name " << dbg.Name()
0495        << " status " << status
0496        << " CustomStatus::Name(status) " << CustomStatus::Name(status)
0497        ;
0498
0499     dbg.add();
0500 #endif
0501
0502     return status ;
0503 }
0504
0505