eic-opticks/u4/U4Material.cc

0001
0002 #include <string>
0003 #include <vector>
0004 #include <sstream>
0005 #include <limits>
0006 #include <csignal>
0007
0008 #include "ssys.h"
0009 #include "spath.h"
0010
0011
0012 #include "SDir.h"
0013 #include "SStr.hh"
0014 #include "sstr.h"
0015 #include "SSim.hh"
0016 #include "SBnd.h"
0017 #include "sdomain.h"
0018 #include "sproplist.h"
0019
0020 #include "NPFold.h"
0021 #include "NP.hh"
0022 #include "SLOG.hh"
0023
0024 #include "G4Material.hh"
0025
0026 #include "G4SystemOfUnits.hh"
0027 #include "G4Version.hh"
0028 #include "U4Material.hh"
0029 #include "U4MaterialPropertyVector.h"
0030 #include "U4NistManager.h"
0031
0032 const plog::Severity U4Material::LEVEL = SLOG::EnvLevel("U4Material", "DEBUG");
0033
0034
0035 std::string U4Material::DescMaterialTable()
0036 {
0037     std::vector<std::string> names ;
0038     GetMaterialNames(names);
0039     std::stringstream ss ;
0040     ss << "U4Material::DescMaterialTable " << names.size() << std::endl ;
0041     for(int i=0 ; i < int(names.size()) ; i++) ss << names[i] << std::endl ;
0042     std::string s = ss.str();
0043     return s ;
0044 }
0045
0046 void U4Material::GetMaterialNames( std::vector<std::string>& names, bool extras )
0047 {
0048     G4MaterialTable* tab =  G4Material::GetMaterialTable();
0049     for(int i=0 ; i < int(tab->size()) ; i++)
0050     {
0051         G4Material* mat = (*tab)[i] ;
0052         const G4String& name = mat->GetName();
0053         names.push_back(name);
0054     }
0055     if(extras)
0056     {
0057         names.push_back(U4Material::SCINTILLATOR);
0058         names.push_back(U4Material::VACUUM);
0059    }
0060 }
0061
0062 /**
0063 U4Material::FindMaterialName
0064 ------------------------------
0065
0066 Look for the names of any of the materials within the volname
0067 string and return the first material name found.
0068 If more than one names are matched this asserts : so be careful
0069 to avoid unwanted material names within volume names.
0070
0071 **/
0072
0073 const char* U4Material::FindMaterialName(const char* volname)
0074 {
0075     std::vector<std::string> names ;
0076     GetMaterialNames(names, true);
0077     const char* matname = nullptr ;
0078     unsigned count = 0 ;
0079     for(unsigned i=0 ; i < names.size() ; i++)
0080     {
0081         const char* mat = names[i].c_str();
0082         if(strstr(volname, mat) != nullptr )
0083         {
0084             if(count == 0 ) matname = strdup(mat);
0085             count+= 1 ;
0086         }
0087     }
0088     LOG_IF(fatal, count > 1) << " count " << count << " volname " << volname << " matname " << matname ;
0089     assert( count < 2 );
0090     return matname ;
0091 }
0092
0093
0094
0095
0096 /**
0097 U4Material::Get
0098 -----------------
0099
0100 Material resolution order:
0101
0102 1. ordinary G4Material::GetMaterial
0103 2. U4Material::Get_ which implements a few specials : Scintillator and Vacuum
0104 3. U4NistManager::GetMaterial for names starting with "G4_"
0105
0106 **/
0107
0108 G4Material* U4Material::Get(const char* name)
0109 {
0110    G4Material* material = G4Material::GetMaterial(name, false);
0111    if( material == nullptr ) material = Get_(name);
0112    if(sstr::StartsWith(name, "G4_")) material = U4NistManager::GetMaterial(name) ;
0113    return material ;
0114 }
0115
0116 G4Material* U4Material::Get_(const char* name)
0117 {
0118     G4Material* material = nullptr ;
0119     if(strcmp(name, SCINTILLATOR)==0) material = MakeScintillator();
0120     if(strcmp(name, VACUUM)==0)       material = Vacuum(name);
0121     return material ;
0122 }
0123
0124 G4Material* U4Material::Vacuum(const char* name)
0125 {
0126     bool dump = false ;
0127     if(dump) std::cout << "[ U4Material::Vacuum " << std::endl ;
0128     G4double z = 1. ;
0129     G4double a = 1.01*CLHEP::g/CLHEP::mole ;
0130     G4double density = 1.00001*CLHEP::universe_mean_density ;  // curious setting to 1. gives a warning
0131     if(dump) std::cout << " U4Material::Vacuum density " << std::scientific << density << std::endl ;
0132     G4Material* material = new G4Material(name, z, a, density );
0133     if(dump) std::cout << "] U4Material::Vacuum " << std::endl ;
0134     return material ;
0135 }
0136
0137 G4MaterialPropertyVector* U4Material::GetProperty(const G4Material* mat, const char* pname)
0138 {
0139     G4MaterialPropertiesTable* mpt = mat->GetMaterialPropertiesTable();
0140     return mpt->GetProperty(pname) ;
0141 }
0142
0143 void U4Material::GetMinMaxValue( double& mn , double& mx, const G4MaterialPropertyVector* prop )
0144 {
0145     size_t len = prop->GetVectorLength() ;
0146     double MAX = std::numeric_limits<double>::max();
0147     mn = MAX ;
0148     mx = -MAX ;
0149     for(unsigned i=0 ; i < len ; i++)
0150     {
0151         double v = (*prop)[i] ;
0152         if( mx < v ) mx = v ;
0153         if( mn > v ) mn = v ;
0154     }
0155 }
0156
0157 std::string U4Material::DescProperty(const G4MaterialPropertyVector* prop)
0158 {
0159     std::stringstream ss ;
0160     if( prop == nullptr )
0161     {
0162         ss << " prop null " ;
0163     }
0164     else
0165     {
0166         G4MaterialPropertyVector* prop_ = const_cast<G4MaterialPropertyVector*>(prop);
0167         double mn, mx ;
0168         GetMinMaxValue(mn, mx, prop);
0169         ss
0170            << std::setw(16) << "GetVectorLength" << std::setw(5) << prop->GetVectorLength()
0171            << std::setw(5) << "mn" << std::fixed << std::setw(10) << std::setprecision(5) << mn
0172            << std::setw(5) << "mx" << std::fixed << std::setw(10) << std::setprecision(5) << mx
0173            << std::setw(12) << "GetMaxValue" << std::fixed << std::setw(10) << std::setprecision(5) << prop_->GetMaxValue()
0174            << std::setw(12) << "GetMinValue" << std::fixed << std::setw(10) << std::setprecision(5) << prop_->GetMinValue()
0175 #if G4VERSION_NUMBER < 1100
0176            << std::setw(23) << "GetMinLowEdgeEnergy/eV" << std::fixed << std::setw(10) << std::setprecision(5) << prop_->GetMinLowEdgeEnergy()/eV
0177            << std::setw(23) << "GetMaxLowEdgeEnergy/eV" << std::fixed << std::setw(10) << std::setprecision(5) << prop_->GetMaxLowEdgeEnergy()/eV
0178 #else
0179            << std::setw(23) << "GetMinLowEdgeEnergy/eV" << std::fixed << std::setw(10) << std::setprecision(5) << prop_->GetMinEnergy()/eV
0180            << std::setw(23) << "GetMaxLowEdgeEnergy/eV" << std::fixed << std::setw(10) << std::setprecision(5) << prop_->GetMaxEnergy()/eV
0181 #endif
0182       ;
0183
0184     }
0185     std::string s = ss.str();
0186     return s ;
0187 }
0188
0189 std::string U4Material::DescProperty(const char* mname, const char* pname)
0190 {
0191     const G4MaterialPropertyVector* prop = GetProperty(mname, pname);
0192     std::stringstream ss ;
0193     ss << std::setw(20) << mname << std::setw(20) << pname ;
0194     ss << DescProperty(prop) ;
0195     std::string s = ss.str();
0196     return s ;
0197 }
0198
0199 std::string U4Material::DescProperty(const char* mname)
0200 {
0201     G4Material* mat = G4Material::GetMaterial(mname, false);
0202     std::vector<std::string> pnames ;
0203     GetPropertyNames(pnames, mat);
0204
0205     std::stringstream ss ;
0206     for(unsigned i=0 ; i < pnames.size() ; i++)  ss << DescProperty(mname, pnames[i].c_str() ) << std::endl ;
0207     std::string s = ss.str();
0208     return s ;
0209 }
0210
0211 std::string U4Material::DescProperty()
0212 {
0213     std::vector<std::string> mnames ;
0214     GetMaterialNames(mnames);
0215     std::stringstream ss ;
0216     for(unsigned i=0 ; i < mnames.size() ; i++)  ss << DescProperty(mnames[i].c_str()) << std::endl ;
0217     std::string s = ss.str();
0218     return s ;
0219 }
0220
0221
0222 G4MaterialPropertyVector* U4Material::GetProperty(const char* mname, const char* pname)
0223 {
0224     bool warning = false ;
0225     const G4Material* mat = G4Material::GetMaterial(mname, warning);
0226     return mat ? GetProperty(mat, pname) : nullptr ;
0227 }
0228
0229
0230 void U4Material::RemoveProperty( const char* key, G4Material* mat )
0231 {
0232     if(mat == nullptr) return ;
0233     const G4String& name = mat->GetName();
0234     G4MaterialPropertiesTable* mpt = mat->GetMaterialPropertiesTable();
0235     G4MaterialPropertyVector* prop = mpt->GetProperty(key);
0236
0237     if( prop == nullptr)
0238     {
0239          LOG(info) << " material " << name << " does not have property " << key ;
0240     }
0241     else
0242     {
0243          LOG(info) << " material " << name << " removing property " << key ;
0244          mpt->RemoveProperty(key);
0245          prop = mpt->GetProperty(key);
0246          assert( prop == nullptr );
0247     }
0248 }
0249
0250 /**
0251 U4Material::GetPropertyNames
0252 -------------------------------
0253
0254 Note difference between Geant4 versions.
0255 With older Geant4 only existing property names are returned.
0256 With newer Geant4 every property name under the sun is returned.
0257
0258 **/
0259
0260
0261 #if G4VERSION_NUMBER < 1100
0262 void U4Material::GetPropertyNames( std::vector<std::string>& names, const G4Material* mat )
0263 {
0264     G4MaterialPropertiesTable* mpt = mat->GetMaterialPropertiesTable();
0265     const G4MaterialPropertiesTable* mpt_ = mat->GetMaterialPropertiesTable();
0266
0267     if( mpt_ == nullptr ) return ;
0268
0269     std::vector<G4String> pnames = mpt_->GetMaterialPropertyNames();
0270
0271     typedef std::map<G4int, G4MaterialPropertyVector*, std::less<G4int> > MIV ;
0272     const MIV* miv =  mpt->GetPropertyMap();
0273     for(MIV::const_iterator it=miv->begin() ; it != miv->end() ; it++ )
0274     {
0275          G4String name = pnames[it->first] ;
0276          names.push_back(name.c_str()) ;
0277     }
0278 }
0279 #else
0280 void U4Material::GetPropertyNames( std::vector<std::string>& names, const G4Material* mat )
0281 {
0282   const G4MaterialPropertiesTable* mpt_ = mat->GetMaterialPropertiesTable();
0283   if( mpt_ == nullptr ) return ;
0284   std::vector<G4String> pnames = mpt_->GetMaterialPropertyNames();
0285   for(std::vector<G4String>::const_iterator iter = pnames.begin(); iter != pnames.end(); iter++)
0286   {
0287       names.push_back(*iter);  // from Ilker Parmaksiz testing with Geant4 11.1.1
0288   }
0289 }
0290 #endif
0291 /**
0292 U4Material::MakePropertyFold_flat
0293 ------------------------------------
0294
0295 Converts the properties of all materials into a
0296 NPFold with keys of form : "materialName"/"propName"
0297
0298
0299 Note that U4Material::MakePropertyFold uses
0300 subfold feature additions to NPFold.h
0301 which more faithfully maps the file system
0302 model into an in memory model. For example
0303 it is then natural to load single materials or
0304 single properties.
0305
0306 **/
0307
0308 NPFold* U4Material::MakePropertyFold_flat()
0309 {
0310     NPFold* fold = new NPFold ;
0311
0312     std::vector<std::string> matnames ;
0313     GetMaterialNames(matnames);
0314     for(unsigned i=0 ; i < matnames.size() ; i++)
0315     {
0316         const char* material = matnames[i].c_str();
0317         const G4Material* mat = G4Material::GetMaterial(material) ;
0318
0319         std::vector<std::string> propnames ;
0320         GetPropertyNames( propnames, mat );
0321
0322         for(unsigned j=0 ; j < propnames.size() ; j++)
0323         {
0324             const char* propname = propnames[j].c_str() ;
0325             G4MaterialPropertyVector* prop = GetProperty(mat, propname );
0326             NP* a = U4MaterialPropertyVector::ConvertToArray(prop);
0327             fold->add( SStr::Format("%s/%s", material, propname), a );
0328         }
0329     }
0330     return fold ;
0331 }
0332
0333 /**
0334 U4Material::MakePropertyFold
0335 ------------------------------
0336
0337 Canonically invoked from U4Tree::convertMaterials
0338 with result going into st.material
0339
0340 Note that NPFold::set_allowempty_r is invoked
0341 after populating the material fold.
0342 This allows the material fold and it subfold recursively
0343 to be empty without supression. This is a workaround
0344 to avoid suppressing some converted Geant4 materials
0345 that lack any properties (for example JUNO Galactic).
0346
0347 **/
0348
0349 NPFold* U4Material::MakePropertyFold()
0350 {
0351     NPFold* fold = new NPFold ;
0352
0353     std::vector<std::string> matnames ;
0354     GetMaterialNames(matnames);
0355     for(unsigned i=0 ; i < matnames.size() ; i++)
0356     {
0357         const char* matname = matnames[i].c_str();
0358         const G4Material* mat = G4Material::GetMaterial(matname) ;
0359         NPFold* matfold = MakePropertyFold(mat) ;
0360         fold->add_subfold( matname, matfold );
0361     }
0362     fold->set_allowempty_r(true);
0363     return fold ;
0364 }
0365
0366
0367 NPFold* U4Material::MakePropertyFold(std::vector<const G4Material*>& mats)
0368 {
0369     NPFold* fold = new NPFold ;
0370     for(unsigned i=0 ; i < mats.size() ; i++)
0371     {
0372         const G4Material* mt = mats[i]  ;
0373         const G4String& mtname = mt->GetName() ;
0374         const char* mtn = mtname.c_str();
0375         NPFold* matfold = MakePropertyFold(mt) ;
0376         fold->add_subfold( mtn, matfold );
0377     }
0378     return fold ;
0379 }
0380
0381
0382
0383
0384
0385 NPFold* U4Material::MakePropertyFold(const G4Material* mat )
0386 {
0387     NPFold* fold = new NPFold ;
0388     std::vector<std::string> propnames ;
0389     GetPropertyNames( propnames, mat );
0390     for(unsigned i=0 ; i < propnames.size() ; i++)
0391     {
0392         const char* propname = propnames[i].c_str() ;
0393         G4MaterialPropertyVector* prop = GetProperty(mat, propname );
0394         if(prop == nullptr) continue ;  // from Ilker Parmaksiz testing with Geant4 11.1.1
0395         NP* a = U4MaterialPropertyVector::ConvertToArray(prop);
0396         fold->add( propname, a );
0397     }
0398     return fold ;
0399 }
0400
0401
0402
0403
0404
0405 bool U4Material::HasMPT(const G4Material* mat) // static
0406 {
0407    return nullptr != mat->GetMaterialPropertiesTable() ;
0408 }
0409
0410 void U4Material::SetMPT(G4Material* mat, G4MaterialPropertiesTable* mpt)
0411 {
0412     mat->SetMaterialPropertiesTable(mpt);
0413 }
0414
0415
0416 std::string U4Material::DescPropertyNames( const G4Material* mat )
0417 {
0418     const G4MaterialPropertiesTable* mpt_ = mat->GetMaterialPropertiesTable();
0419     std::vector<G4String> pnames = mpt_->GetMaterialPropertyNames();
0420     std::vector<G4String> cnames = mpt_->GetMaterialConstPropertyNames();
0421     std::stringstream ss ;
0422     ss << "U4Material::DescPropertyNames " << mat->GetName() << std::endl ;
0423
0424     /*
0425     ss << " MaterialPropertyNames.size " << pnames.size() << std::endl ;
0426     for(int i=0 ; i < int(pnames.size()) ; i++) ss << pnames[i] << std::endl ;
0427
0428     ss << " MaterialConstPropertyNames.size " << cnames.size() << std::endl ;
0429     for(int i=0 ; i < int(cnames.size()) ; i++) ss << cnames[i] << std::endl ;
0430     */
0431
0432     /*
0433     ss << " GetPropertiesMap " << std::endl ;
0434     typedef std::map< G4String, G4MaterialPropertyVector*, std::less<G4String> > MSV ;
0435     MSV* msv = mpt->GetPropertiesMap() ;
0436     for(MSV::const_iterator it=msv->begin() ; it != msv->end() ; it++ ) ss << it->first << std::endl ;
0437
0438     ss << " GetPropertiesCMap " << std::endl ;
0439     typedef std::map< G4String, G4double, std::less<G4String> > MSC ;
0440     MSC* msc = mpt->GetPropertiesCMap();
0441     for(MSC::const_iterator it=msc->begin() ; it != msc->end() ; it++ ) ss << it->first << std::endl ;
0442     */
0443 #if G4VERSION_NUMBER < 1100
0444     G4MaterialPropertiesTable* mpt = mat->GetMaterialPropertiesTable();
0445     ss << " GetPropertyMap " << std::endl ;
0446     typedef std::map<G4int, G4MaterialPropertyVector*, std::less<G4int> > MIV ;
0447     const MIV* miv =  mpt->GetPropertyMap();
0448     for(MIV::const_iterator it=miv->begin() ; it != miv->end() ; it++ ) ss << pnames[it->first] << std::endl ;
0449
0450     ss << " GetConstPropertyMap " << std::endl ;
0451     typedef std::map<G4int, G4double, std::less<G4int> > MIC ;
0452     const MIC* mic =  mpt->GetConstPropertyMap();
0453     for(MIC::const_iterator it=mic->begin() ; it != mic->end() ; it++ ) ss << cnames[it->first] << std::endl ;
0454
0455 #endif
0456     std::string s = ss.str();
0457     return s ;
0458 }
0459
0460
0461
0462
0463
0464 /**
0465 U4Material::MakeMaterialPropertiesTable
0466 -----------------------------------------
0467
0468 The matdir may use envvar tokens, eg::
0469
0470    $HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/stree/material/LS
0471    $HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/stree/material/Water
0472
0473 **/
0474
0475 G4MaterialPropertiesTable* U4Material::MakeMaterialPropertiesTable(const char* matdir_  )
0476 {
0477     const char* ekey = "U4Material__MakeMaterialPropertiesTable_sigint" ;
0478     bool sigint = ssys::getenvbool(ekey) ;
0479
0480
0481     const char* matdir = spath::Resolve(matdir_);
0482
0483     std::vector<std::string> names ;
0484     SDir::List(names, matdir, ".npy", sigint );
0485
0486     LOG_IF(error, names.size() == 0)
0487        << "Failed to read anything from directory "
0488        << " matdir[" << matdir << "] "
0489        << "set envvar to raise SIGINT : " << ekey
0490        ;
0491
0492     G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable();
0493
0494     std::stringstream ss ;
0495     ss << "matdir " << matdir << " names " << names.size() << std::endl ;
0496
0497     for(unsigned i=0 ; i < names.size() ; i++)
0498     {
0499         const char* name = names[i].c_str();
0500         const char* key = SStr::HeadFirst(name, '.');
0501
0502         NP* a = NP::Load(matdir, name  );
0503
0504         char type = Classify(a);
0505         double* values = a->values<double>() ;
0506         ss << Desc(key, a) << std::endl ;
0507
0508
0509 // guessed version fork
0510 #if G4VERSION_NUMBER < 1100
0511         switch(type)
0512         {
0513             case 'C': mpt->AddConstProperty(key, values[1])    ; break ;
0514             case 'P': mpt->AddProperty(key, MakeProperty(a))   ; break ;
0515             case 'F': mpt->AddProperty(key, MakeProperty(a))   ; break ;
0516         }
0517 #else
0518         G4bool createNewKey = true ; // Geant4 11.2.1 throws exception for new keys without this
0519         switch(type)
0520         {
0521             case 'C': mpt->AddConstProperty(key, values[1], createNewKey)    ; break ;
0522             case 'P': mpt->AddProperty(key, MakeProperty(a), createNewKey)   ; break ;
0523             case 'F': mpt->AddProperty(key, MakeProperty(a), createNewKey)   ; break ;
0524         }
0525 #endif
0526     }
0527
0528     std::string s = ss.str();
0529
0530     LOG(LEVEL) << s ;
0531     //std::cout << s << std::endl ;
0532
0533     return mpt ;
0534 }
0535
0536
0537
0538 /**
0539 U4Material::MakeMaterialPropertiesTable
0540 ----------------------------------------
0541
0542 The matdir may use the standard tokens eg::
0543
0544    $IDPath/GScintillatorLib/LS_ori
0545    $IDPath/GMaterialLib/Water_ori
0546
0547 Constructs the properties table by loading the  properties listed in *delim* delimited *keys_*.
0548
0549 **/
0550
0551 G4MaterialPropertiesTable* U4Material::MakeMaterialPropertiesTable(const char* matdir_, const char* keys_, char delim )
0552 {
0553     const char* matdir = spath::Resolve(matdir_);
0554     //const char* matdir = SPath::Resolve(matdir_, NOOP);
0555
0556     std::vector<std::string> keys ;
0557     SStr::Split( keys_, delim, keys );
0558
0559     G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable();
0560     for(unsigned i=0 ; i < keys.size() ; i++)
0561     {
0562         const std::string& key_ = keys[i];
0563         const char* key = key_.c_str();
0564         const char* name = SStr::Format("%s.npy", key );
0565         NP* a = NP::Load(matdir, name );
0566         assert(a);
0567
0568         LOG(info)
0569              << " key " << std::setw(15) << key
0570              << " a.desc " << a->desc()
0571              ;
0572
0573         G4MaterialPropertyVector* v = MakeProperty(a);
0574         mpt->AddProperty(key, v);
0575     }
0576     return mpt ;
0577 }
0578
0579
0580
0581
0582 G4MaterialPropertyVector* U4Material::MakeProperty(const NP* a)  // static
0583 {
0584     std::vector<double> d, v ;
0585     a->psplit<double>(d,v);   // split array into domain and values
0586     assert(d.size() == v.size());
0587     //assert(d.size() > 1 );  // OpticalCONSTANT (scint time fraction property misusage) breaks this
0588     G4MaterialPropertyVector* mpv = new G4MaterialPropertyVector(d.data(), v.data(), d.size() );
0589     return mpv ;
0590 }
0591
0592 G4MaterialPropertyVector* U4Material::MakeProperty( double value ) // static
0593 {
0594     NP* a = MakePropertyArray(value);
0595     return MakeProperty(a);
0596 }
0597
0598 NP* U4Material::MakePropertyArray( double value ) // static
0599 {
0600     sdomain sdom ;
0601     std::vector<double> dom(sdom.energy_eV, sdom.energy_eV+sdom.length) ;
0602     for(unsigned i=0 ; i < dom.size() ; i++)  dom[i] *= eV ;
0603     std::reverse( dom.begin(), dom.end() );
0604
0605     NP* a = NP::Make<double>(sdom.length, 2);
0606     double* vv = a->values<double>();
0607     for(int i=0 ; i < sdom.length ; i++)
0608     {
0609         vv[i*2+0] = dom[i] ;
0610         vv[i*2+1] = value  ;
0611     }
0612     return a ;
0613 }
0614
0615 /**
0616 U4Material::MakeStandardArray
0617 -----------------------------
0618
0619 Canonically invoked from U4Tree::initMaterials
0620
0621 Hmm this is just accessing the MPT of each material and
0622 getting the standard properties or defaults when not present.
0623
0624 But this approach is discrepant for Water which has its
0625 RAYLEIGH property derived from RINDEX with the physics table
0626 of G4OpRayleigh process.
0627
0628 So need to override the source of some props.
0629
0630 **/
0631
0632 NP* U4Material::MakeStandardArray(
0633     std::vector<const G4Material*>& mats,
0634     const std::map<std::string,G4PhysicsVector*>& prop_override ) // static
0635 {
0636     sdomain dom ;
0637     const sproplist* pl = sproplist::Material() ;
0638
0639     int ni = mats.size() ;
0640     int nj = sprop::NUM_PAYLOAD_GRP ;
0641     int nk = dom.length ;
0642     int nl = sprop::NUM_PAYLOAD_VAL ;
0643
0644     NP* a = NP::Make<double>(ni, nj, nk, nl );
0645     double* aa = a->values<double>() ;
0646
0647     std::vector<std::string> names ;
0648
0649     typedef std::map<std::string,int> SI ;
0650     SI override_count ;
0651
0652     for(int i=0 ; i < ni ; i++)
0653     {
0654         const G4Material* mat = mats[i] ;
0655         const G4String& name = mat->GetName() ;
0656         names.push_back(name.c_str()) ;
0657         G4MaterialPropertiesTable* mpt = mat->GetMaterialPropertiesTable();
0658         LOG_IF(LEVEL, mpt == nullptr ) << "U4Material::MakeStandardArray NO MPT " << name ;
0659
0660         std::array<std::string,8> specs ;
0661         assert( nj*nl == 8 );
0662         for(int j=0 ; j < nj ; j++)  // payload groups
0663         {
0664             for(int l=0 ; l < nl ; l++)   // payload values
0665             {
0666                 const sprop* p = pl->get(j,l) ;
0667                 assert( p );
0668                 const char* key = p->name ;
0669
0670                 int prop_idx = j*nl + l ;
0671
0672                 std::stringstream ss ;
0673                 ss << name << "/" << key ;
0674                 std::string spec_ = ss.str() ; // eg "Water/RAYLEIGH"
0675                 const char* spec = spec_.c_str();
0676                 specs[prop_idx] = spec ;
0677             }
0678         }
0679
0680
0681         for(int j=0 ; j < nj ; j++)           // payload groups
0682         {
0683             for(int k=0 ; k < nk ; k++)       // energy or wavelength domain
0684             {
0685                 double energy_eV = dom.energy_eV[k] ;
0686                 double energy = energy_eV * eV ;
0687
0688                 for(int l=0 ; l < nl ; l++)   // payload values
0689                 {
0690                     const sprop* p = pl->get(j,l) ;
0691                     const char* key = p->name ;
0692
0693                     G4PhysicsVector* mpt_prop = ( mpt ? mpt->GetProperty(key) : nullptr );
0694
0695                     int prop_idx = j*nl + l ;
0696                     const char* spec = specs[prop_idx].c_str() ;
0697                     bool has_override = prop_override.count(spec) > 0 ;
0698                     if(has_override) override_count[spec] += 1 ;
0699
0700                     G4PhysicsVector* prop = has_override ? prop_override.at(spec) : mpt_prop ;
0701                     double value = prop ? prop->Value(energy) : p->def ;
0702
0703                     int index = i*nj*nk*nl + j*nk*nl + k*nl + l ;
0704                     aa[index] = value ;
0705                 }
0706             }
0707         }
0708     }
0709     a->set_names(names) ;
0710
0711     std::stringstream oc ;
0712     for(SI::const_iterator it=override_count.begin() ; it != override_count.end() ; it++)
0713         oc << it->first << ":" << it->second << "," ;
0714     std::string oc_report = oc.str();
0715     a->set_meta<std::string>("override_count", oc_report );
0716
0717     return a ;
0718 }
0719
0720
0721
0722
0723
0724 /**
0725 U4Material::Classify
0726 ---------------------
0727
0728 Heuristic to distinguish ConstProperty from Property based on array size and content
0729 and identify fractional property.
0730
0731 C
0732    ConstProperty
0733 P
0734    Property
0735 F
0736    Property that looks like a fractional split with a small number of values summing to 1.
0737
0738
0739 HMM: could do explicitly with metadata keys ?
0740
0741 **/
0742
0743 char U4Material::Classify(const NP* a)
0744 {
0745     assert(a);
0746     assert(a->shape.size() == 2);
0747     assert(a->has_shape(-1,2) && a->ebyte == 8 && a->uifc == 'f' );
0748     const double* values = a->cvalues<double>() ;
0749     char type = a->has_shape(2,2) && values[1] == values[3] ? 'C' : 'P' ;
0750
0751     int numval = a->shape[0]*a->shape[1] ;
0752     double fractot = 0. ;
0753     if(type == 'P' && numval <= 8)
0754     {
0755         for(int i=0 ; i < numval ; i++) if(i%2==1) fractot += values[i] ;
0756         if( fractot == 1. ) type = 'F' ;
0757     }
0758     return type ;
0759 }
0760
0761
0762 std::string U4Material::Desc(const char* key, const NP* a )
0763 {
0764     char type = Classify(a);
0765     const double* values = a->cvalues<double>() ;
0766     int numval = a->shape[0]*a->shape[1] ;
0767     double value = type == 'C' ? values[1] : 0. ;
0768
0769     std::stringstream ss ;
0770
0771     ss << std::setw(20) << key
0772        << " " << std::setw(10) << a->sstr()
0773        << " type " << type
0774        << " value " << std::setw(10) << std::fixed << std::setprecision(3) << value
0775        << " : "
0776        ;
0777
0778     double tot = 0. ;
0779     if(numval <= 8)
0780     {
0781         for(int i=0 ; i < numval ; i++) if(i%2==1) tot += values[i] ;
0782         for(int i=0 ; i < numval ; i++) ss << std::setw(10) << std::fixed << std::setprecision(3) << values[i] << " " ;
0783         ss << " tot: " << std::setw(10) << std::fixed << std::setprecision(3) << tot << " " ;
0784     }
0785     std::string s = ss.str();
0786     return s ;
0787 }
0788
0789
0790
0791
0792
0793
0794 G4MaterialPropertiesTable*  U4Material::MakeMaterialPropertiesTable_FromProp(
0795      const char* a_key, const G4MaterialPropertyVector* a_prop,
0796      const char* b_key, const G4MaterialPropertyVector* b_prop,
0797      const char* c_key, const G4MaterialPropertyVector* c_prop,
0798      const char* d_key, const G4MaterialPropertyVector* d_prop
0799 )
0800 {
0801     G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable();
0802     if(a_key && a_prop) mpt->AddProperty(a_key,const_cast<G4MaterialPropertyVector*>(a_prop));    //  HUH: why not const ?
0803     if(b_key && b_prop) mpt->AddProperty(b_key,const_cast<G4MaterialPropertyVector*>(b_prop));
0804     if(c_key && c_prop) mpt->AddProperty(c_key,const_cast<G4MaterialPropertyVector*>(c_prop));
0805     if(d_key && d_prop) mpt->AddProperty(d_key,const_cast<G4MaterialPropertyVector*>(d_prop));
0806     return mpt ;
0807 }
0808
0809
0810
0811
0812 G4Material* U4Material::MakeWater(const char* name)  // static
0813 {
0814     G4double a, z, density;
0815     G4int nelements;
0816     G4Element* O = new G4Element("Oxygen"  , "O", z=8 , a=16.00*CLHEP::g/CLHEP::mole);
0817     G4Element* H = new G4Element("Hydrogen", "H", z=1 , a=1.01*CLHEP::g/CLHEP::mole);
0818     G4Material* mat = new G4Material(name, density= 1.0*CLHEP::g/CLHEP::cm3, nelements=2);
0819     mat->AddElement(H, 2);
0820     mat->AddElement(O, 1);
0821     return mat ;
0822 }
0823
0824
0825 G4Material* U4Material::MakeMaterial(const G4MaterialPropertyVector* rindex, const char* name)  // static
0826 {
0827     G4Material* mat = MakeWater(name);
0828     G4MaterialPropertiesTable* mpt = MakeMaterialPropertiesTable_FromProp("RINDEX", rindex) ;
0829     mat->SetMaterialPropertiesTable(mpt) ;
0830     return mat ;
0831 }
0832
0833 /**
0834 U4Material::MakeMaterial
0835 --------------------------
0836
0837 ::
0838
0839     G4Material* mat = MakeMaterial("LS", "$IDPath/GScintillatorLib/LS_ori", "RINDEX,FASTCOMPONENT,SLOWCOMPONENT,REEMISSIONPROB,GammaCONSTANT,OpticalCONSTANT");
0840
0841 **/
0842
0843
0844 G4Material* U4Material::MakeMaterial(const char* name, const char* matdir, const char* props )
0845 {
0846     G4Material* mat = MakeWater(name);
0847     G4MaterialPropertiesTable* mpt = MakeMaterialPropertiesTable(matdir, props, ',');
0848     mat->SetMaterialPropertiesTable(mpt) ;
0849     return mat ;
0850 }
0851
0852 G4Material* U4Material::MakeMaterial(const char* name, const char* matdir)
0853 {
0854     G4Material* mat = MakeWater(name);
0855     G4MaterialPropertiesTable* mpt = MakeMaterialPropertiesTable(matdir);
0856     mat->SetMaterialPropertiesTable(mpt) ;
0857     return mat ;
0858 }
0859
0860
0861
0862
0863 /**
0864 U4Material::MakeScintillator : Only used for function tests/debugging
0865 -------------------------------------------------------------------------
0866
0867 Creates material with Water G4Element and density and then loads the properties of LS into its MPT
0868
0869 **/
0870
0871 G4Material* U4Material::MakeScintillator()
0872 {
0873     const char* matdir = "$HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/stree/material/LS" ;
0874     G4Material* mat = MakeMaterial("LS", matdir );
0875     return mat ;
0876 }
0877
0878
0879 /**
0880 U4Material::LoadOri
0881 ---------------------
0882
0883 Currently original material properties are not standardly persisted,
0884 so typically will have to override the IDPath envvar for this to
0885 succeed to load.
0886
0887 **/
0888
0889 G4Material* U4Material::LoadOri(const char* name)
0890 {
0891     const char* matdir = SStr::Format("%s/%s_ori", LIBDIR, name );
0892     G4Material* mat = MakeMaterial(name, matdir );
0893     return mat ;
0894 }
0895
0896 void U4Material::ListOri(std::vector<std::string>& names)
0897 {
0898     //const char* libdir = SPath::Resolve(LIBDIR, NOOP );
0899     const char* libdir = spath::Resolve(LIBDIR);
0900     const char* ext = "_ori" ;
0901     SDir::List(names, libdir, ext );
0902     SDir::Trim(names, ext);
0903     LOG(LEVEL)
0904         << " libdir " << libdir
0905         << " ext " << ext
0906         << " names.size " << names.size()
0907         ;
0908 }
0909
0910 void U4Material::LoadOri()
0911 {
0912     size_t num_mat[2] ;
0913     num_mat[0] = G4Material::GetNumberOfMaterials();
0914     std::vector<std::string> names ;
0915     ListOri(names);
0916
0917     for(int i=0 ; i < int(names.size()) ; i++)
0918     {
0919         const std::string& name = names[i] ;
0920         const char* n = name.c_str();
0921         LOG(LEVEL) << n ;
0922         G4Material* mat = LoadOri(n);
0923         const G4String& name_ = mat->GetName();
0924         bool name_expect = strcmp(name_.c_str(), n ) == 0 ;
0925         assert(name_expect );
0926         if(!name_expect) std::raise(SIGINT);
0927
0928     }
0929     num_mat[1] = G4Material::GetNumberOfMaterials();
0930     LOG(info) << "Increased G4Material::GetNumberOfMaterials from " << num_mat[0] << " to " << num_mat[1] ;
0931
0932 }
0933
0934
0935
0936 /**
0937 U4Material::LoadBnd from $CFBase/CSGFoundry/SSim/bnd.npy so must have already saved it
0938 ----------------------------------------------------------------------------------------
0939
0940 SSim::Load loads from $CFBase/CSGFoundry/SSim where "$CFBase" is an
0941 internal envvar that yield the result of SOpticksResource::CFBase()
0942 The OPTICKS_KEY derived CSG_GGeo directory is returned unless
0943 the CFBASE envvar is defined.
0944
0945 HMM: if the material exists already then need to change its
0946 properties, not scrub the pre-existing material.
0947 Thats needed for scintillators as the standard bnd properties
0948 do not include all the needed scint props.
0949
0950 Load the material properties from the SSim::get_bnd array using SBnd::getPropertyGroup
0951 for each material.
0952
0953 **/
0954
0955 void U4Material::LoadBnd(const char* ssimdir)
0956 {
0957     SSim* sim = SSim::Load(ssimdir);
0958     const SBnd* sb = sim->get_sbnd();
0959     if(sb == nullptr)
0960     {
0961         LOG(fatal) << "failed to load bnd.npy from ssimdir " << ssimdir ;
0962         return ;
0963     }
0964
0965     std::vector<std::string> mnames ;
0966     sb->getMaterialNames(mnames );
0967
0968     const sproplist* pl = sproplist::Material();
0969
0970     std::vector<std::string> pnames ;
0971     pl->getNames(pnames);
0972
0973     sdomain sdom ;
0974     std::vector<double> dom(sdom.energy_eV, sdom.energy_eV+sdom.length) ;
0975     for(unsigned i=0 ; i < dom.size() ; i++)  dom[i] *= eV ;
0976
0977
0978     for(unsigned m=0 ; m < mnames.size() ; m++)
0979     {
0980         const char* mname = mnames[m].c_str() ;
0981         bool warning = false ;
0982         G4Material* prior = G4Material::GetMaterial(mname, warning);
0983         if(prior != nullptr) std::cout << "prior material " << mname << std::endl ;
0984
0985         G4Material* mat = prior ? prior : MakeWater(mname) ;
0986         G4MaterialPropertiesTable* mpt = mat->GetMaterialPropertiesTable()  ;
0987         if(mpt == nullptr)
0988         {
0989             mpt = new G4MaterialPropertiesTable ;
0990             mat->SetMaterialPropertiesTable(mpt);
0991         }
0992
0993         for(unsigned p=0 ; p < pnames.size() ; p++)
0994         {
0995             const char* pname = pnames[p].c_str() ;
0996             std::vector<double> val ;
0997             sb->getProperty(val, mname, pname);
0998             assert( val.size() == dom.size() );
0999
1000             bool reverse = true ;
1001             AddProperty(mpt, pname, dom, val, reverse );
1002         }
1003     }
1004 }
1005
1006
1007 int U4Material::GetIndex(const std::vector<G4String>& nn, const char* key ) // static
1008 {
1009     G4String k(key);
1010     typedef std::vector<G4String> VS ;
1011     typedef VS::const_iterator   VSI ;
1012     VSI b = nn.begin() ;
1013     VSI e = nn.end() ;
1014     VSI p = std::find(b, e, k );
1015     return p == e ? -1 : std::distance(b, p) ;
1016 }
1017
1018 int U4Material::GetPropertyIndex( const G4MaterialPropertiesTable* mpt, const char* key ) // static
1019 {
1020     const std::vector<G4String> names = mpt->GetMaterialPropertyNames() ;
1021     return GetIndex(names, key);
1022 }
1023
1024 G4MaterialPropertyVector* U4Material::AddProperty( G4MaterialPropertiesTable* mpt, const char* key, const std::vector<double>& dom, const std::vector<double>& val , bool reverse )
1025 {
1026      std::vector<double> ldom(dom);
1027      std::vector<double> lval(val);
1028      assert( ldom.size() == lval.size() );
1029      if(reverse)
1030      {
1031          std::reverse(ldom.begin(), ldom.end());
1032          std::reverse(lval.begin(), lval.end());
1033      }
1034      unsigned numval = ldom.size();
1035      double* ldom_v = ldom.data() ;
1036      double* lval_v = lval.data() ;
1037
1038 #if G4VERSION_NUMBER < 1100
1039     G4MaterialPropertyVector* mpv = mpt->AddProperty(key, ldom_v, lval_v, numval );
1040 #else
1041     int keyIdx = GetPropertyIndex(mpt, key);
1042     G4bool createNewKey = keyIdx == -1  ;
1043     G4MaterialPropertyVector* mpv = mpt->AddProperty(key, ldom_v, lval_v, numval, createNewKey );
1044 #endif
1045     return mpv ;
1046 }
1047
1048
1049 /**
1050 U4Material::KludgeRemoveRockRINDEX
1051 -----------------------------------
1052
1053 Removing the Rock RINDEX property is a kludge that makes photons immediately get absorbed
1054 on reaching the Rock. However this kludge is problematic to random aligned running as
1055 the random consumption does not then follow the normal pattern.
1056
1057 In order for the tail random consumption to be amenable to aligning with Opticks
1058 instead add perfect absorber surfaces to have a cleaner termination instead of the "kill"
1059 that happens with NoRINDEX.
1060
1061 For example U4VolumeMaker::RaindropRockAirWater
1062 adds a surface using U4Surface::MakePerfectAbsorberSurface in
1063
1064 **/
1065
1066 void U4Material::KludgeRemoveRockRINDEX() // static
1067 {
1068     RemoveProperty( "RINDEX", G4Material::GetMaterial("Rock") );
1069 }
1070
1071