From ec6f1805b32c4201c7363997fa910cabee8700a7 Mon Sep 17 00:00:00 2001
From: Markus Frank <Markus.Frank@cern.ch>
Date: Mon, 21 Mar 2022 14:06:16 +0100
Subject: [PATCH] Handle material properties (and const properties) beyond the
G4 material table
---
DDCore/src/plugins/StandardPlugins.cpp | 3 +-
DDG4/src/Geant4Converter.cpp | 99 ++++++++++---------
.../ClientTests/src/MaterialTester_geo.cpp | 40 ++++++++
.../compact/ReadMaterialProperties.xml | 37 +++++++
.../scripts/ReadMaterialProperties.py | 75 ++++++++++++++
5 files changed, 204 insertions(+), 50 deletions(-)
create mode 100644 examples/OpticalSurfaces/scripts/ReadMaterialProperties.py
diff --git a/DDCore/src/plugins/StandardPlugins.cpp b/DDCore/src/plugins/StandardPlugins.cpp
index 34211d8bf..ffd0d15cd 100644
--- a/DDCore/src/plugins/StandardPlugins.cpp
+++ b/DDCore/src/plugins/StandardPlugins.cpp
@@ -139,7 +139,8 @@ static long display(Detector& description, int argc, char** argv) {
" -detector <string> Top level DetElement path. Default: '/world' \n"
" -option <string> ROOT Draw option. Default: 'ogl' \n"
" -level <number> Visualization level [TGeoManager::SetVisLevel] Default: 4 \n"
- " -visopt <number> Visualization option [TGeoManager::SetVisOption] Default: 1 \n"
+ " -visopt <number> Visualization option [TGeoManager::SetVisOption] Default: 1\n"
+ " -load Only load the geometry. Do not invoke the display \n"
" -help Print this help output \n"
" Arguments given: " << arguments(argc,argv) << endl << flush;
::exit(EINVAL);
diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp
index a0e44ad0d..2158a9d2a 100644
--- a/DDG4/src/Geant4Converter.cpp
+++ b/DDG4/src/Geant4Converter.cpp
@@ -339,6 +339,7 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
state = kStateUndefined;
break;
}
+ printout(lvl,"Geant4Material","+++ Setting up material %s", name.c_str());
if ( material->IsMixture() ) {
double A_total = 0.0;
double W_total = 0.0;
@@ -354,8 +355,8 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
TGeoElement* e = mix->GetElement(i);
G4Element* g4e = (G4Element*) handleElement(e->GetName(), Atom(e));
if (!g4e) {
- printout(ERROR, "Material",
- "Missing component %s for material %s. A=%f W=%f",
+ printout(ERROR, name,
+ "Missing element component %s for material %s. A=%f W=%f",
e->GetName(), mix->GetName(), A_total, W_total);
}
//mat->AddElement(g4e, (mix->GetAmixt())[i] / A_total);
@@ -379,15 +380,15 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
TGDMLMatrix* matrix = info.manager->GetGDMLMatrix(named->GetTitle());
const char* cptr = ::strstr(matrix->GetName(), GEANT4_TAG_IGNORE);
if ( 0 != cptr ) {
- printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
- matrix->GetName(), matrix->GetTitle());
- continue;
+ printout(INFO,name,"++ Ignore property %s [%s]. Not Suitable for Geant4.",
+ matrix->GetName(), matrix->GetTitle());
+ continue;
}
cptr = ::strstr(matrix->GetTitle(), GEANT4_TAG_IGNORE);
if ( 0 != cptr ) {
- printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
- matrix->GetName(), matrix->GetTitle());
- continue;
+ printout(INFO,name,"++ Ignore property %s [%s]. Not Suitable for Geant4.",
+ matrix->GetName(), matrix->GetTitle());
+ continue;
}
Geant4GeometryInfo::PropertyVector* v =
(Geant4GeometryInfo::PropertyVector*)handleMaterialProperties(matrix);
@@ -412,7 +413,7 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
}
if ( idx < 0 ) {
printout(ERROR, "Geant4Converter",
- "++ UNKNOWN Geant4 CONST Property: %-20s %s [IGNORED]",
+ "++ UNKNOWN Geant4 Property: %-20s %s [IGNORED]",
exc_str.c_str(), named->GetName());
continue;
}
@@ -423,13 +424,13 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
bins[i] *= conv.first, vals[i] *= conv.second;
G4MaterialPropertyVector* vec =
- new G4MaterialPropertyVector(&bins[0], &vals[0], bins.size());
+ new G4MaterialPropertyVector(&bins[0], &vals[0], bins.size());
tab->AddProperty(named->GetName(), vec);
- printout(lvl, "Geant4Converter", "++ Property: %-20s [%ld x %ld] -> %s ",
+ printout(lvl, name, "++ Property: %-20s [%ld x %ld] -> %s ",
named->GetName(), matrix->GetRows(), matrix->GetCols(), named->GetTitle());
for(size_t i=0, count=v->bins.size(); i<count; ++i)
- printout(lvl,named->GetName()," Geant4: %8.3g [MeV] TGeo: %8.3g [GeV] Conversion: %8.3g",
- bins[i], v->bins[i], conv.first);
+ printout(lvl, name, " Geant4: %s %8.3g [MeV] TGeo: %8.3g [GeV] Conversion: %8.3g",
+ named->GetName(), bins[i], v->bins[i], conv.first);
}
/// Attach the material properties if any
TListIter cpropIt(&material->GetConstProperties());
@@ -440,19 +441,19 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
const char* cptr = ::strstr(named->GetName(), GEANT4_TAG_IGNORE);
if ( 0 != cptr ) {
- printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
- named->GetName(), named->GetTitle());
- continue;
+ printout(INFO, name, "++ Ignore CONST property %s [%s].",
+ named->GetName(), named->GetTitle());
+ continue;
}
cptr = ::strstr(named->GetTitle(), GEANT4_TAG_IGNORE);
if ( 0 != cptr ) {
- printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
- named->GetName(), named->GetTitle());
- continue;
+ printout(INFO, name,"++ Ignore CONST property %s [%s].",
+ named->GetName(), named->GetTitle());
+ continue;
}
double v = info.manager->GetProperty(named->GetTitle(),&err);
if ( err != kFALSE ) {
- except("Geant4Converter",
+ except(name,
"++ FAILED to create G4 material %s "
"[Cannot convert const property: %s]",
material->GetName(), named->GetName());
@@ -473,14 +474,14 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
idx = -1;
}
if ( idx < 0 ) {
- printout(ERROR, "Geant4Converter",
+ printout(ERROR, name,
"++ UNKNOWN Geant4 CONST Property: %-20s %s [IGNORED]",
exc_str.c_str(), named->GetName());
continue;
}
// We need to convert the property from TGeo units to Geant4 units
double conv = g4ConstPropertyConversion(idx);
- printout(lvl, "Geant4Converter", "++ CONST Property: %-20s %g ", named->GetName(), v);
+ printout(lvl, name, "++ CONST Property: %-20s %g ", named->GetName(), v);
tab->AddConstProperty(named->GetName(), v * conv);
}
#endif
@@ -491,7 +492,7 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const
str << " log(MEE): " << ionization->GetLogMeanExcEnergy();
else
str << " log(MEE): UNKNOWN";
- printout(lvl, "Geant4Converter", "++ Created G4 material %s", str.str().c_str());
+ printout(lvl, name, "++ Created G4 material %s", str.str().c_str());
info.g4Materials[medium] = mat;
}
return mat;
@@ -562,9 +563,9 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) c
G4Scale3D scal(vals[0], vals[1], vals[2]);
G4VSolid* g4solid = (G4VSolid*)handleSolid(sol->GetName(), sol);
if ( scal.xx()>0e0 && scal.yy()>0e0 && scal.zz()>0e0 )
- solid = new G4ScaledSolid(sh->GetName(), g4solid, scal);
+ solid = new G4ScaledSolid(sh->GetName(), g4solid, scal);
else
- solid = new G4ReflectedSolid(g4solid->GetName()+"_refl", g4solid, scal);
+ solid = new G4ReflectedSolid(g4solid->GetName()+"_refl", g4solid, scal);
#else
except("Geant4Converter","++ TGeoScaledShape are only supported by Geant4 for versions >= 10.3");
#endif
@@ -753,16 +754,16 @@ void* Geant4Converter::handleAssembly(const string& name, const TGeoNode* node)
MyTransform3D transform(tr->GetTranslation(),tr->IsRotation() ? tr->GetRotationMatrix() : s_identity_rot);
if ( is_left_handed(tr) ) {
- G4Scale3D scale;
- G4Rotate3D rot;
- G4Translate3D trans;
- transform.getDecomposition(scale, rot, trans);
- printout(debugReflections ? ALWAYS : lvl, "Geant4Converter",
- "++ Placing reflected ASSEMBLY. dau:%s to mother %s "
- "Tr:x=%8.1f y=%8.1f z=%8.1f Scale:x=%4.2f y=%4.2f z=%4.2f",
+ G4Scale3D scale;
+ G4Rotate3D rot;
+ G4Translate3D trans;
+ transform.getDecomposition(scale, rot, trans);
+ printout(debugReflections ? ALWAYS : lvl, "Geant4Converter",
+ "++ Placing reflected ASSEMBLY. dau:%s to mother %s "
+ "Tr:x=%8.1f y=%8.1f z=%8.1f Scale:x=%4.2f y=%4.2f z=%4.2f",
dau_vol->GetName(), mot_vol->GetName(),
transform.dx(), transform.dy(), transform.dz(),
- scale.xx(), scale.yy(), scale.zz());
+ scale.xx(), scale.yy(), scale.zz());
}
if ( dau_vol->IsA() == TGeoVolumeAssembly::Class() ) {
@@ -817,12 +818,12 @@ void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node)
TGeoMatrix* tr = node->GetMatrix();
if ( !tr ) {
except("Geant4Converter",
- "++ Attempt to handle placement without transformation:%p %s of type %s vol:%p",
- node, node->GetName(), node->IsA()->GetName(), vol);
+ "++ Attempt to handle placement without transformation:%p %s of type %s vol:%p",
+ node, node->GetName(), node->IsA()->GetName(), vol);
}
else if (0 == vol) {
except("Geant4Converter", "++ Unknown G4 volume:%p %s of type %s ptr:%p",
- node, node->GetName(), node->IsA()->GetName(), vol);
+ node, node->GetName(), node->IsA()->GetName(), vol);
}
else {
int copy = node->GetNumber();
@@ -857,9 +858,9 @@ void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node)
printout(lvl, "Geant4Converter", "++ Assembly: makeImprint: dau:%-12s %s in mother %-12s "
"Tr:x=%8.1f y=%8.1f z=%8.1f Scale:x=%4.2f y=%4.2f z=%4.2f",
node->GetName(), node_is_reflected ? "(REFLECTED)" : "",
- mot_vol ? mot_vol->GetName() : "<unknown>",
+ mot_vol ? mot_vol->GetName() : "<unknown>",
transform.dx(), transform.dy(), transform.dz(),
- scale.xx(), scale.yy(), scale.zz());
+ scale.xx(), scale.yy(), scale.zz());
Geant4AssemblyVolume* ass = (Geant4AssemblyVolume*)info.g4AssemblyVolumes[node];
Geant4AssemblyVolume::Chain chain;
chain.emplace_back(node);
@@ -880,12 +881,12 @@ void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node)
copy, // its copy number
checkOverlaps);
printout(debugReflections||debugPlacements ? ALWAYS : lvl, "Geant4Converter",
- "++ Place %svolume %-12s in mother %-12s "
- "Tr:x=%8.1f y=%8.1f z=%8.1f Scale:x=%4.2f y=%4.2f z=%4.2f",
- node_is_reflected ? "REFLECTED " : "", _v.name(),
- mot_vol ? mot_vol->GetName() : "<unknown>",
- transform.dx(), transform.dy(), transform.dz(),
- scale.xx(), scale.yy(), scale.zz());
+ "++ Place %svolume %-12s in mother %-12s "
+ "Tr:x=%8.1f y=%8.1f z=%8.1f Scale:x=%4.2f y=%4.2f z=%4.2f",
+ node_is_reflected ? "REFLECTED " : "", _v.name(),
+ mot_vol ? mot_vol->GetName() : "<unknown>",
+ transform.dx(), transform.dy(), transform.dz(),
+ scale.xx(), scale.yy(), scale.zz());
// First 2 cases can be combined.
// Leave them separated for debugging G4ReflectionFactory for now...
if ( node_is_reflected && !pvPlaced.second )
@@ -1086,13 +1087,13 @@ void* Geant4Converter::handleMaterialProperties(TObject* mtx) const {
if ( 0 != cptr ) { // Check if the property should not be passed to Geant4
printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
- matrix->GetName(), matrix->GetTitle());
+ matrix->GetName(), matrix->GetTitle());
return nullptr;
}
cptr = ::strstr(matrix->GetTitle(), GEANT4_TAG_IGNORE);
if ( 0 != cptr ) { // Check if the property should not be passed to Geant4
printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
- matrix->GetName(), matrix->GetTitle());
+ matrix->GetName(), matrix->GetTitle());
return nullptr;
}
@@ -1109,7 +1110,7 @@ void* Geant4Converter::handleMaterialProperties(TObject* mtx) const {
g4->values.emplace_back(matrix->Get(i,1));
}
printout(lvl, "Geant4Converter",
- "++ Successfully converted material property:%s : %s [%ld rows]",
+ "++ Successfully converted material property:%s : %s [%ld rows]",
matrix->GetName(), matrix->GetTitle(), rows);
info.g4OpticalProperties[matrix] = g4;
}
@@ -1232,7 +1233,7 @@ void* Geant4Converter::handleOpticalSurface(TObject* surface) const {
TGDMLMatrix* matrix = info.manager->GetGDMLMatrix(named->GetTitle());
const char* cptr = ::strstr(matrix->GetName(), GEANT4_TAG_IGNORE);
if ( 0 != cptr ) // Check if the property should not be passed to Geant4
- continue;
+ continue;
if ( 0 == tab ) {
tab = new G4MaterialPropertiesTable();
@@ -1426,7 +1427,7 @@ namespace {
for (typename C::const_iterator i = c.begin(); i != c.end(); ++i) {
const auto& cc = (*i).second;
for (const auto& j : cc) {
- (o->*pmf)(j);
+ (o->*pmf)(j);
}
}
}
diff --git a/examples/ClientTests/src/MaterialTester_geo.cpp b/examples/ClientTests/src/MaterialTester_geo.cpp
index 653ada20e..b84c70023 100644
--- a/examples/ClientTests/src/MaterialTester_geo.cpp
+++ b/examples/ClientTests/src/MaterialTester_geo.cpp
@@ -37,6 +37,15 @@ static Ref_t create_element(Detector& description, xml_h xml_det, SensitiveDetec
Assembly assembly(det_name+"_assembly");
DetElement det(det_name,x_det.typeStr(), x_det.id());
+ if ( x_det.hasChild(_U(box)) ) {
+ xml_det_t xbox = x_det.child(_U(box));
+ Volume box = Volume(det_name+"_box",
+ Box(xbox.x(), xbox.y(), xbox.z()),
+ description.material(xbox.attr<string>(_U(material))));
+ PlacedVolume place = assembly.placeVolume(box);
+ place.addPhysVolID("box",0);
+ }
+
for(xml_coll_t k(x_det,_Unicode(test)); k; ++k) {
xml_comp_t c = k;
Material mat = description.material(c.nameStr());
@@ -66,6 +75,37 @@ static Ref_t create_element(Detector& description, xml_h xml_det, SensitiveDetec
mix->GetAmixt()[i],wmix ? wmix[i] : -1e0);
}
}
+ if ( material->GetNproperties() > 0 ) {
+ printout(INFO,det_name,"+++ Properties: %d", material->GetNproperties());
+ for(Int_t i=0, n=material->GetNproperties(); i<n; ++i) {
+ TGDMLMatrix* m = material->GetProperty(i);
+ printout(INFO,det_name,"+++ \"%s\" [%s] rows:%d cols:%d",
+ m->GetName(), m->GetTitle(), m->GetRows(), m->GetCols());
+ m->Print();
+ }
+ printout(INFO,det_name,"+++ Properties by NAME:");
+ for(Int_t i=0, n=material->GetNproperties(); i<n; ++i) {
+ const auto* name = material->GetProperties().At(i)->GetName();
+ const auto* m = material->GetProperty(name);
+ printout(INFO,det_name,"+++ \"%s\" [%s,%s] cols: %d rows: %d", name,
+ m->GetName(), m->GetTitle(), m->GetCols(), m->GetRows());
+ }
+ }
+ if ( material->GetNconstProperties() > 0 ) {
+ printout(INFO,det_name,"+++ CONST Properties: %d", material->GetNconstProperties());
+ const TList& all = material->GetConstProperties();
+ for(Int_t i=0, n=material->GetNconstProperties(); i<n; ++i) {
+ const TNamed* m = (const TNamed*)all.At(i);
+ double v = material->GetConstProperty(i);
+ printout(INFO,det_name,"+++ \"%s\" [%s] value: %f", m->GetName(), m->GetTitle(), v);
+ }
+ printout(INFO,det_name,"+++ CONST Properties by NAME:");
+ for(Int_t i=0, n=material->GetNconstProperties(); i<n; ++i) {
+ const auto* name = material->GetConstProperties().At(i)->GetName();
+ double v = material->GetConstProperty(name);
+ printout(INFO,det_name,"+++ \"%s\" value: %f", name, v);
+ }
+ }
}
printout(INFO,det_name,"+++ Basic units:");
diff --git a/examples/OpticalSurfaces/compact/ReadMaterialProperties.xml b/examples/OpticalSurfaces/compact/ReadMaterialProperties.xml
index 481291d10..c8c493516 100644
--- a/examples/OpticalSurfaces/compact/ReadMaterialProperties.xml
+++ b/examples/OpticalSurfaces/compact/ReadMaterialProperties.xml
@@ -20,6 +20,10 @@
<comment>Test reading of TGeo's Material Properties</comment>
</info>
+ <debug>
+ <type name="materials" value="1"/>
+ </debug>
+
<define>
<constant name="world_side" value="1*m"/>
<constant name="world_x" value="world_side/2"/>
@@ -200,11 +204,26 @@
4.136*eV
"/>
</ignore>
+
+ <matrix name= "Water__SpecialMine" option="Geant4-ignore" coldim="1" values="
+ 2.034*eV
+ 2.068*eV
+ 2.103*eV
+ 2.139*eV
+ 2.177*eV
+ 2.216*eV
+ "/>
</properties>
<includes>
<gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/>
</includes>
+
+ <properties>
+ <constant name="Birk__Water|Geant4-ignore" value="12.345678"/>
+ <constant name="Water__Mine|Geant4-ignore" value="87.654321"/>
+ </properties>
+
<materials>
<material name="Air">
<D type="density" unit="g/cm3" value="0.0012"/>
@@ -227,6 +246,24 @@
<property name="ABSLENGTH" ref="ABSLENGTH__0x123aff00"/>
<property name="FASTCOMPONENT" ref="FASTCOMPONENT__0x123aff00"/>
<property name="SLOWCOMPONENT" ref="SLOWCOMPONENT__0x123aff00"/>
+ <!-- Properties ignored by Geant4 -->
+ <property name="Property_of_mine" ref="Water__SpecialMine"/>
+ <constant name="BirksConstant" ref="Birk__Water|Geant4-ignore"/>
+ <!-- Constants ignored by Geant4 -->
+ <constant name="Constant_of_mine" ref="Water__Mine|Geant4-ignore"/>
</material>
</materials>
+
+ <readouts>
+ <readout name="Hits">
+ <id>system:8,box:8</id>
+ </readout>
+ </readouts>
+
+ <detectors>
+ <detector id="1" name="MaterialTester" type="MaterialTester" readout="Hits">
+ <box x="50*cm" y="50*cm" z="50*cm" material="Water"/>
+ <test name="Water"/>
+ </detector>
+ </detectors>
</lccdd>
diff --git a/examples/OpticalSurfaces/scripts/ReadMaterialProperties.py b/examples/OpticalSurfaces/scripts/ReadMaterialProperties.py
new file mode 100644
index 000000000..5b4cd3df3
--- /dev/null
+++ b/examples/OpticalSurfaces/scripts/ReadMaterialProperties.py
@@ -0,0 +1,75 @@
+# ==========================================================================
+# AIDA Detector description implementation
+# --------------------------------------------------------------------------
+# Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+# All rights reserved.
+#
+# For the licensing terms see $DD4hepINSTALL/LICENSE.
+# For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+#
+# ==========================================================================
+#
+from __future__ import absolute_import, unicode_literals
+import os
+import sys
+import DDG4
+from DDG4 import OutputLevel as Output
+from g4units import keV
+#
+#
+"""
+
+ dd4hep simulation example setup using the python configuration
+
+ @author M.Frank
+ @version 1.0
+
+"""
+
+
+def run():
+ kernel = DDG4.Kernel()
+ install_dir = os.environ['DD4hepExamplesINSTALL']
+ kernel.loadGeometry(str("file:" + install_dir + "/examples/OpticalSurfaces/compact/ReadMaterialProperties.xml"))
+
+ DDG4.importConstants(kernel.detectorDescription(), debug=False)
+ geant4 = DDG4.Geant4(kernel, tracker='Geant4TrackerCombineAction')
+ geant4.printDetectors()
+ # Configure UI
+ if len(sys.argv) > 1:
+ geant4.setupCshUI(macro=sys.argv[1])
+ else:
+ geant4.setupCshUI()
+
+ # Configure field
+ geant4.setupTrackingField(prt=True)
+ # Configure Event actions
+ prt = DDG4.EventAction(kernel, 'Geant4ParticlePrint/ParticlePrint')
+ prt.OutputLevel = Output.DEBUG
+ prt.OutputType = 3 # Print both: table and tree
+ kernel.eventAction().adopt(prt)
+
+ generator_output_level = Output.INFO
+
+ # Configure G4 geometry setup
+ seq, act = geant4.addDetectorConstruction("Geant4DetectorGeometryConstruction/ConstructGeo")
+ act.DebugMaterials = True
+ act.DebugVolumes = True
+ act.DebugShapes = True
+
+
+ # Setup particle gun
+ gun = geant4.setupGun("Gun", particle='gamma', energy=5 * keV, multiplicity=1)
+ gun.OutputLevel = generator_output_level
+
+ geant4.setupTracker('MaterialTester')
+
+ # Now build the physics list:
+ phys = geant4.setupPhysics('QGSP_BERT')
+ phys.dump()
+
+ geant4.execute()
+
+
+if __name__ == "__main__":
+ run()
--
GitLab