eic-opticks/sysrap/ssolid.h

0001 #pragma once
0002 /**
0003 ssolid.h
0004 ===========
0005
0006 **/
0007
0008 #include <cassert>
0009 #include <iostream>
0010 #include <iomanip>
0011
0012 #include <glm/glm.hpp>
0013 #include "scuda.h"
0014 #include "squad.h"
0015 #include "sgeomdefs.h"
0016
0017 #include "G4VSolid.hh"
0018 #include "G4ThreeVector.hh"
0019
0020 struct ssolid
0021 {
0022     static void GetCenterExtent( float4& ce,             const G4VSolid* solid );
0023     static void GetCenterExtent( glm::tvec4<double>& ce, const G4VSolid* solid );
0024     static G4double Distance_(const G4VSolid* solid, const G4ThreeVector& pos, const G4ThreeVector& dir, EInside& in,
0025              G4ThreeVector* isect=nullptr
0026          );
0027     static G4double Distance(const G4VSolid* solid, const G4ThreeVector& pos, const G4ThreeVector& dir, bool dump );
0028     static void Simtrace( quad4& p, const G4VSolid* solid, bool dump=false);
0029 };
0030
0031 inline void ssolid::GetCenterExtent( glm::tvec4<double>& ce, const G4VSolid* solid ) // static
0032 {
0033     G4ThreeVector pMin ;
0034     G4ThreeVector pMax ;
0035     solid->BoundingLimits(pMin, pMax);
0036     G4ThreeVector center = ( pMin + pMax )/2. ;
0037     G4ThreeVector fulldiag = pMax - pMin ;
0038     G4ThreeVector halfdiag = fulldiag/2.  ;
0039     G4double extent = std::max( std::max( halfdiag.x(), halfdiag.y() ), halfdiag.z() ) ;
0040
0041     ce.x = center.x() ;
0042     ce.y = center.y() ;
0043     ce.z = center.z() ;
0044     ce.w = extent ;
0045 }
0046
0047 inline void ssolid::GetCenterExtent( float4& ce, const G4VSolid* solid ) // static
0048 {
0049     glm::tvec4<double> ce_ ;
0050     GetCenterExtent(ce_, solid );
0051     ce.x = ce_.x ;
0052     ce.y = ce_.y ;
0053     ce.z = ce_.z ;
0054     ce.w = ce_.w ;
0055 }
0056
0057
0058 /**
0059 ssolid::Distance_
0060 ------------------
0061
0062 See u4/tests/G4Orb_Test.cc
0063
0064 **/
0065
0066 inline G4double ssolid::Distance_(const G4VSolid* solid, const G4ThreeVector& pos, const G4ThreeVector& dir, EInside& in,
0067         G4ThreeVector* isect
0068      ) // static
0069 {
0070     in =  solid->Inside(pos) ;
0071     G4double t = kInfinity ;
0072     switch( in )
0073     {
0074         case kInside:  t = solid->DistanceToOut( pos, dir ) ; break ;
0075         case kSurface: t = solid->DistanceToOut( pos, dir ) ; break ;
0076         case kOutside: t = solid->DistanceToIn(  pos, dir ) ; break ;
0077         default:  assert(0) ;
0078     }
0079     if( t != kInfinity && isect != nullptr ) *isect = pos+t*dir ;
0080     return t ;
0081 }
0082
0083
0084
0085 inline G4double ssolid::Distance(const G4VSolid* solid, const G4ThreeVector& pos, const G4ThreeVector& dir, bool dump ) // static
0086 {
0087     EInside in ;
0088     G4double t = Distance_( solid, pos, dir, in );
0089
0090     if(dump && t != kInfinity)
0091     {
0092         std::cout
0093             << " pos "
0094             << "("
0095             << std::fixed << std::setw(10) << std::setprecision(3) << pos.x() << " "
0096             << std::fixed << std::setw(10) << std::setprecision(3) << pos.y() << " "
0097             << std::fixed << std::setw(10) << std::setprecision(3) << pos.z()
0098             << ")"
0099             << " dir "
0100             << "("
0101             << std::fixed << std::setw(10) << std::setprecision(3) << dir.x() << " "
0102             << std::fixed << std::setw(10) << std::setprecision(3) << dir.y() << " "
0103             << std::fixed << std::setw(10) << std::setprecision(3) << dir.z()
0104             << ")"
0105             << " in " << sgeomdefs::EInside_(in )
0106             ;
0107
0108        if( t == kInfinity)
0109        {
0110             std::cout
0111                 << " t " << std::setw(10) << "kInfinity"
0112                 << std::endl
0113                 ;
0114        }
0115        else
0116        {
0117            G4ThreeVector ipos = pos + dir*t ;
0118            std::cout
0119                 << " t " << std::fixed << std::setw(10) << std::setprecision(3) << t
0120                 << " ipos "
0121                 << "("
0122                 << std::fixed << std::setw(10) << std::setprecision(3) << ipos.x() << " "
0123                 << std::fixed << std::setw(10) << std::setprecision(3) << ipos.y() << " "
0124                 << std::fixed << std::setw(10) << std::setprecision(3) << ipos.z()
0125                 << ")"
0126                 << std::endl
0127                 ;
0128        }
0129     }
0130     return t ;
0131 }
0132
0133
0134 /**
0135 ssolid::Simtrace
0136 -------------------
0137
0138 Updates quad4& p in simtrace layout with intersect position onto the solid.
0139
0140 p.q0.f.xyz,w
0141     surface normal at intersect (UNCHECKED) and intersect distance *t*
0142
0143 p.q1.f.xyz, w
0144     intersect position and "tmin" param
0145
0146 p.q2.f.xyz
0147     input position (assumed to be in local frame of solid)
0148
0149 p.q3.f.xyz
0150     input direction (assumed to be in local frame of solid)
0151
0152 **/
0153
0154 inline void ssolid::Simtrace(quad4& p, const G4VSolid* solid, bool dump) // static
0155 {
0156     G4ThreeVector ori(p.q2.f.x, p.q2.f.y, p.q2.f.z);
0157     G4ThreeVector dir(p.q3.f.x, p.q3.f.y, p.q3.f.z);
0158
0159     G4double t = Distance( solid, ori, dir, dump );
0160     //std::cout << "ssolid::Simtrace " << t << std::endl ;
0161
0162     if( t == kInfinity ) return ;   // hmm: perhaps set ipos to ori for MISS ? Currently gets left at origin
0163
0164     G4ThreeVector ipos = ori + dir*t ;
0165     float tmin = 0.f ;
0166
0167     G4ThreeVector inrm = solid->SurfaceNormal(ipos) ; // UNCHECKED
0168     p.q0.f.x = float(inrm.x()) ;
0169     p.q0.f.y = float(inrm.y()) ;
0170     p.q0.f.z = float(inrm.z()) ;
0171     p.q0.f.w = t ;
0172
0173     p.q1.f.x = float(ipos.x()) ;
0174     p.q1.f.y = float(ipos.y()) ;
0175     p.q1.f.z = float(ipos.z()) ;
0176     p.q1.f.w = tmin  ;
0177 }
0178
0179