diff --git a/DDCore/include/DD4hep/Detector.h b/DDCore/include/DD4hep/Detector.h index ee6094646c7e7f91c90258860075ec89aa9ed8e5..42f2a9f6f466d9c6871865fde1434f03ad35e402 100644 --- a/DDCore/include/DD4hep/Detector.h +++ b/DDCore/include/DD4hep/Detector.h @@ -54,6 +54,25 @@ namespace dd4hep { class UriReader; } + /// Helper class to access default temperature and pressure + class STD_Conditions { + public: + enum Conventions { + STP = 1<<0, // Standard temperature and pressure (273.15 Kelvin, 1 ATM) + NTP = 1<<1, // Normal temperature and pressure (293.15 Kelvin, 1 ATM) + USER = 1<<2, // Explicitly set before materials are defined (recommended) + USER_SET = 1<<3, + USER_NOTIFIED = 1<<4 + }; + public: + double pressure; + double temperature; + long convention; + bool is_NTP() const { return (convention&NTP) != 0; } + bool is_STP() const { return (convention&STP) != 0; } + bool is_user_defined() const { return (convention&USER) != 0; } + }; + /// The main interface to the dd4hep detector description package /** * Note: The usage of the factory method: @@ -128,6 +147,14 @@ namespace dd4hep { /// Access the optical surface manager virtual OpticalSurfaceManager surfaceManager() const = 0; + /// Access default conditions (temperature and pressure + virtual const STD_Conditions& stdConditions() const = 0; + /// Set the STD temperature and pressure + virtual void setStdConditions(double temp, double pressure) = 0; + /// Set the STD conditions according to defined types (STP or NTP) + virtual void setStdConditions(const std::string& type) = 0; + + /// Accessor to the map of header entries virtual Header header() const = 0; /// Accessor to the header entry diff --git a/DDCore/include/DD4hep/DetectorImp.h b/DDCore/include/DD4hep/DetectorImp.h index c0c462e94023036cddcc9e241383313f07554d9f..f31ff1a22326eb703cb2927049eb2cf598bff081 100644 --- a/DDCore/include/DD4hep/DetectorImp.h +++ b/DDCore/include/DD4hep/DetectorImp.h @@ -61,6 +61,9 @@ namespace dd4hep { /// Cached map with detector types: typedef std::map<std::string, std::vector<DetElement> > DetectorTypeMap; + /// Standard conditions + mutable STD_Conditions m_std_conditions; + /// Inventory of detector types DetectorTypeMap m_detectorTypes; @@ -207,6 +210,14 @@ namespace dd4hep { virtual OverlayedField field() const override { return m_field; } + + /// Access default conditions (temperature and pressure + virtual const STD_Conditions& stdConditions() const override; + /// Set the STD temperature and pressure + virtual void setStdConditions(double temp, double pressure) override; + /// Set the STD conditions according to defined types (STP or NTP) + virtual void setStdConditions(const std::string& type) override; + /// Accessor to the header entry virtual Header header() const override { return m_header; diff --git a/DDCore/include/XML/UnicodeValues.h b/DDCore/include/XML/UnicodeValues.h index b5dd9437b84b3d699062b39092cfaaef859217e8..7ad2367fa447a08fe70541cc8dd5dd042da79272 100644 --- a/DDCore/include/XML/UnicodeValues.h +++ b/DDCore/include/XML/UnicodeValues.h @@ -454,6 +454,7 @@ UNICODE (starttheta); UNICODE (state); UNICODE (stave); UNICODE (staves); +UNICODE (std_conditions); UNICODE (step); UNICODE (steps); UNICODE (step_x); diff --git a/DDCore/src/DetectorImp.cpp b/DDCore/src/DetectorImp.cpp index f4156e445ba2b506a94006947cd68c0b1c1663a8..c19dbe4520803836a2c34301888af59899963c98 100644 --- a/DDCore/src/DetectorImp.cpp +++ b/DDCore/src/DetectorImp.cpp @@ -26,6 +26,7 @@ #include "DD4hep/detail/VolumeManagerInterna.h" #include "DD4hep/detail/OpticalSurfaceManagerInterna.h" #include "DD4hep/DetectorImp.h" +#include "DD4hep/DD4hepUnits.h" // C/C++ include files #include <iostream> @@ -188,7 +189,10 @@ DetectorImp::DetectorImp(const string& name) gGeoIdentity = new TGeoIdentity(); } m_surfaceManager = new detail::OpticalSurfaceManagerObject(*this); - + m_std_conditions.convention = STD_Conditions::NTP; + m_std_conditions.pressure = NTP_Pressure; + m_std_conditions.temperature = NTP_Temperature; + VisAttr attr("invisible"); attr.setColor(0.5, 0.5, 0.5); attr.setAlpha(1); @@ -316,6 +320,48 @@ Volume DetectorImp::pickMotherVolume(const DetElement& de) const { return 0; } +/// Access default conditions (temperature and pressure +const STD_Conditions& DetectorImp::stdConditions() const { + if ( (m_std_conditions.convention&STD_Conditions::USER_SET) == 0 && + (m_std_conditions.convention&STD_Conditions::USER_NOTIFIED) == 0 ) + { + printout(WARNING,"DD4hep","++ STD conditions NOT defined by client. NTP defaults taken."); + m_std_conditions.convention |= STD_Conditions::USER_NOTIFIED; + } + return m_std_conditions; +} +/// Set the STD temperature and pressure +void DetectorImp::setStdConditions(double temp, double pressure) { + m_std_conditions.temperature = temp; + m_std_conditions.pressure = pressure; + m_std_conditions.convention = STD_Conditions::USER_SET; + if ( std::abs(temp-NTP_Temperature) < 1e-10 && std::abs(pressure-NTP_Pressure) < 1e-10 ) + m_std_conditions.convention |= STD_Conditions::NTP; + else if ( std::abs(temp-STP_Temperature) < 1e-10 && std::abs(pressure-STP_Pressure) < 1e-10 ) + m_std_conditions.convention |= STD_Conditions::STP; + else + m_std_conditions.convention |= STD_Conditions::USER; +} + +/// Set the STD conditions according to defined types (STP or NTP) +void DetectorImp::setStdConditions(const std::string& type) { + if ( type == "STP" ) { + m_std_conditions.temperature = STP_Temperature; + m_std_conditions.pressure = STP_Pressure; + m_std_conditions.convention = STD_Conditions::STP|STD_Conditions::USER_SET; + } + else if ( type == "NTP" ) { + m_std_conditions.temperature = NTP_Temperature; + m_std_conditions.pressure = NTP_Pressure; + m_std_conditions.convention = STD_Conditions::NTP|STD_Conditions::USER_SET; + } + else { + except("DD4hep", + "++ Attempt to set standard conditions to " + "unknown conventions (Only STP and NTP allowed)."); + } +} + /// Retrieve a subdetector element by it's name from the detector description DetElement DetectorImp::detector(const std::string& name) const { HandleMap::const_iterator i = m_detectors.find(name); diff --git a/DDCore/src/plugins/Compact2Objects.cpp b/DDCore/src/plugins/Compact2Objects.cpp index 9522469fdb05923c5352aa72302cc37c72c3fbf5..51060e6c2c25d7402cd04da84bc27a27390adab9 100644 --- a/DDCore/src/plugins/Compact2Objects.cpp +++ b/DDCore/src/plugins/Compact2Objects.cpp @@ -72,7 +72,8 @@ namespace dd4hep { class PropertyConstant; class Parallelworld_Volume; class DetElementInclude; - + class STD_Conditions; + /// Converter instances implemented in this compilation unit template <> void Converter<Debug>::operator()(xml_h element) const; template <> void Converter<World>::operator()(xml_h element) const; @@ -95,6 +96,7 @@ namespace dd4hep { template <> void Converter<PropertyConstant>::operator()(xml_h element) const; #endif template <> void Converter<DetElement>::operator()(xml_h element) const; + template <> void Converter<STD_Conditions>::operator()(xml_h element) const; template <> void Converter<GdmlFile>::operator()(xml_h element) const; template <> void Converter<JsonFile>::operator()(xml_h element) const; template <> void Converter<XMLFile>::operator()(xml_h element) const; @@ -416,14 +418,10 @@ template <> void Converter<Material>::operator()(xml_h e) const { dens_unit = density.attr<double>(_U(unit))/xml::_toDouble(_Unicode(gram/cm3)); } if ( dens_unit != 1.0 ) { - cout << matname << " Density unit:" << dens_unit; - if ( dens_unit != 1.0 ) cout << " " << density.attr<string>(_U(unit)); - cout << " Density Value raw:" << dens_val << " normalized:" << (dens_val*dens_unit) << endl; dens_val *= dens_unit; + printout(s_debug.materials ? ALWAYS : DEBUG, "Compact","Density unit: %.3f [%s] raw: %.3f normalized: %.3f ", + dens_unit, density.attr<string>(_U(unit)).c_str(), dens_val, (dens_val*dens_unit)); } - printout(s_debug.materials ? ALWAYS : DEBUG, "Compact", - "++ Converting material %-16s Density: %.3f.",matname, dens_val); - //throw 1; mat = mix = new TGeoMixture(matname, composites.size(), dens_val); size_t ifrac = 0; vector<double> composite_fractions; @@ -458,6 +456,33 @@ template <> void Converter<Material>::operator()(xml_h e) const { else throw_print("Compact2Objects[ERROR]: Converting material:" + mname + " Element missing: " + nam); } + xml_h temperature = x_mat.child(_U(T), false); + double temp_val = description.stdConditions().temperature; + if ( temperature.ptr() ) { + double temp_unit = _toDouble("kelvin"); + if ( temperature.hasAttr(_U(unit)) ) + temp_unit = temperature.attr<double>(_U(unit)); + temp_val = temperature.attr<double>(_U(value)) * temp_unit; + } + xml_h pressure = x_mat.child(_U(P), false); + double pressure_val = description.stdConditions().pressure; + if ( pressure.ptr() ) { + double pressure_unit = _toDouble("pascal"); + if ( pressure.hasAttr(_U(unit)) ) + pressure_unit = pressure.attr<double>(_U(unit)); + pressure_val = pressure.attr<double>(_U(value)) * pressure_unit; + } +#if 0 + printout(s_debug.materials ? ALWAYS : DEBUG, "Compact", + "++ ROOT raw temperature and pressure: %.3g %.3g", + mat->GetTemperature(),mat->GetPressure()); +#endif + mat->SetTemperature(temp_val); + mat->SetPressure(pressure_val); + printout(s_debug.materials ? ALWAYS : DEBUG, "Compact", + "++ Converting material %-16s Density: %9.7g Temperature:%9.7g Pressure:%9.7g.", + matname, dens_val, temp_val, pressure_val); + mix->SetRadLen(0e0); #if ROOT_VERSION_CODE >= ROOT_VERSION(6,12,0) mix->ComputeDerivedQuantities(); @@ -512,18 +537,6 @@ template <> void Converter<Material>::operator()(xml_h e) const { } } #endif - xml_h temp = x_mat.child(_U(T), false); - if ( temp.ptr() ) { - double temp_val = temp.attr<double>(_U(value)); - double temp_unit = temp.attr<double>(_U(unit), 1.0 /* _toDouble("kelvin") */); - mat->SetTemperature(temp_val*temp_unit); - } - xml_h pressure = x_mat.child(_U(P), false); - if ( pressure.ptr() ) { - double pressure_val = pressure.attr<double>(_U(value)); - double pressure_unit = pressure.attr<double>(_U(unit),1.0 /* _toDouble("pascal") */); - mat->SetPressure(pressure_val*pressure_unit); - } } TGeoMedium* medium = mgr.GetMedium(matname); if (0 == medium) { @@ -626,7 +639,7 @@ template <> void Converter<Atom>::operator()(xml_h e) const { string ref = i.attr<string>(_U(ref)); TGeoIsotope* iso = tab->FindIsotope(ref.c_str()); if ( !iso ) { - except("Compact","Element %s cannot be constructed. Isotope '%s' (fraction:%f) missing!", + except("Compact","Element %s cannot be constructed. Isotope '%s' (fraction: %.3f) missing!", name.c_str(), ref.c_str(), frac); } printout(s_debug.elements ? ALWAYS : DEBUG, "Compact", @@ -650,6 +663,43 @@ template <> void Converter<Atom>::operator()(xml_h e) const { } } +/** Convert compact isotope objects + * + * <std_conditions type="STP or NTP"> // type optional + * <item name="temperature" unit="kelvin" value="273.15"/> + * <item name="pressure" unit="kPa" value="100"/> + * </std_conditions> + */ +template <> void Converter<STD_Conditions>::operator()(xml_h e) const { + xml_dim_t cond(e); + // Create the isotope object in the event it is not yet present from the XML data + if ( cond.ptr() ) { + if ( cond.hasAttr(_U(type)) ) { + description.setStdConditions(cond.typeStr()); + } + xml_h temperature = cond.child(_U(T), false); + double temp_val = description.stdConditions().temperature; + if ( temperature.ptr() ) { + double temp_unit = _toDouble("kelvin"); + if ( temperature.hasAttr(_U(unit)) ) + temp_unit = temperature.attr<double>(_U(unit)); + temp_val = temperature.attr<double>(_U(value)) * temp_unit; + } + xml_h pressure = cond.child(_U(P), false); + double pressure_val = description.stdConditions().pressure; + if ( pressure.ptr() ) { + double pressure_unit = _toDouble("pascal"); + if ( pressure.hasAttr(_U(unit)) ) + pressure_unit = pressure.attr<double>(_U(unit)); + pressure_val = pressure.attr<double>(_U(value)) * pressure_unit; + } + description.setStdConditions(temp_val, pressure_val); + printout(s_debug.materials ? ALWAYS : DEBUG, "Compact", + "+++ Material standard conditions: Temperature: %.3f Kelvin Pressure: %.3f hPa", + temp_val/_toDouble("kelvin"), pressure_val/_toDouble("hPa")); + } +} + #if ROOT_VERSION_CODE >= ROOT_VERSION(6,17,0) /** Convert compact optical surface objects (defines) * @@ -668,7 +718,7 @@ template <> void Converter<OpticalSurface>::operator()(xml_h element) const { if ( e.hasAttr(_U(value)) ) value = e.attr<double>(_U(value)); OpticalSurface surf(description, e.attr<string>(_U(name)), model, finish, type, value); if ( s_debug.surface ) { - printout(ALWAYS,"Compact","+++ Reading optical surface %s Typ:%d model:%d finish:%d value:%f", + printout(ALWAYS,"Compact","+++ Reading optical surface %s Typ:%d model:%d finish:%d value: %.3f", e.attr<string>(_U(name)).c_str(), int(type), int(model), int(finish), value); } for (xml_coll_t props(e, _U(property)); props; ++props) { @@ -1142,7 +1192,7 @@ template <> void Converter<SensitiveDetector>::operator()(xml_h element) const { else if (ecut) { // If no unit is given , we assume the correct Geant4 unit is used! sd.setEnergyCutoff(element.attr<double>(ecut)); } - printout(DEBUG, "Compact", "SensitiveDetector-update: %-18s %-24s Hits:%-24s Cutoff:%f7.3f", sd.name(), + printout(DEBUG, "Compact", "SensitiveDetector-update: %-18s %-24s Hits:%-24s Cutoff:%7.3f", sd.name(), (" [" + sd.type() + "]").c_str(), sd.hitsCollection().c_str(), sd.energyCutoff()); xml_attr_t sequence = element.attr_nothrow(_U(sequence)); if (sequence) { @@ -1457,6 +1507,7 @@ template <> void Converter<Compact>::operator()(xml_h element) const { xml_coll_t(compact, _U(define)).for_each(_U(include), Converter<DetElementInclude>(description)); xml_coll_t(compact, _U(define)).for_each(_U(constant), Converter<Constant>(description)); + xml_coll_t(compact, _U(std_conditions)).for_each(Converter<STD_Conditions>(description)); xml_coll_t(compact, _U(includes)).for_each(_U(gdmlFile), Converter<GdmlFile>(description)); if (element.hasChild(_U(info))) diff --git a/DDCore/src/plugins/StandardPlugins.cpp b/DDCore/src/plugins/StandardPlugins.cpp index 58e7d5613dab130b0ab0dba65e549dda1605b8c7..62e135789574d540be2f25b03e19b8c377b29c5b 100644 --- a/DDCore/src/plugins/StandardPlugins.cpp +++ b/DDCore/src/plugins/StandardPlugins.cpp @@ -500,10 +500,20 @@ static long root_materials(Detector& description, int argc, char** argv) { MaterialPrint(Detector& desc) : description(desc) {} virtual ~MaterialPrint() = default; virtual elt_h print(TGeoMaterial* mat) { - ::printf("%-8s %-32s Aeff=%7.3f Zeff=%7.4f rho=%8.3f [g/mole] radlen=%8.3g [cm] intlen=%8.3g [cm] index=%3d\n", - mat->IsMixture() ? "Mixture" : "Material", - mat->GetName(), mat->GetA(), mat->GetZ(), mat->GetDensity(), + const char* st = "Undefined"; + switch(mat->GetState()) { + case TGeoMaterial::kMatStateSolid: st = "solid"; break; + case TGeoMaterial::kMatStateLiquid: st = "liquid"; break; + case TGeoMaterial::kMatStateGas: st = "gas"; break; + case TGeoMaterial::kMatStateUndefined: + default: st = "Undefined"; break; + } + ::printf("%-32s %-8s\n", mat->GetName(), mat->IsMixture() ? "Mixture" : "Material"); + ::printf(" Aeff=%7.3f Zeff=%7.4f rho=%8.3f [g/mole] radlen=%8.3g [cm] intlen=%8.3g [cm] index=%3d\n", + mat->GetA(), mat->GetZ(), mat->GetDensity(), mat->GetRadLen(), mat->GetIntLen(), mat->GetIndex()); + ::printf(" Temp=%.3g [Kelvin] Pressure=%.5g [pascal] state=%s\n", + mat->GetTemperature(), mat->GetPressure(), st); return elt_h(0); } virtual void print(elt_h, TGeoElement* elt, double frac) { @@ -547,11 +557,34 @@ static long root_materials(Detector& description, int argc, char** argv) { } else if ( mat->GetNelements() == 1 ) { elt.setAttr(_U(formula),mat->GetElement(0)->GetName()); +#if 0 + // We do not need this. It shall be computed by TGeo using the Geant4 formula. + elt_h RL = elt.addChild(_U(RL)); + RL.setAttr(_U(type), "X0"); + RL.setAttr(_U(value), mat->GetRadLen()); + RL.setAttr(_U(unit), "cm"); + elt_h NIL = elt.addChild(_U(NIL)); + NIL.setAttr(_U(type), "lambda"); + NIL.setAttr(_U(value), mat->GetIntLen()); + NIL.setAttr(_U(unit), "cm"); +#endif } elt_h D = elt.addChild(_U(D)); - D.setAttr(_U(type),"density"); - D.setAttr(_U(value),mat->GetDensity()); - D.setAttr(_U(unit),"g/cm3"); + D.setAttr(_U(type), "density"); + D.setAttr(_U(value), mat->GetDensity()); + D.setAttr(_U(unit), "g/cm3"); + if ( mat->GetTemperature() != description.stdConditions().temperature ) { + elt_h T = elt.addChild(_U(T)); + T.setAttr(_U(type), "temperature"); + T.setAttr(_U(value), mat->GetTemperature()); + T.setAttr(_U(unit), "kelvin"); + } + if ( mat->GetPressure() != description.stdConditions().pressure ) { + elt_h P = elt.addChild(_U(P)); + P.setAttr(_U(type), "pressure"); + P.setAttr(_U(value), mat->GetPressure()); + P.setAttr(_U(unit), "pascal"); + } return elt; } virtual void print(elt_h mat, TGeoElement* element, double frac) { diff --git a/DDG4/include/DDG4/Geant4UIManager.h b/DDG4/include/DDG4/Geant4UIManager.h index 82d4f6820436c25ec38f12b75ffa16c3536d4dfe..58502337db75d10f08628ebc9f4d78137f2c296a 100644 --- a/DDG4/include/DDG4/Geant4UIManager.h +++ b/DDG4/include/DDG4/Geant4UIManager.h @@ -79,6 +79,8 @@ namespace dd4hep { Geant4UIManager(Geant4Context* context, const std::string& name); /// Default destructor virtual ~Geant4UIManager(); + /// Install command control messenger to write GDML file from command prompt. + void installCommandMessenger(); /// Start visualization G4VisManager* startVis(); /// Start UI @@ -87,7 +89,9 @@ namespace dd4hep { void start(); /// Stop and release resources void stop(); - /// Run UI + /// Force exiting this process without calling atexit handlers + void forceExit(); + /// Run UI virtual void operator()(void* param); }; diff --git a/DDG4/include/DDG4/Geant4UIMessenger.h b/DDG4/include/DDG4/Geant4UIMessenger.h index 045bc1714d8cbec6bc7ede9330621f1fa6f0790f..391adda88286cf94a6e5b76b03715821e1f82ef6 100644 --- a/DDG4/include/DDG4/Geant4UIMessenger.h +++ b/DDG4/include/DDG4/Geant4UIMessenger.h @@ -38,17 +38,17 @@ namespace dd4hep { typedef std::map<G4UIcommand*, Callback> Actions; protected: /// The UI directory of this component - G4UIdirectory* m_directory; + G4UIdirectory* m_directory; /// Reference to the property manager corresponding to the component PropertyManager* m_properties; /// Component name - std::string m_name; + std::string m_name; /// Path in the UI hierarchy of this component - std::string m_path; + std::string m_path; /// Property update command map - Commands m_propertyCmd; + Commands m_propertyCmd; /// Action map - Actions m_actionCmd; + Actions m_actionCmd; public: /// Initializing constructor @@ -56,12 +56,18 @@ namespace dd4hep { /// Default destructor virtual ~Geant4UIMessenger(); /// Add a new callback structure - void addCall(const std::string& name, const std::string& description, const Callback& cb); - /// Add any callback (without parameters to the messenger + void addCall(const std::string& name, const std::string& description, const Callback& cb, size_t npar=0); + /// Add any callback without parameters to the messenger template <typename Q, typename R, typename T> void addCall(const std::string& name, const std::string& description, Q* p, R (T::*f)()) { CallbackSequence::checkTypes(typeid(Q), typeid(T), dynamic_cast<T*>(p)); - addCall(name, description, Callback(p).make(f)); + addCall(name, description, Callback(p).make(f), 0); + } + /// Add any callback with ONE parameter to the messenger + template <typename Q, typename R, typename T, typename A1> + void addCall(const std::string& name, const std::string& description, Q* p, R (T::*f)(A1)) { + CallbackSequence::checkTypes(typeid(Q), typeid(T), dynamic_cast<T*>(p)); + addCall(name, description, Callback(p).make(f), 1); } /// Export all properties to the Geant4 UI void exportProperties(PropertyManager& mgr); diff --git a/DDG4/plugins/Geant4DetectorGeometryConstruction.cpp b/DDG4/plugins/Geant4DetectorGeometryConstruction.cpp index 0eb974ad81ac88d96cd1717553f24d9168aa546f..d794601633b68b407ccab8ff9f690daba67b1c2e 100644 --- a/DDG4/plugins/Geant4DetectorGeometryConstruction.cpp +++ b/DDG4/plugins/Geant4DetectorGeometryConstruction.cpp @@ -58,15 +58,15 @@ namespace dd4hep { int m_geoInfoPrintLevel; /// Property: G4 GDML dump file name (default: empty. If non empty, dump) std::string m_dumpGDML; - /// Property: DD4hep path to volume to be printed (default: empty) - std::string m_volumePath; /// Write GDML file - int writeGDML(); + int writeGDML(const char* gdml_output); /// Print geant4 volume - int printVolume(); + int printVolume(const char* vol_path); /// Check geant4 volume - int checkVolume(); + int checkVolume(const char* vol_path); + /// Print geant4 material + int printMaterial(const char* mat_name); public: /// Initializing constructor for DDG4 @@ -101,11 +101,14 @@ namespace dd4hep { #include <G4Material.hh> #include <G4Version.hh> #include <G4VSolid.hh> +#include "CLHEP/Units/SystemOfUnits.h" //#ifdef GEANT4_HAS_GDML #include <G4GDMLParser.hh> //#endif +#include <cmath> + using namespace std; using namespace dd4hep; using namespace dd4hep::sim; @@ -129,7 +132,6 @@ Geant4DetectorGeometryConstruction::Geant4DetectorGeometryConstruction(Geant4Con declareProperty("DumpHierarchy", m_dumpHierarchy); declareProperty("DumpGDML", m_dumpGDML=""); - declareProperty("VolumePath", m_volumePath=""); InstanceCount::increment(this); } @@ -164,51 +166,83 @@ void Geant4DetectorGeometryConstruction::constructGeo(Geant4DetectorConstruction dmp.dump("",w); } ctxt->world = w; - if ( !m_dumpGDML.empty() || ::getenv("DUMP_GDML") ) writeGDML(); + if ( !m_dumpGDML.empty() ) writeGDML(m_dumpGDML.c_str()); + else if ( ::getenv("DUMP_GDML") ) writeGDML(::getenv("DUMP_GDML")); enableUI(); } +/// Print geant4 material +int Geant4DetectorGeometryConstruction::printMaterial(const char* mat_name) { + if ( mat_name ) { + auto& g4map = Geant4Mapping::instance().data(); + for ( auto it = g4map.g4Materials.begin(); it != g4map.g4Materials.end(); ++it ) { + const auto* mat = (*it).second; + if ( mat->GetName() == mat_name ) { + const auto* ion = mat->GetIonisation(); + printP2("+++ Dump of GEANT4 material: %s", mat_name); + cout << mat; + if ( ion ) { + cout << " MEE: "; + cout << setprecision(12); + cout << ion->GetMeanExcitationEnergy()/CLHEP::eV; + cout << " [eV]"; + } + else + cout << " MEE: UNKNOWN"; + cout << endl << endl; + return 1; + } + } + warning("+++ printMaterial: FAILED to find the material %s", mat_name); + } + warning("+++ printMaterial: Property materialName not set!"); + return 0; +} + /// Print geant4 volume -int Geant4DetectorGeometryConstruction::printVolume() { - if ( !m_volumePath.empty() ) { +int Geant4DetectorGeometryConstruction::printVolume(const char* vol_path) { + if ( vol_path ) { Detector& det = context()->kernel().detectorDescription(); PlacedVolume top = det.world().placement(); - PlacedVolume pv = detail::tools::findNode(top, m_volumePath); + PlacedVolume pv = detail::tools::findNode(top, vol_path); if ( pv.isValid() ) { auto& g4map = Geant4Mapping::instance().data(); auto it = g4map.g4Volumes.find(pv.volume()); if ( it != g4map.g4Volumes.end() ) { const G4LogicalVolume* vol = (*it).second; - auto* sol = vol->GetSolid(); - string txt; - stringstream str; - str << *(vol->GetMaterial()) << endl; - str << *sol; - stringstream istr(str.str()); - printP2("+++ Dump of ROOT solid: %s", m_volumePath.c_str()); - pv.volume().solid()->InspectShape(); - printP2("+++ Dump of GEANT4 solid: %s", m_volumePath.c_str()); - while(istr.good()) { - getline(istr,txt); - printP2(txt.c_str()); + auto* sol = vol->GetSolid(); + const auto* mat = vol->GetMaterial(); + const auto* ion = mat->GetIonisation(); + printP2("+++ Dump of GEANT4 solid: %s", vol_path); + cout << mat; + if ( ion ) { + cout << " MEE: "; + cout << setprecision(12); + cout << ion->GetMeanExcitationEnergy()/CLHEP::eV; + cout << " [eV]"; } - printP2("Shape: %s cubic volume: %8.3g mm^3 area: %8.3g mm^2", + else + cout << " MEE: UNKNOWN"; + cout << endl << *sol; + printP2("+++ Dump of ROOT solid: %s", vol_path); + pv.volume().solid()->InspectShape(); + printP2("+++ Shape: %s cubic volume: %8.3g mm^3 area: %8.3g mm^2", sol->GetName().c_str(), sol->GetCubicVolume(), sol->GetSurfaceArea()); return 1; } } - warning("+++ printVolume: FAILED to find the volume %s from the top volume",m_volumePath.c_str()); + warning("+++ printVolume: FAILED to find the volume %s from the top volume",vol_path); } warning("+++ printVolume: Property VolumePath not set. [Ignored]"); return 0; } /// Check geant4 volume -int Geant4DetectorGeometryConstruction::checkVolume() { - if ( !m_volumePath.empty() ) { +int Geant4DetectorGeometryConstruction::checkVolume(const char* vol_path) { + if ( vol_path ) { Detector& det = context()->kernel().detectorDescription(); PlacedVolume top = det.world().placement(); - PlacedVolume pv = detail::tools::findNode(top, m_volumePath); + PlacedVolume pv = detail::tools::findNode(top, vol_path); if ( pv.isValid() ) { auto& g4map = Geant4Mapping::instance().data(); auto it = g4map.g4Volumes.find(pv.volume()); @@ -231,15 +265,18 @@ int Geant4DetectorGeometryConstruction::checkVolume() { return 1; } } - warning("+++ checkVolume: FAILED to find the volume %s from the top volume",m_volumePath.c_str()); + warning("+++ checkVolume: FAILED to find the volume %s from the top volume",vol_path); } warning("+++ checkVolume: Property VolumePath not set. [Ignored]"); return 0; } /// Write GDML file -int Geant4DetectorGeometryConstruction::writeGDML() { +int Geant4DetectorGeometryConstruction::writeGDML(const char* output) { G4VPhysicalVolume* w = context()->world(); + if ( output && ::strlen(output) > 0 && output != m_dumpGDML.c_str() ) + m_dumpGDML = output; + //#ifdef GEANT4_HAS_GDML if ( !m_dumpGDML.empty() ) { G4GDMLParser parser; @@ -264,12 +301,15 @@ int Geant4DetectorGeometryConstruction::writeGDML() { /// Install command control messenger to write GDML file from command prompt. void Geant4DetectorGeometryConstruction::installCommandMessenger() { this->Geant4DetectorConstruction::installCommandMessenger(); - m_control->addCall("writeGDML", "GDML: write geometry to file: '"+m_dumpGDML+"' [uses property DumpGDML]", - Callback(this).make(&Geant4DetectorGeometryConstruction::writeGDML)); - m_control->addCall("printVolume", "Print Geant4 volume properties [uses property VolumePath]", - Callback(this).make(&Geant4DetectorGeometryConstruction::printVolume)); - m_control->addCall("checkVolume", "Check Geant4 volume properties [uses property VolumePath]", - Callback(this).make(&Geant4DetectorGeometryConstruction::checkVolume)); + m_control->addCall("writeGDML", "GDML: write geometry to file: '"+m_dumpGDML+ + "' [uses argument - or - property DumpGDML]", + Callback(this).make(&Geant4DetectorGeometryConstruction::writeGDML),1); + m_control->addCall("printVolume", "Print Geant4 volume properties [uses argument]", + Callback(this).make(&Geant4DetectorGeometryConstruction::printVolume),1); + m_control->addCall("checkVolume", "Check Geant4 volume properties [uses argument]", + Callback(this).make(&Geant4DetectorGeometryConstruction::checkVolume),1); + m_control->addCall("printMaterial", "Print Geant4 material properties [uses argument]", + Callback(this).make(&Geant4DetectorGeometryConstruction::printMaterial),1); } diff --git a/DDG4/src/Geant4UIManager.cpp b/DDG4/src/Geant4UIManager.cpp index 37ebc5df82509b624f52bfa40262ac82fa1d2b01..9f6098d5d4580db279a8f6fd979b3ef4a0ad9e71 100644 --- a/DDG4/src/Geant4UIManager.cpp +++ b/DDG4/src/Geant4UIManager.cpp @@ -14,6 +14,7 @@ // Framework include files #include "DDG4/Geant4UIManager.h" #include "DDG4/Geant4Kernel.h" +#include "DDG4/Geant4UIMessenger.h" #include "DD4hep/Primitives.h" #include "DD4hep/Printout.h" @@ -25,6 +26,8 @@ #include "G4UIExecutive.hh" #include "G4RunManager.hh" +#include <cstdlib> + using namespace dd4hep::sim; using namespace std; @@ -47,12 +50,24 @@ Geant4UIManager::Geant4UIManager(Geant4Context* ctxt, const std::string& nam) declareProperty("HaveVIS", m_haveVis=false); declareProperty("HaveUI", m_haveUI=true); declareProperty("Prompt", m_prompt); + enableUI(); } /// Default destructor Geant4UIManager::~Geant4UIManager() { } +/// Install command control messenger to write GDML file from command prompt. +void Geant4UIManager::installCommandMessenger() { + m_control->addCall("exit", "Force exiting this process", + Callback(this).make(&Geant4UIManager::forceExit),0); +} + +/// Force exiting this process without calling atexit handlers +void Geant4UIManager::forceExit() { + std::_Exit(0); +} + /// Start visualization G4VisManager* Geant4UIManager::startVis() { // Initialize visualization diff --git a/DDG4/src/Geant4UIMessenger.cpp b/DDG4/src/Geant4UIMessenger.cpp index 69aecc24056c0c128c22b752f5602b553e7227d8..b875f128b2fbbb3f9c4bc40b61b267f24e40da53 100644 --- a/DDG4/src/Geant4UIMessenger.cpp +++ b/DDG4/src/Geant4UIMessenger.cpp @@ -59,10 +59,21 @@ Geant4UIMessenger::~Geant4UIMessenger() { } /// Add a new callback structure -void Geant4UIMessenger::addCall(const std::string& name, const std::string& description, const Callback& cb) { - G4UIcommand* cmd = new G4UIcmdWithoutParameter((m_path + name).c_str(), this); - cmd->SetGuidance(description.c_str()); - m_actionCmd[cmd] = cb; +void Geant4UIMessenger::addCall(const std::string& name, const std::string& description, const Callback& cb, size_t npar) { + if ( 0 == npar ) { + G4UIcommand* cmd = new G4UIcmdWithoutParameter((m_path + name).c_str(), this); + cmd->SetGuidance(description.c_str()); + m_actionCmd[cmd] = cb; + } + else if ( 1 == npar ) { + G4UIcmdWithAString* cmd = new G4UIcmdWithAString((m_path + name).c_str(), this); + cmd->SetParameterName("p1", true); + cmd->SetGuidance(description.c_str()); + m_actionCmd[cmd] = cb; + } + else { + except("Geant4UIMessenger","+++ Currently only callbacks with one argument are handled! [Contact developers if more are required]"); + } } /// Export all properties to the Geant4 UI @@ -120,7 +131,8 @@ void Geant4UIMessenger::SetNewValue(G4UIcommand *c, G4String v) { Actions::iterator j = m_actionCmd.find(c); if (j != m_actionCmd.end()) { try { - (*j).second.execute(0); + const void* args[] = {v.c_str(), 0}; + (*j).second.execute(args); } catch(const exception& e) { printout(INFO, "Geant4UI", "+++ %s> Exception: Failed to exec action '%s' [%s].", diff --git a/DDParsers/include/Evaluator/DD4hepUnits.h b/DDParsers/include/Evaluator/DD4hepUnits.h index b64e6f48b7ee2d94e5a488238f4ee604d30529f9..6d5a32d4b97158334728bfdf3ddeb2eedcc2bc12 100644 --- a/DDParsers/include/Evaluator/DD4hepUnits.h +++ b/DDParsers/include/Evaluator/DD4hepUnits.h @@ -204,9 +204,9 @@ namespace dd4hep { // // Pressure [E][L^-3] // -#define pascal hep_pascal // a trick to avoid warnings +#define pascal hep_pascal // a trick to avoid warnings static constexpr double hep_pascal = newton / m2; // pascal = 6.24150 e+3 * MeV/mm3 - static constexpr double bar = 100000 * pascal; // bar = 6.24150 e+8 * MeV/mm3 + static constexpr double bar = 100000 * pascal; // bar = 6.24150 e+8 * MeV/mm3 static constexpr double atmosphere = 101325 * pascal; // atm = 6.32420 e+8 * MeV/mm3 // @@ -408,12 +408,16 @@ namespace dd4hep { static constexpr double k_Boltzmann = 8.617343e-11 * MeV / kelvin; // + // IUPAC standard temperature and pressure (STP) + // STP uses 273.15 K (0 °C, 32 °F) and (since 1982) 1 bar (100 kPa) and not 1 atm! + static constexpr double STP_Temperature = 273.15 * kelvin; + static constexpr double STP_Pressure = 1. * bar; // + // NTP uses the NIST convention: 20 °C (293.15 K, 68 °F), 1 atm (14.696 psi, 101.325 kPa) + static constexpr double NTP_Temperature = 293.15 * kelvin; + static constexpr double NTP_Pressure = 1. * atmosphere; // - static constexpr double STP_Temperature = 273.15 * kelvin; - static constexpr double STP_Pressure = 1. * atmosphere; - static constexpr double kGasThreshold = 10. * mg / cm3; - + static constexpr double kGasThreshold = 10. * mg / cm3; // // // diff --git a/examples/ClientTests/CMakeLists.txt b/examples/ClientTests/CMakeLists.txt index 0c8393fc81400d2d278ec8a7061fc96d0c19560f..26f220c4a93171bcbb4326c2c054eaa800614a2f 100644 --- a/examples/ClientTests/CMakeLists.txt +++ b/examples/ClientTests/CMakeLists.txt @@ -208,6 +208,14 @@ dd4hep_add_test_reg( ClientTests_Check_VolumeDivisionTest -plugin DD4hep_VolumeDump REGEX_PASS "Checked 8 physical volume placements") # +# Test Setting temperature and pressure to material +dd4hep_add_test_reg( ClientTests_Check_Temp_Pressure_Air + COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_ClientTests.sh" + EXEC_ARGS geoPluginRun + -destroy -input file:${ClientTestsEx_INSTALL}/compact/Check_Air.xml + -plugin DD4hep_MaterialTable -name Vacuum + REGEX_PASS "Temp=111 \\[Kelvin\\] Pressure=6.2415e-07 \\[pascal\\] state=Undefined" + REGEX_FAIL "Exception;EXCEPTION;ERROR;Error;FATAL" ) # # only if root version > 6.19: MaterialTester # @@ -296,6 +304,15 @@ endforeach() # if (DD4HEP_USE_GEANT4) # + # Test Setting temperature and pressure to material + dd4hep_add_test_reg( ClientTests_sim_Check_Temp_Pressure_Air + COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_ClientTests.sh" + EXEC_ARGS python ${ClientTestsEx_INSTALL}/scripts/Check_Air.py + -geometry ${ClientTestsEx_INSTALL}/compact/Check_Air.xml + REQUIRES DDG4 Geant4 + REGEX_PASS "Imean: 85.538 eV temperature: 333.33 K pressure: 2.22 atm" + REGEX_FAIL "Exception;EXCEPTION;ERROR;Error;FATAL" ) + # # Geant4 test if production cuts are processed dd4hep_add_test_reg( ClientTests_sim_MiniTel_prod_cuts COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_ClientTests.sh" diff --git a/examples/ClientTests/compact/Check_Air.xml b/examples/ClientTests/compact/Check_Air.xml new file mode 100644 index 0000000000000000000000000000000000000000..9da3504f17e9be00d11b2fef23db06cda7b15b83 --- /dev/null +++ b/examples/ClientTests/compact/Check_Air.xml @@ -0,0 +1,66 @@ +<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0" + xmlns:xs="http://www.w3.org/2001/XMLSchema" + xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> + + <info name= "ShapeCheck" + title= "Checking individual shapes in D4hep by comparing the mesh vertices" + author= "Markus Frank" + url= "http://www.cern.ch/lhcb" + status= "development" + version="1.0"> + <comment>Shape tester</comment> + </info> + <steer open="false" close="false"/> + + <debug> + <type name="materials" value="1"/> + </debug> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + </includes> + + <define> + <constant name="world_side" value="300*cm"/> + <constant name="world_x" value="world_side"/> + <constant name="world_y" value="world_side"/> + <constant name="world_z" value="world_side"/> + </define> + + <std_conditions type="STP"> + <T unit="kelvin" value="333.33"/> + <P unit="atmosphere" value="2.22"/> + </std_conditions> + + <display> + <vis name="Box_vis" alpha="1.0" r="1" g="0" b="0" showDaughters="true" visible="true"/> + </display> + + <materials> + <material name="Air"> + <D type="density" unit="g/cm3" value="0.001214"/> + <fraction n="0.7494" ref="N"/> + <fraction n="0.2369" ref="O"/> + <fraction n="0.0129" ref="Ar"/> + <fraction n="0.0008" ref="H"/> + </material> + + <!-- We model vakuum just as very thin air --> + <material name="Vacuum"> + <D type="density" unit="g/cm3" value="0.0000000001" /> + <T unit="kelvin" value="111.11"/> + <P unit="pascal" value="1e-10"/> + <composite n="1" ref="Air"/> + </material> + </materials> + + <detectors> + <detector id="3" name="B1" type="DD4hep_BoxSegment" vis="Box_vis"> + <comment>Vertical box</comment> + <material name="Vacuum"/> + <box x="10" y="20" z="30"/> + <position x="-10" y="30" z="10"/> + <rotation x="0" y="0" z="0"/> + </detector> + </detectors> +</lccdd> diff --git a/examples/ClientTests/scripts/Check_Air.py b/examples/ClientTests/scripts/Check_Air.py new file mode 100644 index 0000000000000000000000000000000000000000..dac842bae35be4a2c4b2b5d93990646bf821cc6d --- /dev/null +++ b/examples/ClientTests/scripts/Check_Air.py @@ -0,0 +1,70 @@ +""" + dd4hep example setup using the python configuration + + \author M.Frank + \version 1.0 + +""" +from __future__ import absolute_import, unicode_literals +import logging +import sys + + +logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO) + + +def help(): + logging.info("Check_shape.py -option [-option] ") + logging.info(" -geometry <geometry file> Geometry file ") + logging.info(" -vis Enable visualization ") + logging.info(" -batch Batch execution ") + + +def run(): + geo = None + vis = False + batch = False + for i in xrange(len(sys.argv)): + c = sys.argv[i].upper() + if c.find('BATCH') < 2 and c.find('BATCH') >= 0: + batch = True + if c[:4] == '-GEO': + geo = sys.argv[i + 1] + if c[:4] == '-VIS': + vis = True + + if not geo: + help() + sys.exit(1) + + import DDG4 + kernel = DDG4.Kernel() + # Configure UI + geant4 = DDG4.Geant4(kernel, tracker='Geant4TrackerCombineAction') + if batch: + kernel.UI = '' + ui = geant4.setupCshUI(ui=None, vis=vis) + else: + ui = geant4.setupCshUI(vis=vis) + kernel.loadGeometry(geo) + # Configure field + geant4.setupTrackingField(prt=True) + # Now build the physics list: + geant4.setupPhysics('') + kernel.physicsList().enableUI() + DDG4.setPrintLevel(DDG4.OutputLevel.DEBUG) + # + ui.Commands = [ + '/ddg4/ConstructGeometry/printVolume /world_volume_1' + ,'/ddg4/ConstructGeometry/printMaterial Air' + ,'/ddg4/ConstructGeometry/printMaterial Vacuum' + ,'/ddg4/UI/exit' + ] + kernel.NumEvents = 0 + kernel.configure() + kernel.initialize() + kernel.run() + kernel.terminate() + +if __name__ == "__main__": + run() diff --git a/examples/ClientTests/scripts/Check_shape.py b/examples/ClientTests/scripts/Check_shape.py index a8c46415bcdc20d4b47f4c02375ae3b6d4833dc4..55bad194eba1676b0f21c870b0a70ddb50e2f9ac 100644 --- a/examples/ClientTests/scripts/Check_shape.py +++ b/examples/ClientTests/scripts/Check_shape.py @@ -41,10 +41,11 @@ def run(): kernel = DDG4.Kernel() # Configure UI geant4 = DDG4.Geant4(kernel, tracker='Geant4TrackerCombineAction') - ui = geant4.setupCshUI(ui=None, vis=vis) if batch: kernel.UI = '' - + ui = geant4.setupCshUI(ui=None, vis=vis) + else: + ui = geant4.setupCshUI(vis=vis) kernel.loadGeometry(geo) # Configure field geant4.setupTrackingField(prt=True) @@ -53,12 +54,10 @@ def run(): kernel.physicsList().enableUI() DDG4.setPrintLevel(DDG4.OutputLevel.DEBUG) # - # '/ddg4/ConstructGeometry/DumpGDML test.gdml', - # '/ddg4/ConstructGeometry/writeGDML', + # '/ddg4/ConstructGeometry/writeGDML test.gdml', ui.Commands = [ - '/ddg4/ConstructGeometry/VolumePath /world_volume_1/Shape_Test_0/Shape_Test_vol_0_0', - '/ddg4/ConstructGeometry/printVolume', - 'exit' + ,'/ddg4/ConstructGeometry/printVolume /world_volume_1/Shape_Test_0/Shape_Test_vol_0_0' + ,'exit' ] kernel.NumEvents = 0 kernel.configure() @@ -66,6 +65,5 @@ def run(): kernel.run() kernel.terminate() - if __name__ == "__main__": run()