eic-opticks/u4/U4Navigator.h

0001 #pragma once
0002 /**
0003 U4Navigator.h
0004 ===============
0005
0006
0007 **/
0008
0009 #include <cassert>
0010 #include <iostream>
0011
0012 #include "scuda.h"
0013 #include "squad.h"
0014
0015 #include "G4Navigator.hh"
0016 #include "G4TransportationManager.hh"
0017 #include "U4Tree.h"
0018
0019
0020 struct U4Navigator_Stats
0021 {
0022     int num_call = 0 ;
0023     int num_nullpv0 = 0 ;
0024     int num_withpv0 = 0 ;
0025
0026     int num_isect = 0 ;
0027     int num_nullpv1 = 0 ;
0028     int num_withpv1 = 0 ;
0029
0030     std::string desc() const ;
0031 };
0032
0033 inline std::string U4Navigator_Stats::desc() const
0034 {
0035     double frac_nullpv0 = double(num_nullpv0)/double(num_call) ;
0036     double frac_nullpv1 = double(num_nullpv1)/double(num_isect) ;
0037     std::stringstream ss ;
0038     ss
0039         << " num_call " << num_call
0040         << " num_isect " << num_isect
0041         << " num_nullpv0 " << num_nullpv0
0042         << " num_withpv0 " << num_withpv0
0043         << " num_nullpv1 " << num_nullpv1
0044         << " num_withpv1 " << num_withpv1
0045         << " frac_nullpv0 " << std::setw(10) << std::setprecision(3) << std::fixed << frac_nullpv0
0046         << " frac_nullpv1 " << std::setw(10) << std::setprecision(3) << std::fixed << frac_nullpv1
0047         ;
0048
0049     std::string str = ss.str() ;
0050     return str ;
0051 }
0052
0053
0054 struct U4Navigator_Intersect
0055 {
0056     G4ThreeVector ipos = {} ;
0057     G4VPhysicalVolume* pv0 = nullptr ;
0058     G4VPhysicalVolume* pv1 = nullptr ;
0059     G4LogicalVolume* lv1 = nullptr ;
0060     G4VSolid* so1 = nullptr ;
0061     G4double t = 0. ;
0062     G4bool valid = false ;
0063     int nidx = -10 ;
0064     int prim = -10 ;
0065     const char* prn = nullptr ;
0066
0067     void zero();
0068     std::string desc() const ;
0069 };
0070
0071
0072 inline void U4Navigator_Intersect::zero()
0073 {
0074     ipos = {} ;
0075     pv0 = nullptr ;
0076     pv1 = nullptr ;
0077     lv1 = nullptr ;
0078     so1 = nullptr ;
0079     t = 0 ;
0080     valid = false ;
0081     nidx = -10 ;
0082     prim = -10 ;
0083     prn = nullptr ;
0084 }
0085
0086 inline std::string U4Navigator_Intersect::desc() const
0087 {
0088     std::stringstream ss ;
0089     ss
0090         << " t " << std::setw(10) << std::fixed << std::setprecision(3) << t
0091         << " valid "  << ( valid ? "YES" : "NO " )
0092         << " nidx "  << std::setw(8) << nidx
0093         << " prim "  << std::setw(6) << prim
0094         << " prn " << ( prn ? prn : "-" )
0095         << " pv0 " << ( pv0 ? pv0->GetName() : "-" )
0096         << " pv1 " << ( pv1 ? pv1->GetName() : "-" )
0097         << " so1 " << ( so1 ? so1->GetName() : "-" )
0098         ;
0099
0100     std::string str = ss.str() ;
0101     return str ;
0102 }
0103
0104
0105
0106 struct U4Navigator
0107 {
0108     static G4Navigator* GetNav();
0109     static double Distance(const G4ThreeVector& ori, const G4ThreeVector& dir );
0110
0111
0112     U4Navigator(const U4Tree* tree);
0113
0114     const U4Tree* tree ;
0115     G4Navigator* nav ;
0116     U4Navigator_Stats stats = {} ;
0117     U4Navigator_Intersect isect = {} ;
0118
0119     void simtrace(quad4& p );
0120     void getIntersect( const G4ThreeVector& ori, const G4ThreeVector& dir );
0121 };
0122
0123 inline G4Navigator* U4Navigator::GetNav()
0124 {
0125     G4Navigator* _nav = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
0126     return _nav ;
0127 }
0128
0129
0130 inline double U4Navigator::Distance( const G4ThreeVector& ori, const G4ThreeVector& dir )
0131 {
0132     U4Navigator nav(nullptr);
0133     nav.getIntersect( ori, dir);
0134     return nav.isect.t ;
0135 }
0136
0137
0138
0139 inline U4Navigator::U4Navigator(const U4Tree* _tree)
0140     :
0141     tree(_tree),
0142     nav(nullptr)  // GetNav when needed
0143 {
0144 }
0145
0146
0147 /**
0148 U4Navigator::getIntersect
0149 ---------------------------
0150
0151 Canonical caller U4Navigator::simtrace
0152
0153 **/
0154
0155
0156 inline void U4Navigator::getIntersect( const G4ThreeVector& ori, const G4ThreeVector& dir )
0157 {
0158     isect.zero();
0159
0160     if( nav == nullptr ) nav = GetNav();
0161
0162     G4VPhysicalVolume* pv0 = nullptr ;
0163     {
0164         const G4bool pRelativeSearch=false ;
0165         const G4bool ignoreDirection=false ;
0166         pv0 = nav->LocateGlobalPointAndSetup(ori, &dir, pRelativeSearch, ignoreDirection );
0167     }
0168     stats.num_call++ ;
0169     if(pv0 == nullptr) stats.num_nullpv0++ ;
0170     if(pv0 != nullptr) stats.num_withpv0++ ;
0171
0172
0173     G4double t ;
0174     {
0175         const G4double pCurrentProposedStepLength = kInfinity ;
0176         G4double pNewSafety ;
0177         t = nav->ComputeStep(ori, dir, pCurrentProposedStepLength, pNewSafety);
0178     }
0179
0180     G4bool valid = t > 0. && t != kInfinity ;
0181
0182     isect.pv0 = pv0 ;
0183     isect.t = t ;
0184     isect.valid = valid ;
0185
0186     assert(tree);
0187
0188     if( valid )
0189     {
0190         G4ThreeVector ipos = ori + t*dir ;
0191         G4VPhysicalVolume* pv1 = nav->LocateGlobalPointAndSetup(ipos, nullptr, false, true );
0192         G4LogicalVolume* lv1 = pv1 ? pv1->GetLogicalVolume() : nullptr ;
0193         G4VSolid* so1 = lv1 ? lv1->GetSolid() : nullptr ;
0194
0195         isect.ipos = ipos ;
0196         isect.pv1 = pv1 ;
0197         isect.lv1 = lv1 ;
0198         isect.so1 = so1 ;
0199
0200
0201         // try to identify the intersected volume
0202         int nidx        = pv1  == nullptr ? -1      : tree->get_nidx( pv1 )           ; // look for pv1 in stree::pvs vector
0203         int prim        = nidx == -1      ? -1      : tree->get_prim_for_nidx( nidx ) ;
0204         const char* prn = prim == -1      ? nullptr : tree->get_prname( prim )        ;
0205
0206         isect.nidx = nidx ;
0207         isect.prim = prim ;
0208         isect.prn = prn ;
0209
0210         stats.num_isect += 1;
0211         if(isect.pv1 == nullptr) stats.num_nullpv1++ ;
0212         if(isect.pv1 != nullptr) stats.num_withpv1++ ;
0213     }
0214
0215 }
0216
0217 /**
0218 U4Navigator::simtrace
0219 -----------------------
0220
0221 Canonical callstack::
0222
0223     U4Simtrace::EndOfRunAction > U4Simtrace::Scan
0224
0225 Assuming standard simtrace layout (see sevent::add_simtrace and below)
0226 this uses U4Navigator::Distance to populate the intersect position.
0227 HMM: how to get the surface normal at the intersect position ?
0228
0229 * will need to get the G4VSolid ?
0230
0231 * q0.f.xyz normal at intersect (not implemented) q0.f.z distance to intersect
0232 * q1.f.xyz intersect position    q1.f.w tmin
0233 * q2.f.xyz trace origin
0234 * q3.f.xyz trace direction
0235
0236 The simtrace layout is documented in sevent::add_simtrace
0237
0238 **/
0239
0240 inline void U4Navigator::simtrace( quad4& p )
0241 {
0242     G4ThreeVector ori(p.q2.f.x, p.q2.f.y, p.q2.f.z);
0243     G4ThreeVector dir(p.q3.f.x, p.q3.f.y, p.q3.f.z);
0244
0245     getIntersect(ori, dir);
0246
0247     p.q0.f.x = 0.f ; // not yet implemented : surface normal at intersect
0248     p.q0.f.y = 0.f ;
0249     p.q0.f.z = 0.f ;
0250     p.q0.f.w = isect.t ;
0251
0252     p.q1.f.x = float(isect.ipos.x()) ;
0253     p.q1.f.y = float(isect.ipos.y()) ;
0254     p.q1.f.z = float(isect.ipos.z()) ;
0255     p.q1.f.w = 0.f  ; // tmin
0256
0257     unsigned globalPrimIdx = isect.prim == -1 ? 0xffffu : ( isect.prim & 0xffffu ) ;
0258     // Formerly isect.nidx but as that is not equivalent globalPrimIdx
0259     // so have added U4Tree/stree functionality to get the globalPrimIdx
0260     // rather than CSGFoundry which is the natural place to get globalPrimIdx from.
0261
0262     unsigned boundary = 0u & 0xffffu ;
0263     unsigned globalPrimIdx_boundary = ( globalPrimIdx << 16 ) | boundary  ;
0264     unsigned identity = 0u ;
0265
0266     p.q2.u.w = globalPrimIdx_boundary ;
0267     p.q3.u.w = identity ;
0268 }
0269
0270