Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-10 07:50:36 

0001 #include <sstream>
0002 #include <csignal>
0003 #include <array>
0004 
0005 #include "G4Material.hh"
0006 #include "G4ThreeVector.hh"
0007 #include "G4LogicalVolume.hh"
0008 #include "G4PVPlacement.hh"
0009 #include "G4Box.hh"
0010 #include "G4Orb.hh"
0011 
0012 #include "spath.h"
0013 #include "sstr.h"
0014 #include "ssys.h"
0015 #include "s_bb.h"
0016 
0017 #include "SPlace.h"
0018 
0019 #include "SLOG.hh"
0020 
0021 #include "U4.hh"
0022 #include "U4Material.hh"
0023 #include "U4Surface.h"
0024 #include "U4SolidMaker.hh"
0025 #include "U4VolumeMaker.hh"
0026 #include "U4RotationMatrix.h"
0027 #include "U4Volume.h"
0028 #include "U4GDML.h"
0029 #include "U4ThreeVector.h"
0030 
0031 #ifdef WITH_PMTSIM
0032 #include "PMTSim.hh"
0033 #endif
0034 
0035 
0036 
0037 const plog::Severity U4VolumeMaker::LEVEL = SLOG::EnvLevel("U4VolumeMaker", "DEBUG");
0038 const char* U4VolumeMaker::GEOM = ssys::getenvvar("GEOM", "BoxOfScintillator");
0039 const char* U4VolumeMaker::METH = nullptr ;
0040 
0041 std::string U4VolumeMaker::Desc() // static
0042 {
0043     std::stringstream ss ;
0044     ss << "U4VolumeMaker::Desc" ;
0045     ss << " GEOM " << ( GEOM ? GEOM : "-" ) ;
0046     ss << " METH " << ( METH ? METH : "-" ) ;
0047 #ifdef WITH_PMTSIM
0048     ss << " WITH_PMTSIM " ;
0049 #else
0050     ss << " not-WITH_PMTSIM " ;
0051 #endif
0052     std::string s = ss.str();
0053     return s ;
0054 }
0055 
0056 std::string U4VolumeMaker::Desc( const G4ThreeVector& tla )
0057 {
0058     std::stringstream ss ;
0059     ss << " tla ("
0060        <<  " " << std::fixed << std::setw(10) << std::setprecision(3) << tla.x()
0061        <<  " " << std::fixed << std::setw(10) << std::setprecision(3) << tla.y()
0062        <<  " " << std::fixed << std::setw(10) << std::setprecision(3) << tla.z()
0063        <<  ") "
0064        ;
0065 
0066     std::string s = ss.str();
0067     return s ;
0068 }
0069 
0070 
0071 
0072 
0073 
0074 /**
0075 U4VolumeMaker::PV
0076 ---------------------
0077 
0078 Invokes several PV getter methods the first to provide a PV wins.
0079 
0080 PVS_
0081     Small number of special cased geometries created locally, eg RaindropRockAirWater
0082 PVG_
0083     load geometry from a gdmlpath using spath::GDMLPathFromGEOM method
0084 PVP_
0085     PMTSim getter, requiring WITH_PMTSIM macro to be set meaning that PMTSim pkg was found by CMake
0086 PVL_
0087     For names starting with List, multiple volumes are created and arranged eg into a grid.
0088     The names to create are extracted using the name argument as a comma delimited list
0089 PV1_
0090     Places single lv once or duplicates it several times depending on
0091     strings found in the name such as Grid,Cube,Xoff,Yoff,Zoff
0092 
0093 **/
0094 
0095 const G4VPhysicalVolume* U4VolumeMaker::PV(){ return PV(GEOM); }
0096 const G4VPhysicalVolume* U4VolumeMaker::PV(const char* name)
0097 {
0098     std::cerr << "U4VolumeMaker::PV name " << name << std::endl ;
0099     LOG(LEVEL) << "[" << name ;
0100     const G4VPhysicalVolume* pv = nullptr ;
0101     if(pv == nullptr) pv = PVS_(name);
0102     if(pv == nullptr) pv = PVG_(name);
0103     if(pv == nullptr) pv = PVP_(name);
0104     if(pv == nullptr) pv = PVL_(name);
0105     if(pv == nullptr) pv = PV1_(name);
0106     LOG_IF(error, pv == nullptr) << "returning nullptr for name [" << name << "]" ;
0107     LOG(LEVEL) << "]" << name ;
0108     return pv ;
0109 }
0110 
0111 /**
0112 U4VolumeMaker::PVG_
0113 ---------------------
0114 
0115 Attempts to load geometry from a gdmlpath using the spath::GDMLPathFromGEOM
0116 resolution method.
0117 
0118 Additionally using spath::GDMLSub resolution when hello_GDMLSub
0119 envvar exists the string obtained is used to select a PV from within the loaded
0120 tree of volumes.  For example with GEOM of J000 the below envvar is checked::
0121 
0122     export J000_GDMLSub=HamamatsuR12860sMask_virtual0x:0:1000
0123 
0124 Note that while loading GDML and then selecting a PV from it is a useful capability
0125 it is generally only useful for a first look at an issue or to isolate an issue.
0126 As typically to understand what is going wrong with a geometry it is necessary
0127 to iterate making changes to the geometry. In order to do that it is necessary
0128 to take control of the geometry defining code for example in j/PMTSim.
0129 
0130 **/
0131 
0132 const char* U4VolumeMaker::PVG_WriteNames     = "U4VolumeMaker_PVG_WriteNames" ;
0133 const char* U4VolumeMaker::PVG_WriteNames_Sub = "U4VolumeMaker_PVG_WriteNames_Sub" ;
0134 
0135 const G4VPhysicalVolume* U4VolumeMaker::PVG_(const char* name)
0136 {
0137     METH = "PVG_" ;
0138     const char* gdmlpath = spath::GDMLPathFromGEOM(name) ;
0139     bool exists = gdmlpath && spath::Exists(gdmlpath) ;
0140 
0141     const char* sub = spath::GEOMSub(name);
0142 
0143     std::cerr
0144         << "U4VolumeMaker::PVG_"
0145         << " name " << name
0146         << " gdmlpath " << ( gdmlpath ? gdmlpath : "-" )
0147         << " sub " << ( sub ? sub : "-" )
0148         << " exists " << exists
0149         << std::endl
0150         ;
0151 
0152 
0153     const G4VPhysicalVolume* loaded = exists ? U4GDML::Read(gdmlpath) : nullptr ;
0154 
0155     if(loaded && ssys::getenvbool(PVG_WriteNames))
0156         U4Volume::WriteNames( loaded, spath::Resolve("$TMP", PVG_WriteNames));
0157 
0158     const G4VPhysicalVolume* pv = loaded ;
0159 
0160     if( loaded && sub )
0161     {
0162         const G4VPhysicalVolume* pv_sub = U4Volume::FindPVSub( loaded, sub ) ;
0163         G4LogicalVolume* lv_sub = pv_sub->GetLogicalVolume();
0164         std::vector<G4LogicalVolume*> lvs = {lv_sub} ;
0165         pv = Wrap( name, lvs );
0166 
0167         if(ssys::getenvbool(PVG_WriteNames_Sub))
0168             U4Volume::WriteNames( pv, spath::Resolve("$TMP", PVG_WriteNames_Sub));
0169     }
0170 
0171     LOG(LEVEL)
0172           << " name " << name
0173           << " gdmlpath " << gdmlpath
0174           << " exists " << exists
0175           << " loaded " << loaded
0176           << " sub " << ( sub ? sub : "-" )
0177           << " pv " << pv
0178           ;
0179 
0180     return pv ;
0181 }
0182 
0183 
0184 /**
0185 U4VolumeMaker::PVP_ : Get PMTSim PV
0186 ------------------------------------
0187 
0188 PMTSim::HasManagerPrefix
0189     returns true for names starting with one of: hama, nnvt, hmsk, nmsk, lchi
0190 
0191 PMTSim::GetLV PMTSim::GetPV
0192      these methods act as go betweens to the underlying managers with prefixes
0193      that identify the managers offset
0194 
0195 
0196 TODO: need to generalize the wrapping
0197 
0198 **/
0199 
0200 const G4VPhysicalVolume* U4VolumeMaker::PVP_(const char* name)
0201 {
0202     METH = "PVP_" ;
0203     const G4VPhysicalVolume* pv = nullptr ;
0204 #ifdef WITH_PMTSIM
0205     const char* geomlist = spath::GEOMList(name);   // consult envvar name_GEOMList
0206     std::vector<std::string> names ;
0207     if( geomlist == nullptr )
0208     {
0209         names.push_back(name);
0210     }
0211     else
0212     {
0213         sstr::Split(geomlist, ',', names );
0214     }
0215     int num_names = names.size();
0216 
0217     LOG(LEVEL)
0218          << "[ WITH_PMTSIM"
0219          << " geomlist [" << ( geomlist ? geomlist : "-" ) << "] "
0220          << " name [" << name << "] "
0221          << " num_names " << num_names
0222          ;
0223 
0224     std::vector<G4LogicalVolume*> lvs ; ;
0225 
0226     for(int i=0 ; i < num_names ; i++)
0227     {
0228         const char* n = names[i].c_str();
0229         bool has_manager_prefix = PMTSim::HasManagerPrefix(n) ;
0230         LOG(LEVEL) << "[ WITH_PMTSIM n [" << n << "] has_manager_prefix " << has_manager_prefix ;
0231 
0232         if(!has_manager_prefix) return nullptr ;
0233         //assert( has_manager_prefix );
0234 
0235         G4LogicalVolume* lv = PMTSim::GetLV(n) ;
0236         LOG_IF(fatal, lv == nullptr ) << "PMTSim::GetLV returned nullptr for n [" << n << "]" ;
0237         assert( lv );
0238 
0239         lvs.push_back(lv);
0240     }
0241 
0242     pv = Wrap( name, lvs ) ;
0243 
0244     LOG(LEVEL) << "]" ;
0245 #else
0246     LOG(info) << " not-WITH_PMTSIM name [" << name << "]" ;
0247 #endif
0248     return pv ;
0249 }
0250 
0251 
0252 /**
0253 U4VolumeMaker::PVS_
0254 ----------------------
0255 
0256 Special names with local C++ implementations of geometry
0257 
0258 **/
0259 
0260 const G4VPhysicalVolume* U4VolumeMaker::PVS_(const char* name)
0261 {
0262     METH = "PVS_" ;
0263     const G4VPhysicalVolume* pv = nullptr ;
0264     if(     strcmp(name,"BoxOfScintillator" ) == 0)      pv = BoxOfScintillator(1000.);
0265     else if(strcmp(name,"RaindropRockAirWater" ) == 0)   pv = RaindropRockAirWater(false);
0266     else if(strcmp(name,"RaindropRockAirWaterSD" ) == 0) pv = RaindropRockAirWater(true);
0267     else if(strcmp(name,"BigWaterPool" ) == 0)           pv = BigWaterPool();
0268     else if(sstr::StartsWith(name,"LocalFastenerAcrylicConstruction" )) pv = LocalFastenerAcrylicConstruction(name);
0269     else if(sstr::StartsWith(name,"Local" ))                            pv = Local(name);
0270     return pv ;
0271 }
0272 const G4VPhysicalVolume* U4VolumeMaker::PVL_(const char* name)
0273 {
0274     METH = "PVL_" ;
0275     if(!sstr::StartsWith(name, "List")) return nullptr  ;
0276     std::vector<G4LogicalVolume*> lvs ;
0277     LV(lvs, name + strlen("List") );
0278     const G4VPhysicalVolume* pv = WrapLVGrid(lvs, 1, 1, 1 );
0279     return pv ;
0280 }
0281 
0282 /**
0283 U4VolumeMaker::PV1_
0284 -------------------
0285 
0286 name specifies single LV for wrapping to form the PV
0287 
0288 **/
0289 const G4VPhysicalVolume* U4VolumeMaker::PV1_(const char* name)
0290 {
0291     METH = "PV1_" ;
0292     const char* matname_ = nullptr ;  // look for material in name or default to Vacuum
0293     G4LogicalVolume* lv = LV(name, matname_) ;
0294     LOG_IF(error, lv == nullptr) << " failed to access lv for name " << name ;
0295     if(lv == nullptr) return nullptr ;
0296 
0297     const G4VPhysicalVolume* pv = nullptr ;
0298     bool grid = strstr(name, "Grid") != nullptr ;
0299     bool cube = strstr(name, "Cube") != nullptr ;
0300     bool xoff = strstr(name, "Xoff") != nullptr ;
0301     bool yoff = strstr(name, "Yoff") != nullptr ;
0302     bool zoff = strstr(name, "Zoff") != nullptr ;
0303     bool simp = strstr(name, "Simple") != nullptr ;
0304 
0305     if(grid)      pv =   WrapLVGrid(     lv, 1, 1, 1 );
0306     else if(cube) pv =   WrapLVCube(     lv, 100., 100., 100. );
0307     else if(xoff) pv =   WrapLVOffset(   lv, 200.,   0.,   0. );
0308     else if(yoff) pv =   WrapLVOffset(   lv,   0., 200.,   0. );
0309     else if(zoff) pv =   WrapLVOffset(   lv,   0.,   0., 200. );
0310     else if(simp) pv =   WrapLVSimple(   lv );
0311     else          pv =   WrapLVCube(     lv,   0.,   0.,   0. );
0312     return pv ;
0313 }
0314 
0315 
0316 /**
0317 U4VolumeMaker::LV
0318 ---------------------------
0319 
0320 Note the generality:
0321 
0322 * U4SolidMaker::Make can create many different solid shapes based on the start of the name
0323 * U4Material::FindMaterialName extracts the material name by looking for material names within *name*
0324 
0325 HMM: maybe provide PMTSim LV this way too ? Based on some prefix ?
0326 
0327 **/
0328 
0329 G4LogicalVolume*  U4VolumeMaker::LV(const char* name, const char* matname_)
0330 {
0331     const G4VSolid* solid_  = U4SolidMaker::Make(name);
0332     LOG_IF(error, solid_==nullptr) << " failed to access solid for name " << name ;
0333     if(solid_ == nullptr) return nullptr ;
0334 
0335     G4VSolid* solid = const_cast<G4VSolid*>(solid_);
0336     const char* matname = matname_ ? matname_ : U4Material::FindMaterialName(name) ;
0337     G4Material* material = U4Material::Get( matname ? matname : U4Material::VACUUM );
0338     G4LogicalVolume* lv = new G4LogicalVolume( solid, material, name, 0,0,0,true );
0339     return lv ;
0340 }
0341 
0342 
0343 
0344 
0345 
0346 /**
0347 U4VolumeMaker::LV vector interface : creates multiple LV using delimited names string
0348 --------------------------------------------------------------------------------------
0349 **/
0350 
0351 void U4VolumeMaker::LV(std::vector<G4LogicalVolume*>& lvs , const char* names_, char delim )
0352 {
0353     std::vector<std::string> names ;
0354     sstr::Split(names_ , delim, names );
0355     unsigned num_names = names.size();
0356     LOG(LEVEL) << " names_ " << names_ << " num_names " << num_names ;
0357     assert( num_names > 1 );
0358 
0359     for(unsigned i=0 ; i < num_names ; i++)
0360     {
0361         const char* name = names[i].c_str();
0362         const char* matname_ = nullptr ;
0363         G4LogicalVolume* lv = LV(name, matname_ ) ;
0364         assert(lv);
0365         lvs.push_back(lv);
0366     }
0367 }
0368 
0369 
0370 /**
0371 U4VolumeMaker::Wrap
0372 -----------------------
0373 
0374 Consults envvar ${GEOM}_GEOMWrap for the wrap config,
0375 which when present must be one of : AroundCylinder, AroundSphere, AroundCircle
0376 
0377 **/
0378 
0379 const G4VPhysicalVolume* U4VolumeMaker::Wrap( const char* name, std::vector<G4LogicalVolume*>& items_lv )
0380 {
0381     const char* wrap = spath::GEOMWrap(name);
0382     LOG(LEVEL) << "[ name " << name << " GEOMWrap " << ( wrap ? wrap : "-" ) ;
0383     std::cerr << "U4VolumeMaker::Wrap "  << "[ name " << name << " GEOMWrap " << ( wrap ? wrap : "-" ) << std::endl ;
0384     const G4VPhysicalVolume* pv = wrap == nullptr ? WrapRockWater( items_lv ) : WrapAroundItem( name, items_lv, wrap );
0385     LOG(LEVEL) << "] name " << name << " wrap " << ( wrap ? wrap : "-" ) << " pv " << ( pv ? "YES" : "NO" ) ;
0386     return pv ;
0387 }
0388 
0389 /**
0390 U4VolumeMaker::WrapRockWater
0391 -------------------------------
0392 
0393     +---------------Rock------------+
0394     |                               |
0395     |                               |
0396     |     +--------Water------+     |
0397     |     |                   |     |
0398     |     |                   |     |
0399     |     |      +-Item-+     |     |
0400     |     |      |      |     |     |
0401     |     |      |      |     |     |
0402     |     |      +------+     |     |
0403     |     |                   |     |
0404     |     |                   |     |
0405     |     +-------------------+     |
0406     |                               |
0407     |                               |
0408     +-------------------------------+
0409 
0410 
0411 **/
0412 
0413 const G4VPhysicalVolume* U4VolumeMaker::WrapRockWater( std::vector<G4LogicalVolume*>& items_lv )
0414 {
0415     assert( items_lv.size() >= 1 );
0416     G4LogicalVolume* item_lv = items_lv[items_lv.size()-1] ;
0417 
0418     LOG(LEVEL) << "[ items_lv.size " << items_lv.size()   ;
0419     std::cerr << "[ items_lv.size " << items_lv.size() << std::endl  ;
0420 
0421     double rock_halfside  = ssys::getenv_<double>(U4VolumeMaker_WrapRockWater_Rock_HALFSIDE , 1000.);
0422     double water_halfside = ssys::getenv_<double>(U4VolumeMaker_WrapRockWater_Water_HALFSIDE,  900. );
0423 
0424 
0425 
0426     std::vector<double>* _boxscale = ssys::getenv_vec<double>(U4VolumeMaker_WrapRockWater_BOXSCALE, "1,1,1" );
0427     if(_boxscale) assert(_boxscale->size() == 3 );
0428     const double* boxscale = _boxscale ? _boxscale->data() : nullptr ;
0429 
0430     LOG(LEVEL) << U4VolumeMaker_WrapRockWater_Rock_HALFSIDE << " " << rock_halfside ;
0431     LOG(LEVEL) << U4VolumeMaker_WrapRockWater_Water_HALFSIDE << " " << water_halfside ;
0432 
0433     std::cerr << U4VolumeMaker_WrapRockWater_Rock_HALFSIDE << " " << rock_halfside << std::endl ;
0434     std::cerr << U4VolumeMaker_WrapRockWater_Water_HALFSIDE << " " << water_halfside << std::endl ;
0435 
0436 
0437     G4LogicalVolume*  rock_lv  = Box_(rock_halfside,  "Rock",  nullptr, boxscale );
0438     G4LogicalVolume*  water_lv = Box_(water_halfside, "Water", nullptr, boxscale );
0439 
0440     //const char* flip_axes = "Z" ;
0441     const char* flip_axes = nullptr ;
0442     const G4VPhysicalVolume* item_pv  = Place(item_lv,  water_lv, flip_axes );  assert( item_pv );
0443 
0444 
0445     LOG(LEVEL) << std::endl << U4Volume::Traverse( item_pv );
0446 
0447     std::vector<std::string>* vbs1 = nullptr ;
0448     vbs1 = ssys::getenv_vec<std::string>(U4VolumeMaker_WrapRockWater_BS1, "pyrex_vacuum_bs:hama_body_phys:hama_inner1_phys", ':' );
0449 
0450     G4LogicalBorderSurface* bs1 = nullptr ;
0451     if(vbs1 && vbs1->size() == 3 )
0452     {
0453          const std::string& bs1n = (*vbs1)[0] ;
0454          const std::string& pv1 = (*vbs1)[1] ;
0455          const std::string& pv2 = (*vbs1)[2] ;
0456          bs1 = U4Surface::MakePerfectDetectorBorderSurface(bs1n.c_str(), pv1.c_str(), pv2.c_str(), item_pv );
0457 
0458         LOG(error)
0459             << " attempt to add  PerfectDetectorBorderSurface between volumes "
0460             << " bs1n " << bs1n
0461             << " pv1 " << pv1
0462             << " pv2 " << pv2
0463             << " bs1 " << ( bs1 ? "YES" : "NO" )
0464             << " using config from " << U4VolumeMaker_WrapRockWater_BS1
0465             ;
0466     }
0467     //if(bs1 == nullptr) std::raise(SIGINT);
0468 
0469 
0470     const G4VPhysicalVolume* water_pv = Place(water_lv,  rock_lv);  assert( water_pv );
0471     const G4VPhysicalVolume* rock_pv  = Place(rock_lv,  nullptr );
0472 
0473     G4LogicalBorderSurface* water_rock_bs = U4Surface::MakePerfectAbsorberBorderSurface("water_rock_bs", water_pv, rock_pv );
0474     assert( water_rock_bs );
0475     if(!water_rock_bs) std::raise(SIGINT);
0476 
0477     LOG(LEVEL) << "]"  ;
0478     return rock_pv ;
0479 }
0480 
0481 /**
0482 U4VolumeMaker::WrapAroundItem
0483 -------------------------------
0484 
0485 The *item_lv* is repeated many times using transforms from U4VolumeMaker::MakeTransforms
0486 controlled by the *prefix* string which must be one of::
0487 
0488     AroundSphere
0489     AroundCylinder
0490     AroundCircle
0491 
0492 which makes use of SPlace::AroundSphere SPlace::AroundCylinder SPlace::AroundCircle.
0493 All those repeats have a "Water" box mother volume which is contained within "Rock".
0494 
0495 **/
0496 
0497 NP* U4VolumeMaker::TRS = nullptr ;
0498 
0499 const G4VPhysicalVolume* U4VolumeMaker::WrapAroundItem( const char* name, std::vector<G4LogicalVolume*>& items_lv, const char* prefix )
0500 {
0501     assert( items_lv.size() >= 1 );
0502 
0503     NP* trs = MakeTransforms(name, prefix) ;
0504     TRS = trs ;
0505 
0506     LOG(LEVEL)
0507         << " items_lv.size " << items_lv.size()
0508         << " prefix " << prefix
0509         << " trs " << ( trs ? trs->sstr() : "-" )
0510         ;
0511 
0512 
0513     double rock_halfside   = ssys::getenv_<double>(U4VolumeMaker_WrapAroundItem_Rock_HALFSIDE,  20000.);
0514     double water_halfside  = ssys::getenv_<double>(U4VolumeMaker_WrapAroundItem_Water_HALFSIDE, 19000.);
0515 
0516     std::vector<double>* _rock_boxscale = ssys::getenv_vec<double>(U4VolumeMaker_WrapAroundItem_Rock_BOXSCALE, "1,1,1" );
0517     std::vector<double>* _water_boxscale = ssys::getenv_vec<double>(U4VolumeMaker_WrapAroundItem_Water_BOXSCALE, "1,1,1" );
0518 
0519     if(_rock_boxscale) assert(_rock_boxscale->size() == 3 );
0520     if(_water_boxscale) assert(_water_boxscale->size() == 3 );
0521 
0522     const double* rock_boxscale = _rock_boxscale ? _rock_boxscale->data() : nullptr ;
0523     const double* water_boxscale = _water_boxscale ? _water_boxscale->data() : nullptr ;
0524 
0525     G4LogicalVolume*  rock_lv  = Box_(rock_halfside,  "Rock" , nullptr, rock_boxscale );
0526     G4LogicalVolume*  water_lv = Box_(water_halfside, "Water", nullptr, water_boxscale );
0527 
0528     WrapAround(prefix, trs, items_lv, water_lv );
0529     // item_lv placed inside water_lv once for each transform
0530 
0531     const G4VPhysicalVolume* water_pv = Place(water_lv,  rock_lv);
0532     assert( water_pv );
0533     if(!water_pv) std::raise(SIGINT);
0534 
0535     const G4VPhysicalVolume* rock_pv  = Place(rock_lv,  nullptr );
0536 
0537     return rock_pv ;
0538 }
0539 
0540 
0541 
0542 NP* U4VolumeMaker::GetTransforms() // static
0543 {
0544     return TRS ;
0545 }
0546 
0547 /**
0548 U4VolumeMaker::SaveTransforms
0549 -----------------------------
0550 
0551 **/
0552 
0553 
0554 void U4VolumeMaker::SaveTransforms( const char* savedir ) // static
0555 {
0556     NP* TRS = U4VolumeMaker::TRS ;
0557     if(TRS) TRS->save(savedir, "TRS.npy") ;
0558 }
0559 
0560 
0561 
0562 /**
0563 U4VolumeMaker::WrapVacuum
0564 ----------------------------
0565 
0566 The LV provided is placed within a WorldBox of halfside extent and the world PV is returned.
0567 
0568 **/
0569 
0570 const G4VPhysicalVolume* U4VolumeMaker::WrapVacuum( G4LogicalVolume* item_lv )
0571 {
0572     double halfside = ssys::getenv_<double>(U4VolumeMaker_WrapVacuum_HALFSIDE, 1000.);
0573     LOG(LEVEL) << U4VolumeMaker_WrapVacuum_HALFSIDE << " " << halfside ;
0574 
0575     G4LogicalVolume*   vac_lv  = Box_(halfside, U4Material::VACUUM );
0576 
0577     const G4VPhysicalVolume* item_pv = Place(item_lv, vac_lv);
0578     assert( item_pv );
0579     if(!item_pv) std::raise(SIGINT);
0580 
0581     const G4VPhysicalVolume* vac_pv  = Place(vac_lv, nullptr );
0582 
0583     return vac_pv ;
0584 }
0585 
0586 
0587 
0588 const G4VPhysicalVolume* U4VolumeMaker::WrapLVSimple( G4LogicalVolume* item_lv )
0589 {
0590     double halfside = ssys::getenv_<double>(U4VolumeMaker_WrapLVSimple_HALFSIDE, 1000.);
0591     LOG(LEVEL) << U4VolumeMaker_WrapLVSimple_HALFSIDE << " " << halfside ;
0592 
0593     G4LogicalVolume*   vac_lv  = Box_(halfside, U4Material::VACUUM );
0594 
0595     const G4VPhysicalVolume* item_pv = Place(item_lv, vac_lv);
0596     assert( item_pv );
0597     if(!item_pv) std::raise(SIGINT);
0598 
0599     const G4VPhysicalVolume* vac_pv  = Place(vac_lv, nullptr );
0600 
0601     return vac_pv ;
0602 }
0603 
0604 
0605 
0606 
0607 
0608 
0609 /**
0610 U4VolumeMaker::WrapLVGrid
0611 ---------------------------
0612 
0613 Returns a physical volume with the argument lv placed multiple times
0614 in a grid specified by (nx,ny,nz) integers. (1,1,1) yields 3x3x3 grid.
0615 
0616 The vector argument method places the lv in a grid within a box.
0617 Grid ranges::
0618 
0619     -nx:nx+1
0620     -ny:ny+1
0621     -nz:nz+1
0622 
0623 Example (nx,ny,nz):
0624 
0625 1,1,1
0626      yields a grid with 3 elements on each side, for 3*3*3=27
0627 0,0,0
0628      yields a single element
0629 
0630 **/
0631 
0632 const G4VPhysicalVolume* U4VolumeMaker::WrapLVGrid( G4LogicalVolume* lv, int nx, int ny, int nz  )
0633 {
0634     std::vector<G4LogicalVolume*> lvs ;
0635     lvs.push_back(lv);
0636     return WrapLVGrid(lvs, nx, ny, nz ) ;
0637 }
0638 
0639 const G4VPhysicalVolume* U4VolumeMaker::WrapLVGrid( std::vector<G4LogicalVolume*>& lvs, int nx, int ny, int nz  )
0640 {
0641     unsigned num_lv = lvs.size();
0642 
0643     int extent = std::max(std::max(nx, ny), nz);
0644     G4double sc = 500. ;
0645     G4double halfside = sc*extent*3. ;
0646 
0647     LOG(LEVEL)
0648         << " num_lv " << num_lv
0649         << " (nx,ny,nz) (" << nx << "," << ny << " " << nz << ")"
0650         << " extent " << extent
0651         << " halfside " << halfside
0652         ;
0653 
0654     const G4VPhysicalVolume* world_pv = WorldBox(halfside);
0655     G4LogicalVolume* world_lv = world_pv->GetLogicalVolume();
0656 
0657     unsigned count = 0 ;
0658     for(int ix=-nx ; ix < nx+1 ; ix++)
0659     for(int iy=-ny ; iy < ny+1 ; iy++)
0660     for(int iz=-nz ; iz < nz+1 ; iz++)
0661     {
0662         G4LogicalVolume* ulv = lvs[count%num_lv] ;
0663         count += 1 ;
0664         assert( ulv );
0665 
0666          const G4String& ulv_name_ = ulv->GetName();
0667          std::string ulv_name = ulv_name_ ;
0668          ulv_name += "_plc_" ;
0669 
0670         const char* iname = GridName( ulv_name.c_str(), ix, iy, iz, "" );
0671         G4ThreeVector tla( sc*double(ix), sc*double(iy), sc*double(iz) );
0672 
0673         LOG(LEVEL)
0674            << " G4PVPlacement "
0675            << " count " << std::setw(3) << count
0676            << Desc(tla)
0677            <<  iname
0678            ;
0679         const G4VPhysicalVolume* pv_n = new G4PVPlacement(0, tla, ulv ,iname,world_lv,false,0);
0680         assert( pv_n );
0681         if(!pv_n) std::raise(SIGINT);
0682     }
0683     return world_pv ;
0684 }
0685 
0686 
0687 
0688 
0689 /**
0690 U4VolumeMaker::WrapLVOffset
0691 -----------------------------
0692 
0693 This is used from PV1_ for Xoff Yoff Zoff names
0694 
0695 1. use maximum of tx,ty,tz to define world box halfside
0696 2. places lv within world volume
0697 3. adds BoxMinusOrb at origin
0698 
0699 **/
0700 
0701 const G4VPhysicalVolume* U4VolumeMaker::WrapLVOffset( G4LogicalVolume* lv, double tx, double ty, double tz )
0702 {
0703     double halfside = 3.*std::max( std::max( tx, ty ), tz );
0704     assert( halfside > 0. );
0705 
0706     const G4VPhysicalVolume* world_pv = WorldBox(halfside);
0707     G4LogicalVolume* world_lv = world_pv->GetLogicalVolume();
0708 
0709     AddPlacement(world_lv, lv, tx, ty, tz );
0710 
0711     bool bmo = true ;
0712     if(bmo) AddPlacement( world_lv, "BoxMinusOrb", 0., 0., 0. );
0713 
0714     return world_pv ;
0715 }
0716 
0717 /**
0718 U4VolumeMaker::WrapLVCube
0719 -------------------------------
0720 
0721 This if used from PV1_ for "Cube" string.
0722 Places the lv at 8 positions at the corners of a cube.
0723 
0724   ZYX
0725 0 000
0726 1 001
0727 2 010
0728 3 011
0729 4 100
0730 5 101
0731 6 110
0732 7 111
0733 
0734               110             111
0735                 +-----------+
0736                /|          /|
0737               / |         / |
0738              /  |        /  |
0739             +-----------+   |
0740             |   +-------|---+ 011
0741             |  /        |  /
0742             | /         | /
0743             |/          |/
0744             +-----------+
0745           000          001
0746 
0747    Z   Y
0748    | /
0749    |/
0750    0-- X
0751 
0752 **/
0753 
0754 const G4VPhysicalVolume* U4VolumeMaker::WrapLVCube( G4LogicalVolume* lv, double tx, double ty, double tz )
0755 {
0756     double halfside = 3.*std::max( std::max( tx, ty ), tz );
0757 
0758     const G4VPhysicalVolume* world_pv = WorldBox(halfside);
0759     G4LogicalVolume* world_lv = world_pv->GetLogicalVolume();
0760     G4String name = lv->GetName();
0761 
0762     for(unsigned i=0 ; i < 8 ; i++)
0763     {
0764         bool px = ((i & 0x1) != 0) ;
0765         bool py = ((i & 0x2) != 0) ;
0766         bool pz = ((i & 0x4) != 0) ;
0767         G4ThreeVector tla( px ? tx : -tx ,  py ? ty : -ty,  pz ? tz : -tz );
0768         const char* iname = GridName(name.c_str(), int(px), int(py), int(pz), "" );
0769         const G4VPhysicalVolume* pv = new G4PVPlacement(0,tla,lv,iname,world_lv,false,0);
0770         assert( pv );
0771         if(!pv) std::raise(SIGINT);
0772     }
0773     return world_pv ;
0774 }
0775 
0776 
0777 
0778 
0779 
0780 
0781 
0782 /**
0783 U4VolumeMaker::AddPlacement : Translation places *lv* within *mother_lv*
0784 ---------------------------------------------------------------------------
0785 
0786 Used from U4VolumeMaker::WrapLVOffset
0787 
0788 **/
0789 
0790 const G4VPhysicalVolume* U4VolumeMaker::AddPlacement( G4LogicalVolume* mother_lv,  G4LogicalVolume* lv,  double tx, double ty, double tz )
0791 {
0792     G4ThreeVector tla(tx,ty,tz);
0793     const char* pv_name = spath::Name( lv->GetName().c_str(), "_placement" );
0794     LOG(LEVEL) << Desc(tla) << " " << pv_name ;
0795     const G4VPhysicalVolume* pv = new G4PVPlacement(0,tla,lv,pv_name,mother_lv,false,0);
0796     return pv ;
0797 }
0798 const G4VPhysicalVolume* U4VolumeMaker::AddPlacement( G4LogicalVolume* mother, const char* name,  double tx, double ty, double tz )
0799 {
0800     const char* matname_ = nullptr ;
0801     G4LogicalVolume* lv = LV(name, matname_);
0802     return AddPlacement( mother, lv, tx, ty, tz );
0803 }
0804 
0805 const char* U4VolumeMaker::GridName(const char* prefix, int ix, int iy, int iz, const char* suffix)
0806 {
0807     std::stringstream ss ;
0808     ss << prefix << ix << "_" << iy << "_" << iz << suffix ;
0809     std::string s = ss.str();
0810     return strdup(s.c_str());
0811 }
0812 
0813 const char* U4VolumeMaker::PlaceName(const char* prefix, int ix, const char* suffix)
0814 {
0815     std::stringstream ss ;
0816     if(prefix) ss << prefix ;
0817     ss << ix ;
0818     if(suffix) ss << suffix ;
0819     std::string s = ss.str();
0820     return strdup(s.c_str());
0821 }
0822 
0823 
0824 
0825 
0826 
0827 /**
0828 U4VolumeMaker::WorldBox
0829 --------------------------
0830 
0831 Used from U4VolumeMaker::WrapLVOffset U4VolumeMaker::WrapLVCube
0832 
0833 **/
0834 
0835 const G4VPhysicalVolume* U4VolumeMaker::WorldBox( double halfside, const char* mat )
0836 {
0837     if(mat == nullptr) mat = U4Material::VACUUM ;
0838     return Box(halfside, mat, "World", nullptr );
0839 }
0840 const G4VPhysicalVolume* U4VolumeMaker::BoxOfScintillator( double halfside )
0841 {
0842     return BoxOfScintillator(halfside, "BoxOfScintillator", nullptr );
0843 }
0844 const G4VPhysicalVolume* U4VolumeMaker::BoxOfScintillator( double halfside, const char* prefix, G4LogicalVolume* mother_lv )
0845 {
0846     return Box(halfside, U4Material::SCINTILLATOR, "BoxOfScintillator", mother_lv);
0847 }
0848 
0849 
0850 
0851 const G4VPhysicalVolume* U4VolumeMaker::Box(double halfside, const char* mat, const char* prefix, G4LogicalVolume* mother_lv )
0852 {
0853     if(prefix == nullptr) prefix = mat ;
0854     G4LogicalVolume* lv = Box_(halfside, mat, prefix);
0855     return Place(lv, mother_lv);
0856 }
0857 
0858 
0859 const G4VPhysicalVolume* U4VolumeMaker::Place( G4LogicalVolume* lv, G4LogicalVolume* mother_lv, const char* flip_axes )
0860 {
0861     const char* lv_name = lv->GetName().c_str() ;
0862     const char* pv_name = spath::Name(lv_name, "_pv") ;
0863 
0864     U4RotationMatrix* flip = flip_axes ? U4RotationMatrix::Flip(flip_axes) : nullptr ;
0865     return new G4PVPlacement(flip,G4ThreeVector(), lv, pv_name, mother_lv, false, 0);
0866 }
0867 
0868 
0869 const G4VPhysicalVolume* U4VolumeMaker::LocalFastenerAcrylicConstruction(const char* name)
0870 {
0871     const char* PREFIX = "LocalFastenerAcrylicConstruction" ;
0872     assert( sstr::StartsWith(name,PREFIX ));
0873     int num = strlen(name) > strlen(PREFIX) ? std::atoi( name + strlen(PREFIX) ) : 8 ;
0874 
0875     LOG(info)
0876         << " name " <<  ( name ? name : "-" )
0877         << " num " << num
0878         ;
0879 
0880     G4double universe_halfside = 300.*mm ;
0881     const char* universe_material = "G4_AIR" ;
0882 
0883     G4LogicalVolume* universe_lv  = Box_(universe_halfside, universe_material );
0884     G4LogicalVolume* object_lv = LV(name, "G4_Pb" ) ;
0885 
0886     const G4VPhysicalVolume* universe_pv = new G4PVPlacement(0,G4ThreeVector(),  universe_lv ,  "universe_pv", nullptr,false,0);
0887     [[maybe_unused]] const G4VPhysicalVolume* object_pv = new G4PVPlacement(0,G4ThreeVector(), object_lv ,"object_pv", universe_lv,false,0);
0888     assert( object_pv );
0889 
0890     return universe_pv ;
0891 }
0892 
0893 
0894 /**
0895 U4VolumeMaker::Local
0896 ----------------------
0897 
0898 Local prefixed names like "LocalOuterReflectorOrbSubtraction" are passed to U4SolidMaker
0899 to construct solids without the prefix eg: "OuterReflectorOrbSubtraction".
0900 Then the bounding limits are used to get the extent of the solid which is used to
0901 configure the size of the universe volume. Note that this functionality can
0902 be used from the below script that creates the Geant4 geometry and then translates
0903 into an Opticks one which is persisted.
0904 
0905     ~/o/g4cx/tests/G4CX_U4TreeCreateCSGFoundryTest.sh
0906 
0907 Thence can use the below tools to visualize in various ways::
0908 
0909     cxr_min.sh  # ray trace render
0910     cxt_min.sh  # simtrace cross sections
0911     ssst.sh     # triangulated SScene viz
0912 
0913 
0914 Consider phicut applied to a shape, which results in
0915 geometry offset from the origin::
0916 
0917 
0918             +
0919            /   +
0920           /   /
0921          /   /
0922         +   /
0923            +
0924 
0925 
0926                        +
0927 
0928 
0929 Using s_bb::Extent to define the universe size can result
0930 in too small a universe, so use s_bb::AbsMax
0931 
0932 
0933 **/
0934 
0935 const G4VPhysicalVolume* U4VolumeMaker::Local(const char* name)
0936 {
0937     const char* PREFIX = "Local" ;
0938     assert( sstr::StartsWith(name,PREFIX ));
0939     const char* name_after_PREFIX = name + strlen(PREFIX);
0940 
0941     G4LogicalVolume* object_lv = LV(name_after_PREFIX, "G4_Pb" ) ;
0942     const G4VSolid* const solid = object_lv->GetSolid();
0943 
0944     G4ThreeVector pMin ;
0945     G4ThreeVector pMax ;
0946     solid->BoundingLimits(pMin, pMax);
0947 
0948 
0949     std::array<double,6> bb = {pMin.x(), pMin.y(), pMin.z(), pMax.x(), pMax.y(), pMax.z() };
0950     double extent = s_bb::Extent<double>(bb.data());
0951     double absmax = s_bb::AbsMax<double>(bb.data());
0952     LOG(info)
0953         << " name " <<  ( name ? name : "-" )
0954         << " name_after_PREFIX " <<  ( name_after_PREFIX ? name_after_PREFIX : "-" )
0955         << " bb " << s_bb::Desc<double>(bb.data())
0956         << " extent " << extent
0957         << " absmax " << absmax
0958         ;
0959 
0960     G4double universe_halfside = 1.5*absmax*mm ;  // extent not appropriate for geometry offset from origin
0961     const char* universe_material = "G4_AIR" ;
0962 
0963     G4LogicalVolume* universe_lv  = Box_(universe_halfside, universe_material );
0964 
0965     const G4VPhysicalVolume* universe_pv = new G4PVPlacement(0,G4ThreeVector(),  universe_lv ,  "universe_pv", nullptr,false,0);
0966     [[maybe_unused]] const G4VPhysicalVolume* object_pv = new G4PVPlacement(0,G4ThreeVector(), object_lv ,"object_pv", universe_lv,false,0);
0967     assert( object_pv );
0968 
0969     return universe_pv ;
0970 }
0971 
0972 
0973 
0974 
0975 
0976 
0977 
0978 /**
0979 U4VolumeMaker::RaindropRockAirWater
0980 -------------------------------------
0981 
0982 cf CSG/CSGMaker.cc CSGMaker::makeBoxedSphere
0983 
0984 
0985    +------------------------+
0986    | Rock/Container         |
0987    |    +-----------+       |
0988    |    | Air/Medium|       |
0989    |    |    . .    |       |
0990    |    |   .Wa/Drop|       |
0991    |    |    . .    |       |
0992    |    |           |       |
0993    |    +-----|-----+       |
0994    |                        |
0995    +----------|--rock_halfs-+
0996   -X                        +X
0997 
0998 
0999 
1000 Defaults::
1001 
1002     HALFSIDE: 100.
1003     FACTOR: 1.
1004 
1005     drop_radius        :  halfside/2.         : 50.
1006     medium_halfside    :  halfside*factor     : 100.
1007     container_halfside :  1.1*halfside*factor : 110.
1008 
1009     top of drop   (0,0, 50)
1010     bot of drop   (0,0,-50)
1011     left of drop  (-50,0,0)
1012     right of drop ( 50,0,0)
1013 
1014 
1015 An easy way to get some scattering and absorption to happen
1016 is to increase U4VolumeMaker_RaindropRockAirWater_FACTOR to 10.
1017 for example.
1018 
1019 **/
1020 
1021 void U4VolumeMaker::RaindropRockAirWater_Configure(
1022     std::vector<std::string>& mats,
1023     std::vector<double>& rindex,
1024     double& universe_halfside,    // VACUUM
1025     double& container_halfside,   // Rock
1026     double& medium_halfside,      // Air
1027     double& drop_radius,          // Water
1028     std::string& dropshape        // Orb
1029     )
1030 {
1031     ssys::fill_evec(mats,   U4VolumeMaker_RaindropRockAirWater_MATS, "VACUUM,G4_Pb,G4_AIR,G4_WATER", ',' ) ;
1032     ssys::fill_evec(rindex, U4VolumeMaker_RaindropRockAirWater_RINDEX, "0,0,1,1.333", ',' ) ;
1033 
1034     double halfside = ssys::getenv_<double>(U4VolumeMaker_RaindropRockAirWater_HALFSIDE, 100.);
1035     double factor   = ssys::getenv_<double>(U4VolumeMaker_RaindropRockAirWater_FACTOR,   1.);
1036     dropshape = ssys::getenvvar(U4VolumeMaker_RaindropRockAirWater_DROPSHAPE, "Orb");
1037 
1038     LOG(LEVEL) << U4VolumeMaker_RaindropRockAirWater_MATS << " " << ssys::desc_vec(&mats) ;
1039     LOG(LEVEL) << U4VolumeMaker_RaindropRockAirWater_HALFSIDE << " " << halfside ;
1040     LOG(LEVEL) << U4VolumeMaker_RaindropRockAirWater_FACTOR   << " " << factor ;
1041     LOG(LEVEL) << U4VolumeMaker_RaindropRockAirWater_DROPSHAPE << " " << dropshape ;
1042 
1043     universe_halfside = 1.2*halfside*factor ;    // VACUUM
1044     container_halfside = 1.1*halfside*factor ;   // rock  : G4_Pb
1045     medium_halfside = halfside*factor ;          // air   : G4_AIR
1046     drop_radius = halfside/2. ;                  // water : G4_WATER
1047 }
1048 
1049 /**
1050 U4VolumeMaker::RaindropRockAirWater
1051 ------------------------------------
1052 
1053 +-----------+--------------+
1054 | original  |  generalized |
1055 +===========+==============+
1056 | Rock      | container    |
1057 +-----------+--------------+
1058 | Air       | medium       |
1059 +-----------+--------------+
1060 | Water     | drop         |
1061 +-----------+--------------+
1062 
1063 
1064 Notice that so long as all the LV are created prior to creating the PV,
1065 which need the LV for placement and mother logical, there is no need to be
1066 careful with creation order of the volumes.
1067 
1068 This is suggestive of how to organize the API, instead of focussing on methods
1069 to create PV it is more flexible to have API that create LV that are then put
1070 together by the higher level methods that make less sense to generalize.
1071 
1072 **/
1073 
1074 const G4VPhysicalVolume* U4VolumeMaker::RaindropRockAirWater(bool sd)
1075 {
1076     std::vector<std::string> mats ;
1077     std::vector<double> rindex ;
1078 
1079     double universe_halfside ;
1080     double container_halfside ;
1081     double medium_halfside ;
1082     double drop_radius ;
1083     std::string dropshape ;
1084 
1085     RaindropRockAirWater_Configure( mats, rindex, universe_halfside, container_halfside, medium_halfside, drop_radius, dropshape );
1086 
1087     static const int N = 4 ;
1088     bool mats_expect = mats.size() == N  ;
1089     bool rindex_expect = rindex.size() == N ;
1090 
1091     LOG_IF(error, !mats_expect)   << " unexpected mats.size " << mats.size() << " expecting " << N ;
1092     LOG_IF(error, !rindex_expect) << " unexpected rindex.size " << rindex.size() << " expecting " << N ;
1093     assert( mats_expect );
1094     assert( rindex_expect );
1095 
1096     std::array<G4Material*, N> materials ;
1097 
1098     for(int i=0 ; i < N ; i++)
1099     {
1100         const char* mat = mats[i].c_str() ;
1101         G4Material* material = U4Material::Get(mat) ;
1102         if( material == nullptr ) std::cerr
1103             << "U4Material::Get"
1104             << " FAILED FOR [" << mat << "]"
1105             << std::endl
1106             ;
1107         materials[i] = material ;
1108     }
1109 
1110     G4Material* universe_material = materials[0] ;  // default "VACUUM"
1111     G4Material* container_material = materials[1] ;  // default "G4_Pb"
1112     G4Material* medium_material  = materials[2] ;    // default "G4_AIR"
1113     G4Material* drop_material  = materials[3] ;      // default "G4_WATER"
1114 
1115     // fix up materials that dont have RINDEX properties, if that is configured
1116     for(int i=0 ; i < N ; i++)
1117     {
1118         G4Material* mat = materials[i] ;
1119         assert(mat);
1120 
1121         if(rindex[i] > 0. && !U4Material::HasMPT(mat) )
1122         {
1123             U4Material::SetMPT(mat,U4MaterialPropertiesTable::Create("RINDEX", rindex[i])) ;
1124         }
1125     }
1126 
1127 
1128     G4LogicalVolume* universe_lv  = Box_(universe_halfside, universe_material );
1129     G4LogicalVolume* container_lv  = Box_(container_halfside, container_material );
1130     G4LogicalVolume* medium_lv   = Box_(medium_halfside, medium_material );
1131 
1132     G4LogicalVolume* drop_lv = nullptr ;
1133     if(     strcmp(dropshape.c_str(), "Orb") == 0 ) drop_lv = Orb_(drop_radius, drop_material );
1134     else if(strcmp(dropshape.c_str(), "Box") == 0 ) drop_lv = Box_(drop_radius, drop_material );
1135     LOG_IF(fatal, drop_lv == nullptr ) << " drop_lv NULL : dropshape must be one of Orb,Box not [" << dropshape  << "]" ;
1136     assert( drop_lv );
1137 
1138     const G4VPhysicalVolume* drop_pv = new G4PVPlacement(0,G4ThreeVector(), drop_lv ,"drop_pv", medium_lv,false,0);
1139     const G4VPhysicalVolume* medium_pv = new G4PVPlacement(0,G4ThreeVector(),   medium_lv ,  "medium_pv",  container_lv,false,0);
1140     const G4VPhysicalVolume* container_pv = new G4PVPlacement(0,G4ThreeVector(),  container_lv ,  "container_pv", universe_lv,false,0);
1141     const G4VPhysicalVolume* universe_pv = new G4PVPlacement(0,G4ThreeVector(),  universe_lv ,  "universe_pv", nullptr,false,0);
1142 
1143     assert( drop_pv );
1144     assert( medium_pv );
1145     assert( container_pv );
1146     assert( universe_pv );
1147 
1148     G4LogicalBorderSurface* medium_container_bs = U4Surface::MakePerfectAbsorberBorderSurface(
1149        "medium_container_bs", medium_pv, container_pv );
1150     assert( medium_container_bs );
1151     if(!medium_container_bs) std::raise(SIGINT);
1152 
1153     if(sd)
1154     {
1155         //drop_lv->SetSensitiveDetector(new U4SensitiveDetector("drop_sd"));
1156         // not needed, it is the below surface that is necessary
1157 
1158         G4LogicalBorderSurface* medium_drop_bs = U4Surface::MakePerfectDetectorBorderSurface("medium_drop_bs", medium_pv, drop_pv );
1159         assert( medium_drop_bs );
1160         if(!medium_drop_bs) std::raise(SIGINT);
1161 
1162     }
1163     return universe_pv ;
1164 }
1165 
1166 
1167 
1168 const G4VPhysicalVolume* U4VolumeMaker::BigWaterPool()
1169 {
1170 
1171     double m_heightWP = 43500 ;
1172     double m_deadWaterThickness = 100 ;
1173     double m_tyvekThickness = 2 ;
1174     double m_radWP = 21750 ;
1175 
1176     G4Material* vetoWater = U4Material::Get("G4_WATER") ;
1177     G4VSolid* solidWaterPool = nullptr ;
1178     G4LogicalVolume* logicWaterPool = nullptr ;
1179 
1180     const G4VPhysicalVolume* pTyvekFilm = nullptr ;
1181     const G4VPhysicalVolume* pWaterPool = nullptr ;
1182 
1183     G4VSolid* solidTyvekFilm = nullptr ;
1184     G4LogicalVolume* logicAirGap = nullptr ;
1185     G4LogicalVolume* logicDeadWater = nullptr ;
1186     G4LogicalVolume* logicTyvekFilm = nullptr ;
1187 
1188     G4Material* DeadWater = U4Material::Get("G4_WATER") ;
1189     G4Material* Tyvek = U4Material::Get("G4_WATER") ;
1190     G4Material* BufferMaterials = U4Material::Get("G4_WATER") ;
1191 
1192     std::array<G4Material*,3> mats = {DeadWater, Tyvek, BufferMaterials };
1193     // ALL THOSE ARE SAME POINTER
1194     for(int i=0 ; i < int(mats.size()) ; i++)
1195     {
1196         G4Material* mat = mats[i];
1197         std::cout
1198              << " i " << i
1199              << " mat " << std::hex << size_t(mat)  << std::dec
1200              << " mpt " << U4Material::HasMPT(mat)
1201              << "\n"
1202              ;
1203 
1204         if(!U4Material::HasMPT(mat)) U4Material::SetMPT(mat,U4MaterialPropertiesTable::Create("RINDEX", 1.333)) ;
1205     }
1206 
1207 
1208 
1209    // DeadWater:logicDeadWater
1210     {
1211         G4VSolid* solidDeadWater = new G4Tubs("sDeadWater_shell", 0., m_radWP, m_heightWP/2., 0.*deg, 360*deg);
1212         logicDeadWater = new G4LogicalVolume(solidDeadWater, DeadWater, "lDeadWater", 0, 0, 0);
1213     }
1214 
1215 
1216    // Tyvek:logicTyvekFilm
1217     {
1218         double tyvek_btmZ = -m_heightWP/2 + m_deadWaterThickness - m_tyvekThickness;
1219         double tyvek_topZ = -tyvek_btmZ;
1220         double tyvek_outR = m_radWP - m_deadWaterThickness + m_tyvekThickness;
1221         //double tyvek_inR  = m_radWP - m_deadWaterThickness;
1222 
1223         static G4double  zPlane[] = {tyvek_btmZ, tyvek_topZ}; //-21.65m-2mm, -21.65m, 21.65m, 21.65m+2mm
1224         static G4double  rInner[] = {0., 0.};
1225         static G4double  rOuter[] = {tyvek_outR, tyvek_outR}; //21.65m+2mm
1226 
1227         G4VSolid* solidTyvek_shell = new G4Polycone("sTyvek_shell", //const G4String& pName,
1228                                          0, //G4double  phiStart,
1229                                          360*deg, //G4double  phiTotal,
1230                                          2, //G4int         numZPlanes,
1231                                          zPlane, //const G4double  zPlane[],
1232                                          rInner, //const G4double  rInner[],
1233                                          rOuter //const G4double  rOuter[])
1234                                         );
1235 
1236 
1237          solidTyvekFilm = solidTyvek_shell ;
1238          logicTyvekFilm = new G4LogicalVolume(solidTyvekFilm, Tyvek, "lTyvekFilm", 0, 0, 0);
1239      }
1240 
1241 
1242 
1243     // vetoWater:logicWaterPool
1244      {
1245         double vetoWater_btmZ = -m_heightWP/2 + m_deadWaterThickness;  // 2025/7/2 added m_deadWaterThickness
1246         double vetoWater_topZ =  m_heightWP/2 - m_deadWaterThickness;
1247         double vetoWater_outR = m_radWP - m_deadWaterThickness ;
1248 
1249         static G4double  zPlane[] = {vetoWater_btmZ, vetoWater_topZ};
1250         static G4double  rInner[] = {0., 0.};
1251         static G4double  rOuter[] = {vetoWater_outR, vetoWater_outR};
1252         solidWaterPool = new G4Polycone("sOuterWaterPool", //const G4String& pName,
1253                                          0, //G4double  phiStart,
1254                                          360*deg, //G4double  phiTotal,
1255                                          2, //G4int         numZPlanes,
1256                                          zPlane, //const G4double  zPlane[],
1257                                          rInner, //const G4double  rInner[],
1258                                          rOuter //const G4double  rOuter[])
1259                                         );
1260 
1261         logicWaterPool = new G4LogicalVolume(solidWaterPool,
1262             vetoWater,
1263             "lOuterWaterPool",
1264             0, 0, 0);
1265      }
1266 
1267 
1268 
1269 
1270     G4LogicalVolume* logicMaskVirtual = nullptr ;
1271     {
1272         // jcv R12860OnlyFrontMaskManager
1273 
1274         G4double htop_thickness=8.*mm ;
1275         G4double htop_gap=2.*mm ;
1276 
1277         //G4double MAGIC_virtual_thickness = 0.05*mm;
1278         G4double MAGIC_virtual_thickness = 0.10*mm;
1279 
1280         G4double requator_thickness=8.*mm ;
1281         G4double requator_gap=2.*mm ;
1282 
1283         G4double mask_radiu_in = 254.*mm + requator_gap;
1284         G4double mask_radiu_out = mask_radiu_in + requator_thickness;
1285         G4double mask_radiu_virtual = mask_radiu_out + MAGIC_virtual_thickness;
1286 
1287 
1288         G4double htop_in = 184.*mm + htop_gap;
1289         G4double htop_out = htop_in + htop_thickness;
1290 
1291 
1292 
1293         G4double height_in = 10*mm + 5.*mm;
1294         G4double height_out = height_in + 10*mm;
1295         G4double height_virtual = height_out + MAGIC_virtual_thickness;
1296 
1297         G4double zPlane[] = {
1298                             -height_virtual,
1299                             htop_out + MAGIC_virtual_thickness
1300                             };
1301         G4double rInner[] = {0.,
1302                              0.};
1303         G4double rOuter[] = {mask_radiu_virtual,
1304                              mask_radiu_virtual};
1305 
1306         G4VSolid* SolidMaskVirtual = new G4Polycone(
1307                                     "sMask_virtual",
1308                                     0,
1309                                     360*deg,
1310                                     2,
1311                                     zPlane,
1312                                     rInner,
1313                                     rOuter
1314                                     );
1315 
1316         logicMaskVirtual = new G4LogicalVolume(
1317             SolidMaskVirtual,
1318             BufferMaterials,
1319             "lMaskVirtual",
1320             0,
1321             0,
1322             0);
1323     }
1324 
1325 
1326     // place logicDeadWater within logicAirGap
1327     new G4PVPlacement(0, G4ThreeVector(0, 0, 0),logicDeadWater,"pDeadWater",logicAirGap , false, 0);
1328 
1329     // place logicTyvekFilm within logicDeadWater
1330     pTyvekFilm = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),logicTyvekFilm,"pTyvekFilm",logicDeadWater , false, 0);
1331 
1332     // place logicWaterPool within logicTyvekFilm
1333     pWaterPool = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),logicWaterPool,"pWaterPool",logicTyvekFilm , false, 0);
1334 
1335     U4Surface::MakePerfectAbsorberBorderSurface( "WaterPool_Tyvek_bs", pWaterPool, pTyvekFilm );
1336 
1337 
1338     {
1339         // grabbed from gdb VetoPmtPosBall::initialize    "b 82 if i == 51006"
1340 
1341         G4double x = -994.96197509765625 ;
1342         G4double y = 20447.890625 ;
1343         G4double z = 730.7440185546875 ;
1344         G4double theta = 1.53511662609056 ;
1345         G4double phi = 1.6194159277617648 ;
1346 
1347         G4ThreeVector pos(x, y, z);
1348         G4RotationMatrix rot;
1349         rot.rotateY(theta);
1350         rot.rotateZ(phi);
1351         G4Transform3D trans(rot, pos);
1352         G4Transform3D& tr = trans ;
1353 
1354         // WaterPoolConstruction::inject
1355         auto _rot = tr.getRotation();
1356         auto _translation = tr.getTranslation();
1357         _translation.setR(20.502*m);
1358 
1359         G4Transform3D tr_ideal(_rot, _translation);
1360 
1361         // place logicMaskVirtual within logicWaterPool
1362         new G4PVPlacement(tr_ideal,logicMaskVirtual,"pMask",logicWaterPool , false, 0);
1363     }
1364 
1365     return pTyvekFilm ;
1366 
1367 }
1368 
1369 
1370 
1371 G4Material* U4VolumeMaker::GetMaterial(const char* mat)
1372 {
1373     G4Material* material  = U4Material::Get(mat);
1374     LOG_IF(fatal, material==nullptr) << " U4Material::Get failed for mat[" << ( mat ? mat : "-" ) << "]" ;
1375     assert(material);
1376     return material ;
1377 }
1378 
1379 G4LogicalVolume* U4VolumeMaker::Orb_( double radius, const char* mat, const char* prefix )
1380 {
1381     if( prefix == nullptr ) prefix = mat ;
1382     G4Material* material = GetMaterial(mat) ;
1383     return Orb_(radius, material, prefix );
1384 }
1385 G4LogicalVolume* U4VolumeMaker::Orb_( double radius, G4Material* material, const char* prefix )
1386 {
1387     if( prefix == nullptr ) prefix = material->GetName().c_str() ;
1388     G4Orb* solid = new G4Orb( spath::Name(prefix,"_solid"), radius );
1389     G4LogicalVolume* lv = new G4LogicalVolume( solid, material, spath::Name(prefix,"_lv"));
1390     return lv ;
1391 }
1392 G4LogicalVolume* U4VolumeMaker::Box_( double halfside, const char* mat, const char* prefix, const double* scale )
1393 {
1394     if( prefix == nullptr ) prefix = mat ;
1395     G4Material* material = GetMaterial(mat) ;
1396     return Box_( halfside, material, prefix, scale );
1397 }
1398 
1399 G4LogicalVolume* U4VolumeMaker::Box_( double halfside, G4Material* material, const char* prefix, const double* scale )
1400 {
1401     if( prefix == nullptr ) prefix = material->GetName().c_str() ;
1402     double hx = scale ? scale[0]*halfside : halfside ;
1403     double hy = scale ? scale[1]*halfside : halfside ;
1404     double hz = scale ? scale[2]*halfside : halfside ;
1405 
1406     LOG(LEVEL) << " hx " << hx << " hy " << hy << " hz " << hz << " halfside " << halfside ;
1407 
1408     G4Box* solid = new G4Box( spath::Name(prefix,"_solid"), hx, hy, hz );
1409     G4LogicalVolume* lv = new G4LogicalVolume( solid, material, spath::Name(prefix,"_lv"));
1410     return lv ;
1411 }
1412 
1413 
1414 NP* U4VolumeMaker::MakeTransforms( const char* name, const char* prefix )
1415 {
1416     const char* opts = "TR,tr,R,T,r,t" ;
1417     NP* trs = nullptr ;
1418 
1419     if(strcmp(prefix, "AroundSphere")==0)
1420     {
1421         double radius = 17000. ;
1422         double item_arclen = 600. ;   // 400. has lots of overlap, 1000. too spaced
1423         unsigned num_ring = 8 ;
1424         trs = SPlace::AroundSphere(  opts, radius, item_arclen, num_ring );
1425     }
1426     else if(strcmp(prefix, "AroundCylinder")==0)
1427     {
1428         double radius = 17000. ;
1429         double halfheight = radius ;
1430         unsigned num_ring = 8 ;
1431         unsigned num_in_ring = 16 ;
1432         trs = SPlace::AroundCylinder(opts, radius, halfheight, num_ring, num_in_ring  );
1433     }
1434     else if(strcmp(prefix, "AroundCircle")==0)
1435     {
1436         double radius = ssys::getenv_<double>(U4VolumeMaker_MakeTransforms_AroundCircle_radius, 1000.);
1437         unsigned numInRing = ssys::getenv_<unsigned>(U4VolumeMaker_MakeTransforms_AroundCircle_numInRing, 4u );
1438         double fracPhase = ssys::getenv_<double>(U4VolumeMaker_MakeTransforms_AroundCircle_fracPhase, 0.);
1439         trs = SPlace::AroundCircle(opts, radius, numInRing, fracPhase  );
1440     }
1441     return trs ;
1442 }
1443 
1444 /**
1445 U4VolumeMaker::WrapAround
1446 ----------------------------
1447 
1448 Place *lv* multiple times inside *mother_lv* using the *trs*
1449 tranforms to configure the placement rotations and translations.
1450 
1451 NB there is no check that the placements are within the mother_lv
1452 so it is up to the user to ensure that the mother volume is sufficiently
1453 large to accommodate the lv with all the transforms applied to it.
1454 
1455 **/
1456 
1457 void U4VolumeMaker::WrapAround( const char* prefix, const NP* trs, std::vector<G4LogicalVolume*>& lvs, G4LogicalVolume* mother_lv )
1458 {
1459     unsigned num_place = trs->shape[0] ;
1460     unsigned place_tr = trs->shape[1] ;
1461     unsigned place_values = place_tr*4*4 ;
1462     unsigned num_lv = lvs.size();
1463 
1464 
1465     assert( trs->has_shape(num_place,place_tr,4,4) );
1466     assert( place_tr == 6 );  // expected number of different options from "TR,tr,R,T,r,t"
1467     enum { _TR, _tr, _R, _T, _r, _t } ;  // order must match opts "TR,tr,R,T,r,t"
1468 
1469     const double* tt = trs->cvalues<double>();
1470 
1471     for(unsigned i=0 ; i < num_place ; i++)
1472     {
1473         const double* T = tt + place_values*i + _T*16 ;
1474         const double* R = tt + place_values*i + _R*16 ;
1475 
1476         // TODO: get these from a single matrix, not 6
1477 
1478         G4ThreeVector tla( T[12], T[13], T[14] );
1479 
1480         LOG(LEVEL) << " i " << std::setw(7) << " tla " << U4ThreeVector::Desc(tla) ;
1481 
1482         bool transpose = true ;
1483         U4RotationMatrix* rot = new U4RotationMatrix( R, transpose );  // ISA G4RotationMatrix
1484 
1485         LOG(LEVEL) << " i " << std::setw(7) << " rot " << U4RotationMatrix::Desc(rot) ;
1486 
1487         const char* iname = PlaceName(prefix, i, nullptr);
1488 
1489         G4bool pMany_unused = false ;
1490         G4int  pCopyNo = (i+1)*10 ;
1491         G4LogicalVolume* lv = lvs[i%num_lv] ;
1492 
1493         const G4VPhysicalVolume* pv_n = new G4PVPlacement(rot, tla, lv, iname, mother_lv, pMany_unused, pCopyNo ); // 1st ctor
1494         assert( pv_n );
1495         if(!pv_n) std::raise(SIGINT);
1496 
1497     }
1498 }
1499
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!