eic-opticks/u4/LocalCustomART.h

0001 #pragma once
0002 /**
0003 CustomART : Modifies G4OpBoundaryProcess with Custom calc of A/R/T coeffs
0004 ===========================================================================
0005
0006 * CustomART<JPMT> is instanciated in G4OpBoundaryProcess ctor
0007 * aims to provides customization with minimal code change to G4OpBoundaryProcess
0008 * aims to make few detector specific assumptions by moving such specifics into
0009   the template type that provides PMT parameter access
0010
0011 What detector/Geant4 geometry specific assumptions are made ?
0012 ---------------------------------------------------------------
0013
0014 1. ART customized for surfaces with names starting with '@' on local_z > 0.
0015 2. pmtid obtained from Geant4 volume ReplicaNumber
0016 3. templated PMT Parameter accessor type fulfils an API::
0017
0018     get_pmtcat(int pmtid) const ;
0019     get_pmtid_qe( int pmtid, double energy) const ;
0020     get_stackspec( std::array<double, 16>& a_spec, int pmtcat, double energy_eV ) const ;
0021     // HMM: thats assuming 4 layer stack
0022
0023
0024 Overview
0025 ----------
0026
0027 Trying to do less than CustomBoundary, instead just calculate::
0028
0029     theTransmittance
0030     theReflectivity
0031     theEfficiency
0032
0033 TODO: should also probably be setting::
0034
0035    type = dielectric_dielectric
0036    theFinish = polished
0037
0038 With everything else (deciding on ARTD, changing mom, pol, theStatus)
0039 done by the nearly "standard" G4OpBoundaryProcess.
0040
0041
0042 Is 2-layer (Pyrex,Vacuum) polarization direction calc applicable to 4-layer (Pyrex,ARC,PHC,Vacuum) situation ?
0043 -----------------------------------------------------------------------------------------------------------------
0044
0045 The Geant4 calculation of polarization direction for a simple
0046 boundary between two layers is based on continuity of E and B fields
0047 in S and P directions at the boundary, essentially Maxwells Eqn boundary conditions.
0048 Exactly the same thing yields Snells law::
0049
0050     n1 sin t1  = n2 sin t2
0051
0052 My thinking on this is that Snell's law with a bunch of layers would be::
0053
0054     n1 sin t1 = n2 sin t2 =  n3 sin t3 =  n4 sin t4
0055
0056 So the boundary conditions from 1->4 where n1 and n4 are real still gives::
0057
0058     n1 sin t1 = n4 sin t4
0059
0060 Even when n2,t2,n3,t3 are complex.
0061
0062 So by analogy that makes me think that the 2-layer polarization calculation
0063 between layers 1 and 4 (as done by G4OpBoundaryProcess)
0064 should still be valid even when there is a stack of extra layers
0065 inbetween layer 1 and 4.
0066
0067 Essentially the stack calculation changes A,R,T so it changes
0068 how much things happen : but it doesnt change what happens.
0069 So the two-layer polarization calculation from first and last layer
0070 should still be valid to the situation of the stack.
0071
0072 Do you agree with this argument ?
0073
0074 **/
0075
0076 #include "G4ThreeVector.hh"
0077
0078 #include "IPMTAccessor.h"
0079 #include "Layr.h"
0080 #include "U4Touchable.h"
0081
0082 struct CustomART
0083 {
0084     int    count ;
0085     double zlocal ;
0086     double lposcost ;
0087
0088     const IPMTAccessor* accessor ;
0089
0090     G4double& theTransmittance ;
0091     G4double& theReflectivity ;
0092     G4double& theEfficiency ;
0093
0094     const G4ThreeVector& theGlobalPoint ;
0095     const G4ThreeVector& OldMomentum ;
0096     const G4ThreeVector& OldPolarization ;
0097     const G4ThreeVector& theRecoveredNormal ;
0098     const G4double& thePhotonMomentum ;
0099
0100     CustomART(
0101         const IPMTAccessor* accessor,
0102         G4double& theTransmittance,
0103         G4double& theReflectivity,
0104         G4double& theEfficiency,
0105         const G4ThreeVector& theGlobalPoint,
0106         const G4ThreeVector& OldMomentum,
0107         const G4ThreeVector& OldPolarization,
0108         const G4ThreeVector& theRecoveredNormal,
0109         const G4double& thePhotonMomentum
0110     );
0111
0112     //char maybe_doIt(const char* OpticalSurfaceName, const G4Track& aTrack, const G4Step& aStep) ;
0113     double local_z( const G4Track& aTrack );
0114     void doIt(const G4Track& aTrack, const G4Step& aStep );
0115 };
0116
0117 inline CustomART::CustomART(
0118     const IPMTAccessor* accessor_,
0119           G4double& theTransmittance_,
0120           G4double& theReflectivity_,
0121           G4double& theEfficiency_,
0122     const G4ThreeVector& theGlobalPoint_,
0123     const G4ThreeVector& OldMomentum_,
0124     const G4ThreeVector& OldPolarization_,
0125     const G4ThreeVector& theRecoveredNormal_,
0126     const G4double&      thePhotonMomentum_
0127     )
0128     :
0129     count(0),
0130     zlocal(-1.),
0131     lposcost(-2.),
0132     accessor(accessor_),
0133     theTransmittance(theTransmittance_),
0134     theReflectivity(theReflectivity_),
0135     theEfficiency(theEfficiency_),
0136     theGlobalPoint(theGlobalPoint_),
0137     OldMomentum(OldMomentum_),
0138     OldPolarization(OldPolarization_),
0139     theRecoveredNormal(theRecoveredNormal_),
0140     thePhotonMomentum(thePhotonMomentum_)
0141 {
0142 }
0143
0144 /**
0145 CustomART::is_upper_z
0146 ------------------------
0147
0148 Q:What is lposcost for ?
0149
0150 A:Preparing for doing this on GPU, as lposcost is available there already but zlocal is not,
0151   so want to check the sign of lposcost is following that of zlocal. It looks
0152   like it should::
0153
0154     157 inline double Hep3Vector::cosTheta() const {
0155     158   double ptot = mag();
0156     159   return ptot == 0.0 ? 1.0 : dz/ptot;
0157     160 }
0158
0159
0160 **/
0161
0162 inline double CustomART::local_z( const G4Track& aTrack )
0163 {
0164     const G4AffineTransform& transform = aTrack.GetTouchable()->GetHistory()->GetTopTransform();
0165     G4ThreeVector localPoint = transform.TransformPoint(theGlobalPoint);
0166     zlocal = localPoint.z() ;
0167     lposcost = localPoint.cosTheta() ;
0168     return zlocal  ;
0169 }
0170
0171
0172 /**
0173 CustomART::maybe_doIt
0174 ------------------------
0175
0176 As need to also handle traditional POM it avoids duplication
0177 for the maybe_doIt branching to happen at a higher level inside CustomG4OpBoundaryProcess
0178
0179 inline char CustomART::maybe_doIt(const char* OpticalSurfaceName, const G4Track& aTrack, const G4Step& aStep )
0180 {
0181     if( OpticalSurfaceName == nullptr || OpticalSurfaceName[0] != '@') return 'N' ;
0182
0183     return doIt(aTrack, aStep) ;
0184 }
0185
0186 **/
0187
0188
0189
0190 /**
0191 CustomART<J>::doIt
0192 --------------------
0193
0194 NB stack is flipped for minus_cos_theta > 0. so:
0195
0196 * stack.ll[0] always incident side
0197 * stack.ll[3] always transmission side
0198
0199 **/
0200
0201 inline void CustomART::doIt(const G4Track& aTrack, const G4Step& aStep )
0202 {
0203     G4double minus_cos_theta = OldMomentum*theRecoveredNormal ;
0204
0205     G4double energy = thePhotonMomentum ;
0206     G4double wavelength = twopi*hbarc/energy ;
0207     G4double energy_eV = energy/eV ;
0208     G4double wavelength_nm = wavelength/nm ;
0209
0210     const G4VTouchable* touch = aTrack.GetTouchable();
0211     int pmtid = U4Touchable::ReplicaNumber(touch);
0212
0213     int pmtcat = accessor->get_pmtcat( pmtid ) ;
0214     double _qe = minus_cos_theta > 0. ? 0.0 : accessor->get_pmtid_qe( pmtid, energy ) ;
0215     // following the old junoPMTOpticalModel with "backwards" _qe always zero
0216
0217     std::array<double,16> a_spec ;
0218     accessor->get_stackspec(a_spec, pmtcat, energy_eV );
0219     StackSpec<double,4> spec ;
0220     spec.import( a_spec );
0221
0222     bool dump = false ;
0223
0224     if(dump) std::cerr << "CustomART::doIt" << " spec " << std::endl << spec << std::endl ;
0225
0226     // SUSPECT: there will be significant repetition of get_stackspec
0227     // calls with same pmtcat and energy_eV : so could try caching ?
0228
0229     Stack<double,4> stack(wavelength_nm, minus_cos_theta, spec );
0230
0231     const double _si = stack.ll[0].st.real() ;
0232     const double _si2 = sqrtf( 1. - minus_cos_theta*minus_cos_theta );
0233
0234     double E_s2 = _si > 0. ? (OldPolarization*OldMomentum.cross(theRecoveredNormal))/_si : 0. ;
0235     E_s2 *= E_s2;
0236
0237     // E_s2 : S-vs-P power fraction : signs make no difference as squared
0238     // E_s2 matches E1_perp*E1_perp see sysrap/tests/stmm_vs_sboundary_test.cc
0239
0240     double one = 1.0 ;
0241     double S = E_s2 ;
0242     double P = one - S ;
0243
0244     if(dump) std::cerr
0245         << "CustomART::doIt"
0246         << " count " << count
0247         << " pmtid " << pmtid
0248         << " _si " << std::fixed << std::setw(10) << std::setprecision(5) << _si
0249         << " _si2 " << std::fixed << std::setw(10) << std::setprecision(5) << _si2
0250         << " theRecoveredNormal " << theRecoveredNormal
0251         << " OldPolarization*OldMomentum.cross(theRecoveredNormal) " << OldPolarization*OldMomentum.cross(theRecoveredNormal)
0252         << " E_s2 " << std::fixed << std::setw(10) << std::setprecision(5) << E_s2
0253         << std::endl
0254         ;
0255
0256     double T = S*stack.art.T_s + P*stack.art.T_p ;  // matched with TransCoeff see sysrap/tests/stmm_vs_sboundary_test.cc
0257     double R = S*stack.art.R_s + P*stack.art.R_p ;
0258     double A = one - (T+R);
0259
0260     theTransmittance = T ;
0261     theReflectivity  = R ;
0262
0263     if(dump) std::cerr
0264         << "CustomART::doIt"
0265         << " count " << count
0266         << " S " << std::fixed << std::setw(10) << std::setprecision(5) << S
0267         << " P " << std::fixed << std::setw(10) << std::setprecision(5) << P
0268         << " T " << std::fixed << std::setw(10) << std::setprecision(5) << T
0269         << " R " << std::fixed << std::setw(10) << std::setprecision(5) << R
0270         << " A " << std::fixed << std::setw(10) << std::setprecision(5) << A
0271         << std::endl
0272         ;
0273
0274     // stackNormal is not flipped (as minus_cos_theta is fixed at -1.) presumably this is due to _qe definition
0275     Stack<double,4> stackNormal(wavelength_nm, -1. , spec );
0276
0277     double An = one - (stackNormal.art.T + stackNormal.art.R) ;
0278     double D = _qe/An;   // LOOKS STRANGE TO DIVIDE BY An  : HMM MAYBE NEED TO DIVIDE BY A TOO ?
0279
0280     theEfficiency = D ;
0281
0282     bool expect = D <= 1. ;
0283     if(!expect) std::cerr
0284         << "CustomART::doIt"
0285         << " FATAL "
0286         << " ERR: D > 1. : " << D
0287         << " _qe " << _qe
0288         << " An " << An
0289         << std::endl
0290         ;
0291
0292     assert( expect );
0293 }
0294
0295