eic-opticks/sysrap/SVolume.h

0001 #pragma once
0002 /**
0003 SVolume.h
0004 =============
0005
0006 Developed within j/PMTSim but as generally useful moved to sysrap header-only.
0007
0008 An early user of this was x4/tests/X4IntersectVolumeTest.cc as used by x4/xxv.sh
0009
0010 HUH: thats wierd, lots of passing around vectors of solids just to get their names
0011
0012 **/
0013
0014 #include <vector>
0015 #include <string>
0016 #include <array>
0017
0018 struct NP ;
0019 class G4VSolid ;
0020 class G4LogicalVolume ;
0021 class G4VPhysicalVolume ;
0022
0023 struct SVolume
0024 {
0025     static void SaveTransforms( std::vector<double>* tr, std::vector<G4VSolid*>* names, const char* path );
0026     static void SaveTransforms( std::vector<double>* tr, std::vector<G4VSolid*>* names, const char* fold, const char* name );
0027
0028     static void GetSolidNames(std::vector<std::string>& names,  std::vector<G4VSolid*>* solids );
0029
0030     static NP* MakeArray( std::vector<double>* tr, std::vector<G4VSolid*>* solids );
0031     static void DumpTransforms( std::vector<double>* tr, std::vector<G4VSolid*>* names, const char* msg );
0032
0033     static void Traverse(const G4VPhysicalVolume* const pv, std::vector<double>* tr, std::vector<G4VSolid*>* names );
0034     static void Traverse_r(const G4VPhysicalVolume* const pv, int depth, std::vector<double>* tr, std::vector<G4VSolid*>* names);
0035     static bool IsIdentityRotation(const std::array<double, 16>& a, double epsilon ) ;
0036     static void GetObjectTransform(std::array<double, 16>& a, const G4VPhysicalVolume* const pv);
0037     static void DumpSolids();
0038 };
0039
0040
0041 #include <cassert>
0042 #include <cstdlib>
0043 #include <cstring>
0044 #include <iostream>
0045 #include <iomanip>
0046 #include <array>
0047
0048 #include "NP.hh"
0049
0050 #include "G4LogicalVolume.hh"
0051 #include "G4PVPlacement.hh"
0052 #include "G4SolidStore.hh"
0053 #include "G4DisplacedSolid.hh"
0054 #include "G4Material.hh"
0055
0056
0057 inline void SVolume::SaveTransforms( std::vector<double>* tr, std::vector<G4VSolid*>* solids, const char* fold, const char* name ) // static
0058 {
0059     NP* a = MakeArray(tr, solids);
0060     a->save(fold, name);
0061 }
0062
0063 inline void SVolume::SaveTransforms( std::vector<double>* tr, std::vector<G4VSolid*>* solids, const char* path ) // static
0064 {
0065     NP* a = MakeArray(tr, solids);
0066     a->save(path);
0067 }
0068
0069 inline void SVolume::GetSolidNames(std::vector<std::string>& names,  std::vector<G4VSolid*>* solids ) // static
0070 {
0071     unsigned num = solids->size() ;
0072     for(unsigned i=0 ; i < num ; i++)
0073     {
0074         G4VSolid* solid = (*solids)[i] ;
0075         G4String name = solid->GetName() ;
0076         names.push_back(name);
0077     }
0078 }
0079
0080 /**
0081 SVolume::MakeArray
0082 ---------------------
0083
0084 Creates array of transforms with the names of the solids as metadata.
0085
0086
0087
0088 **/
0089 inline NP* SVolume::MakeArray( std::vector<double>* tr, std::vector<G4VSolid*>* solids ) // static
0090 {
0091     unsigned num_solid = solids->size() ;
0092     unsigned num_value = tr->size();
0093     assert( num_value % 16 == 0 );
0094     unsigned num_transform = num_value/num_solid/16 ;
0095
0096     bool expect = num_value == num_solid*num_transform*16 ;
0097
0098     if(!expect) std::cerr
0099         << "SVolume::MakeArray"
0100         << " num_solid " << num_solid
0101         << " num_value " << num_value
0102         << " num_transform " << num_transform
0103         << " num_solid*num_transform*16 " << num_solid*num_transform*16
0104         << std::endl
0105         ;
0106     assert(expect);
0107
0108     std::vector<std::string> names ;
0109     GetSolidNames(names, solids) ;
0110
0111     NP* a = NP::Make<double>( num_solid, num_transform, 4, 4 );   // NB changed shape to allow multiple transforms per solid
0112     a->read( tr->data() );
0113     a->set_names(names) ;  // NB changed, to set_names formerly used a->meta string
0114
0115     return a ;
0116 }
0117
0118 inline void SVolume::DumpTransforms( std::vector<double>* tr, std::vector<G4VSolid*>* solids, const char* msg ) // static
0119 {
0120     std::cout << msg << std::endl ;
0121     unsigned num = solids->size() ;
0122     assert( tr->size() == num*16 );
0123     double* trd = tr->data() ;
0124
0125     for(unsigned i=0 ; i < num ; i++)
0126     {
0127         G4VSolid* solid = (*solids)[i] ;
0128         G4String soname = solid->GetName() ;
0129         std::cout
0130             << std::setw(5) << i
0131             << " : "
0132             << std::setw(25) << soname
0133             ;
0134
0135         std::cout << " tr ( " ;
0136         for(int j=0 ; j < 16 ; j++ )
0137             std::cout << std::fixed << std::setw(7) << std::setprecision(2) << trd[i*16+j] << " " ;
0138         std::cout << ")" << std::endl ;
0139     }
0140 }
0141
0142 /**
0143 SVolume::Traverse
0144 -------------------
0145
0146 This is used for example by PMTFastSim::GetPV
0147
0148 **/
0149
0150 inline void SVolume::Traverse(const G4VPhysicalVolume* const pv, std::vector<double>* tr , std::vector<G4VSolid*>* solids ) // static
0151 {
0152     Traverse_r( pv, 0, tr, solids);
0153 }
0154
0155 /**
0156 SVolume::Traverse_r
0157 ---------------------
0158
0159 Recursively traverses the volume hierarchy collecting single level of transforms and solids.
0160 NB Assumes no multi-level transforms, so in this form its only appropriate for very simple volume trees.
0161
0162 **/
0163
0164 inline void SVolume::Traverse_r(const G4VPhysicalVolume* const pv, int depth, std::vector<double>* tr, std::vector<G4VSolid*>* solids ) // static
0165 {
0166     const G4LogicalVolume* lv = pv->GetLogicalVolume() ;
0167     G4VSolid* so = lv->GetSolid();
0168     const G4Material* const mt = lv->GetMaterial() ;
0169     G4String soname = so->GetName();
0170
0171     std::array<double, 16> a ;
0172     GetObjectTransform(a, pv);
0173     if(tr) for(unsigned i=0 ; i < 16 ; i++) tr->push_back( a[i] );
0174     if(solids) solids->push_back(so);
0175
0176     std::cout
0177         << "SVolume::Traverse_r"
0178         << " depth " << std::setw(2) << depth
0179         << " pv " << std::setw(22) << pv->GetName()
0180         << " lv " << std::setw(22) << lv->GetName()
0181         << " so " << std::setw(22) << so->GetName()
0182         << " mt " << std::setw(15) << mt->GetName()
0183         ;
0184
0185     bool identity_rot = IsIdentityRotation( a, 1e-6 );
0186     std::cout << " tr ( " ;
0187     unsigned i0 = identity_rot ? 4*3 : 0 ;
0188     for(unsigned i=i0 ; i < 16 ; i++) std::cout << std::fixed << std::setw(7) << std::setprecision(2) << a[i] << " " ;
0189     std::cout << ") " << std::endl ;
0190
0191     for (size_t i=0 ; i < size_t(lv->GetNoDaughters()) ;i++ )
0192     {
0193         const G4VPhysicalVolume* const daughter_pv = lv->GetDaughter(i);
0194         Traverse_r( daughter_pv, depth+1, tr, solids );
0195     }
0196 }
0197
0198 inline bool SVolume::IsIdentityRotation(const std::array<double, 16>& a, double epsilon )
0199 {
0200     unsigned not_identity=0 ;
0201     for(unsigned i=0 ; i < 3 ; i++)
0202     for(unsigned j=0 ; j < 3 ; j++)
0203     {
0204         unsigned idx = i*4 + j ;
0205         double expect = i == j ? 1. : 0. ;
0206         if( std::abs( a[idx] - expect) > epsilon ) not_identity += 1 ;
0207     }
0208     return not_identity == 0 ;
0209 }
0210
0211
0212 inline void SVolume::GetObjectTransform(std::array<double, 16>& a, const G4VPhysicalVolume* const pv)
0213 {
0214    // preferred for interop with glm/Opticks : obj relative to mother
0215     G4RotationMatrix rot = pv->GetObjectRotationValue() ;
0216     G4ThreeVector    tla = pv->GetObjectTranslation() ;
0217     G4Transform3D    t(rot,tla);
0218
0219     a[ 0] = t.xx() ; a[ 1] = t.yx() ; a[ 2] = t.zx() ; a[ 3] = 0.f    ;
0220     a[ 4] = t.xy() ; a[ 5] = t.yy() ; a[ 6] = t.zy() ; a[ 7] = 0.f    ;
0221     a[ 8] = t.xz() ; a[ 9] = t.yz() ; a[10] = t.zz() ; a[11] = 0.f    ;
0222     a[12] = t.dx() ; a[13] = t.dy() ; a[14] = t.dz() ; a[15] = 1.f    ;
0223 }
0224
0225
0226 inline void SVolume::DumpSolids() // static
0227 {
0228     G4SolidStore* store = G4SolidStore::GetInstance() ;
0229     std::cout << "SVolume::DumpSolids G4SolidStore.size " << store->size() << std::endl ;
0230     for( unsigned i=0 ; i < store->size() ; i++)
0231     {
0232         G4VSolid* solid = (*store)[i] ;
0233
0234         G4DisplacedSolid* disp = dynamic_cast<G4DisplacedSolid*>(solid) ;
0235         G4VSolid* moved = disp ? disp->GetConstituentMovedSolid() : nullptr ;
0236
0237         std::cout
0238             << " i " << std::setw(5) << i
0239             << " solid.name "
0240             << std::setw(30) << solid->GetName()
0241             << " moved.name "
0242             << std::setw(30) << ( moved ? moved->GetName() : "-" )
0243             << std::endl
0244             ;
0245     }
0246 }
0247